How write the 'for' loop in mathematica

[heinricht]

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

 

[alphysicist]

Hi heinricht,

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:

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:

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,

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.

 

[heinricht]

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

 

[alphysicist]

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

 

[heinricht]

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?

 

[alphysicist]

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

p[x_]:= If [ EvenQ [ Quotient [ x , step ] ] , avg + amp/2 , avg - amp/2 ]

Once that is defined, if I calculate:

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

 

[heinricht]

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

 

[alphysicist]

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

 

[heinricht]

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

 

[alphysicist]

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.

 

[heinricht]

Hi in new year :),

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

 

[heinricht]

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