Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

How write the 'for' loop in mathematica

  1. Dec 10, 2008 #1
    I have a problem how write the 'for' loop in mathematica:

    my previous code looks like:
    count=0;
    flag=0;
    step=100;

    For i=1:1000
    count=count+1;
    if count>step
    if flag== 0;
    flag=1;
    else
    flag=0;
    end
    count=1;
    end
    if flag==0
    p(i)=g
    else
    p(i)=-g
    end
    end

    I need to rewrite this in mathematica but I do not know how.
    Is it possible to write as a function, module or block?
    thanks
     
  2. jcsd
  3. Dec 10, 2008 #2

    alphysicist

    User Avatar
    Homework Helper

    Re: matlab-mathematica

    Hi heinricht,

    I think something like this would closely follow your program:

    Code (Text):

    count=0
    flag=0
    step=100

    p=Table[0,{i,1,1000}]  (* add a semicolon after this to not display the zero matrix*)


    Do[{
    count=count+1,
    If[ count > step , {If[ flag==0, flag=1, flag=0], count=1} ],
    If[flag==0, p[[i]]=g, p[[i]]=-g]} , {i,1,1000}]

    p   (* this is just displaying the result*)
     

    It looks like your program is creating alternating sequences of +g and -g of 100 terms each. If that is your goal, you can more directly create it by using something like:

    Code (Text):

    pos=Table[g,{i,1,100}]

    neg=Table[-g,{i,1,100}]

    p={}

    Do[ {p = Append[p,pos],
            p=Append[p,neg]},{i,1,5}]

    p=Flatten[p]   (*removing nested brackets *)

    p   (*display result*)
     
    and you can easily shorten this quite a bit (but this form is understandable, I think).



    As an example of a shorter program, you might try,

    Code (Text):

    p=g  Flatten[
          Table[ Append[ Table[1,{i,1,100}],Table[-1,{i,1,100}]],{i,1,5}]]
     
    which does the same as the previous two programs.
     
  4. Dec 11, 2008 #3
    Re: matlab-mathematica

    Thanks a lot of. I have one more question:
    I have equation: NDSolve[m (q[t] - M[t] - g (2 M[t] + 1) k) == M'[t], {t,0,1000}]
    and I need to use p from previous program as a function of q[t]
    I used p[[t]] which take t element from the list p and wrote
    q[t_]=p[[t]]
    is this correct?

    or can I write sth like:

    p[t_]:= Module[{step, avg, amp, freq, per, t1, t0, time},
    amp = 0.1;
    avg = 0.25;
    freq = 10^8;
    per = 1/freq;
    t1 = 100/10^9;
    t0 = 0;
    time = (t1 - t0)/t;
    step = per/time;

    count = 0;
    state = 0;
    Do[{count = count + 1,
    If[count > step, {If[state == 0, state = 1, state = 0],
    count = 1}],
    If[state == 0, p[] = avg + amp/2, p[] = avg - amp/2]}, {i,
    1, t}];
    p[[t]]];

    and then use this function in ODE equation:
    NDSolve[m (p[t] - M[t] - g (2 M[t] + 1) k) == M'[t], {t,0,1000}]
     
    Last edited: Dec 11, 2008
  5. Dec 11, 2008 #4

    alphysicist

    User Avatar
    Homework Helper

    Re: matlab-mathematica

    Don't you need an initial condition for this?

    I don't think this will work. As NDSolve is solving, you need to be able to get values for all t between 0 and 1000. However if you try to access q[1.5], for example, you'll get an error, since p[[t]] requires t to be an integer.




    Well, I'm not really certain what you are doing here, so I'll just make a few comments. You need to call the module something beside p, since p is the list inside the module (right now you have p[t_] on the left hand side of the equation, but also use p as the list called with p[]).

    Also, you need to make sure the the list p is created before this module is called.

    And here again there is the issue with non-integer values of t, which will cause an error with the module the way it is written. For example, let's say that you know that

    p[[25]]=.2
    p[[26]]=.3

    Now when NDSolve gets to t=25.678, what value do you want it to call? There is no p[[256.678]]. Do you want it to be equal to p[[25]] or p[[26]], or something else?

    For example, if you always want the index rounded down, you can change the last line of your module to this:


    p[[ Floor[t] ]]];

    which will just truncate the part of t after the decimal. Then when t=25.678, you would always get the value for p[[25]]. (Instead of Floor, you can use Round or Ceiling to get different behaviors.)



    (Again, I'm not sure of your goal here, but you might want to consider not creating a list p, and instead just creating a function.)
     
  6. Dec 12, 2008 #5
    Re: matlab-mathematica

    my previuos program is below: I used p as a constant-number(0.25)

    constant number: p,g0,ep,gamma,tau,m
    number = 50000; (*number_steps*)

    g10 = g0/(1 + ep F0);
    int= NSolve[{(g10/g0 (2 M0 + 1) - 1) F0 == 0,
    m (p - M0 - (g0 (2 M0 + 1)) F0) == 0} , {M0, F0}];
    x0 = M0 /. int;
    y0 = F0 /. int;
    (*Initial conditions*)
    m00 = x0[[3]]
    f00 = y0[[3]]

    g = g0/(1 + ep F[t]);
    s = NDSolve [{(g/g0 (2 M[t] + 1) - 1) F[t] == F'[t],
    m (p - M[t] - g (2 M[t] + 1) F[t]) == M'[t],
    F[0] == f00, M[0] == m00},
    {F[t], M[t]}, {t, 1, number}] ;
    time = Table[{t/\[Gamma]}, {t, number}];
    solF = Table[F[t] /. s, {t, number}];
    solM = Table[M[t] /. s, {t, number}];
    fo = Flatten[solF];
    mo = Flatten[solM];
    tg = Flatten[time];
    final = Thread[{tg, fo}];
    final2 = Thread[{tg, mo}];
    ListLinePlot[final, PlotRange -> All]
    ListLinePlot[final2, PlotRange -> All]

    and this program is working fine and now I need use p in NDSolve as a function of p[t] using formula
    amp = 0.1;
    avg = 0.25;
    f = 10^8;
    per = 1/f;
    t1 = 100/10^9;
    t0 = 0;
    time = (t1 - t0)/number;
    step = per/time;

    count = 0;
    state = 0;
    Do[{count = count + 1,
    If[count > step, {If[state == 0, state = 1, state = 0],
    count = 1}],
    If[state == 0, p[] = avg + amp/2, p[] = avg - amp/2]}, {i,
    1, t}]

    beacuse p is also in NSolve where I have initial conditions from physical point of view I can use p0=0.01 so I do not have to calculate exactly value.

    Can I change also in NDSolve time step=1? I gues that now the time step is 0.001 that I have funny thing p[[25.687]] if not could you help me with these function?
     
  7. Dec 12, 2008 #6

    alphysicist

    User Avatar
    Homework Helper

    Re: matlab-mathematica



    I don't think you would want to do that; I believe mathematica needs to adjust the step size itself so as to get a good result.

    The function needs to be a step function, right? There are many ways to do this, but I think a good way that never returns a zero would be something like (after you have already calculated step, avg, and amp):

    Code (Text):

    p[x_]:= If [ EvenQ [ Quotient [ x , step ] ] , avg + amp/2 , avg - amp/2 ]
     
    Once that is defined, if I calculate:

    Code (Text):

    step=3;

    avg=12;

    amp=6;

    Plot[p[x],{x,0,10} , PlotRange->{0,20}]
     
    I get this:

    http://img154.imageshack.us/img154/4494/squarewavedq2.jpg

    So as long as you calculate step, avg, and amp before you run NDSolve I believe it will work.




    (I also think you could use your module, but change the last line of the module to

    p[[ Ceiling[t] ]] ];

    to get the same behavior.)
     
  8. Dec 15, 2008 #7
    Re: matlab-mathematica

    I used your function-see below- in my program. The step is fine (picture1) but when I used these function in NDSolve it does not work (picture 2): see the attachment


    g0 = 1;
    ep = 1/10^2;
    Gamma = 0.98/10^-12;
    Tau = 0.816/10^9;
    m = 1/(Gamma*Tau);
    number = 100 000;

    g10 = g0/(1 + ep F0);
    intensity = NSolve[{(g10/g0 (2 M0 + 1) - 1) F0 == 0,
    m (avg - M0 - (g0 (2 M0 + 1)) F0) == 0} , {M0, F0}];
    x0 = M0 /. intensity;
    y0 = F0 /. intensity;
    (*Initial conditions*)
    Print["I0: ", m00 = x0[[3]]]
    Print["N0: ", f00 = y0[[3]]]

    amp = 0.1;
    avg = 0.25;
    freq = 10^8; (*Hz*)
    per = 1/freq; (*period of modulated current*)
    step = per/10^-12;

    p[t_] := If[EvenQ[Quotient[t, step]], avg + amp/2, avg - amp/2];
    Plot[p[t], {t, 0, 50000}, PlotRange -> {0, 0.4}]

    g = g0/(1 + ep F[t]);
    s = NDSolve [{(g/g0 (2 M[t] + 1) - 1) F[t] == F'[t],

    m (p[t] - M[t] - g (2 M[t] + 1) F[t]) == M'[t],
    F[0] == f00, M[0] == m00},
    {F[t], M[t]}, {t, 1, number}] ;

    Plot[F[t] /. s, {t, 0, number}, PlotRange -> All]
    Plot[M[t] /. s, {t, 0, number}, PlotRange -> All]

    I also used previous function: p[t_] := If[t > step, avg + amp/2, avg - amp/2];
    but the result is only one step: see the picture 3 and 4
     

    Attached Files:

    Last edited: Dec 15, 2008
  9. Dec 15, 2008 #8

    alphysicist

    User Avatar
    Homework Helper

    Re: matlab-mathematica

    What kind of result were you expecting? If you plot from 0 to 1000 (instead of from 0 to number) I believe you can see there are smooth oscillations; should there be something different?


    (I noticed a few minor things, some of which I had to change to get your program to work for me. There was an extra space in the assignment for number (it was 100 000 instead of 100000), and I had to make the G in Gamma lowercase. Were those just typos?

    Also, you had NDSolve go from 1 to number, and plotted the results from 0 to number. Was that on purpose?)
     
  10. Dec 16, 2008 #9
    Re: matlab-mathematica

    I have done these what I need exactly. I used simple function
    p[t_] := If[t > step, If[t > 2*step, avg + amp/2, avg - amp/2], avg + amp/2];
    and the result are below(data1,jpg)
    Now I need to multiply steps only.

    I tried used p[t_] as module and obtain the same result but it does not work. Can you tell me where might be the problem? Thanks.

    p[t_] := Module[{step, avg, amp, freq, per, t1, t0, time, number,q},
    amp = 0.1;
    avg = 0.25;
    freq = 10^8;
    per = 1/freq;
    t1 = 100/10^9;
    t0 = 0;
    number = 100000;
    time = (t1 - t0)/number;
    step = per/time;
    count = 0;
    state = 0;
    Do[{count = count + 1,
    If[count > step, {If[state == 0, state = 1, state = 0],
    count = 1}],
    If[state == 0, q[[t]] = avg + amp/2, q[[t]] = avg - amp/2]}, {i,
    1, t}];
    q[[Ceiling[t]]]
    ];
     

    Attached Files:

    Last edited: Dec 16, 2008
  11. Dec 18, 2008 #10

    alphysicist

    User Avatar
    Homework Helper

    Re: matlab-mathematica


    What is picture 2 from your last post? Is it a plot of the module p[t]? Or is it a plot of the results after you run NDSolve?


    If it is just a plot of the module p[t], how did you plot it? Did you use

    Plot[p[t],{t,0,50000}]

    or something different? When I plot that for the module p, I get an identical result to your picture 1 (although there were a few errors about the list q that were reported, so I had to modify). It does not look to me like you need to build a list q inside the module, since only its last result is used.
     
  12. Jan 8, 2009 #11
    Re: matlab-mathematica

    Hi in new year :),

    the second picture is a result of NDSolve.
    I did not plot the module p[t].
     
  13. Jan 8, 2009 #12
    Re: matlab-mathematica

    Now I have a problem with number of steps in NDSolve.
    My program:
    avg = 0.5;
    g10 = 1;
    dg = 10^-3;
    g20 = g10 - dg;
    alpha = 2.6;
    i = Sqrt[-1];
    ep = 1/10^2;
    Gamma = 0.98/10^-12;
    Tau = 0.816/10^9;
    m = 1/(Gamma*Tau);
    beta = 2/3;
    number = 30000;

    amp = 0.2;
    freq = 10^8;
    tran = 10^-9;
    per = 1/freq;
    step = per/10^-12;
    steptran = tran/10^-12;
    phigh = avg + amp/2;
    plow = avg - amp/2;
    slope = (phigh - plow)/steptran;

    c[t_] := If[t > step,
    If[t > 2*step,
    If[t < 3*step - steptran, phigh, phigh - (t - (3*step - steptran))slope ],
    If[t < 2*step - steptran, plow, plow + (t - (2*step - steptran)) slope ] ],
    If[t < step - steptran, phigh, phigh - (t - (1*step - steptran)) slope ] ];

    Plot[c[t], {t, 0, number}, PlotRange -> {0.3, 0.65}]

    Clear[M, F1, F2, t];
    Mu = 1.0163;
    Sigma = 0.75;
    g17 = g10/(1 + ep F1[t] + ep beta F2[t]);
    g27 = g20/(1 + ep beta F1[t] + ep F2[t]);
    g1 = g17 1/(
    1 - Sigma (1 - Mu)) (1 - Sigma] (2 M[t] + 1 - Mu));
    g2 = g27 1/(
    1 - Sigma (1 - Mu)) (1 + Sigma (2 M[t] + 1 - Mu));

    s = NDSolve [{(g1/g10 (2 M[t] + 1) - 1) F1[t] == F1'[t],
    (g2/g10 (2 M[t] + 1) - 1) F2[t] == F2'[t],
    m (c[t] - M[t] - (2 M[t] + 1) (g1 F1[t]) + g2 F2[t]) == M'[t],
    F2[0] == 0.01, F1[0] == 0.02, M[0] == 0.03},
    {F1[t], F2[t], M[t]}, {t, 0, number}] ;
    Plot[F1[t] /. s, {t, 0, number}, PlotRange -> All]
    Plot[F2[t] /. s, {t, 0, number}, PlotRange -> All]

    When I run I got message:
    NDSolve::mxst: Maximum number of 10000 steps reached at the point t == \
    1942.86059216402`. >>

    is this connected with MaxSteps in NDSolve? or sth different
     
    Last edited: Jan 8, 2009
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: How write the 'for' loop in mathematica
  1. Mathematica looping (Replies: 3)

Loading...