Mathematica How write the 'for' loop in mathematica

  • Thread starter Thread starter heinricht
  • Start date Start date
  • Tags Tags
    Loop Mathematica
Click For Summary
The discussion focuses on converting a 'for' loop from a previous programming language into Mathematica. The original code alternates values of a variable based on a counter and a flag, and users suggest using Mathematica's `Do` or `Table` functions for this purpose. A user seeks clarification on using a list as a function in an ODE solved by `NDSolve`, raising concerns about accessing non-integer indices. Suggestions include defining a step function and ensuring proper handling of time steps in the context of numerical solutions. The conversation emphasizes the importance of correctly implementing functions and modules in Mathematica for effective computation.
heinricht
Messages
7
Reaction score
0
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
 
Physics news on Phys.org


Hi heinricht,

heinricht said:
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

I think something like this would closely follow your program:

Code: 
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: 
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: 
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.
 


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:


heinricht said:
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}]

Don't you need an initial condition for this?

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]]

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.

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}]



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.)
 


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?
 


heinricht said:
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 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.

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?

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: 
p[x_]:= If [ EvenQ [ Quotient [ x , step ] ] , avg + amp/2 , avg - amp/2 ]

Once that is defined, if I calculate:

Code: 
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.)
 
Last edited by a moderator:


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
 

Attachments

  • data.JPG
    data.JPG
    15 KB · Views: 593
Last edited:


heinricht said:
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

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?)
 


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]]]
];
 

Attachments

  • data1.JPG
    data1.JPG
    8.5 KB · Views: 591
Last edited:
  • #10


heinricht said:
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]]]
];


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.
 
  • #11


Hi in new year :),

the second picture is a result of NDSolve.
I did not plot the module p[t].
 
  • #12


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:

Similar threads

Mathematica Bar legend for a different function in mathematica
  • · Replies 1 ·
Replies
1
Views
2K
Mathematica Problems solving this differential equation for a Pendulum with Mathematica
  • · Replies 4 ·
Replies
4
Views
4K
MATLAB Conversion from Mathematica to Matlab
  • · Replies 1 ·
Replies
1
Views
2K
Mathematica How to check a particular solution of System of Linear ODEs?
  • · Replies 1 ·
Replies
1
Views
3K
Mathematica Plot gradient vector in ContourPlot
  • · Replies 3 ·
Replies
3
Views
3K
Mathematica Solve complex equation analytically
  • · Replies 19 ·
Replies
19
Views
2K
Mathematica Implementation of tessarines in Mathematica
  • · Replies 1 ·
Replies
1
Views
2K
Mathematica How to simplify expressions in Mathematica with log and exp in them?
  • · Replies 6 ·
Replies
6
Views
4K
Mathematica Numerical integration over a Green's function
  • · Replies 13 ·
Replies
13
Views
2K
Mathematica [Mathematica] How to use for loop to process a list of data
  • · Replies 6 ·
Replies
6
Views
4K