Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

NDSolve, Loop, Table, Plot

146 views
Skip to first unread message

Solver

unread,
Jun 13, 2007, 7:54:49 AM6/13/07
to
Hi,

I have a system of 3 differential equations and I would like to
numerically solve them for different values of one parameter (e.g., p1
= 0.5, 0.6, 0.7, 0.8, etc). Then, I want to 1) make a table for the
results of B0[200], B1[200], and B2[200] for all values of p1 and 2)
make plots of B0, B1, B2 through time for the different values of p1.
Below is the code I tried. The table seems to work sometimes but when
I run the code multiple times, the output is the same for all values
of p1 (i.e., {{3.99315, 3.99315, 3.99315, 3.99315, 3.99315, 3.99315},
{
0.0813405, 0.0813405, 0.0813405, 0.0813405, 0.0813405,
0.0813405}, {0.021349, 0.021349, 0.021349, 0.021349, 0.021349,
0.021349}}. I don't understand why?
The Plot function gives me multiple plots but they are always the
same, so it is not drawing plots for the different values of p1.

J = 1;
r0 = 0.2;
r1 = 0.2;
r2 = 0.2;
p1 = 0.1;

dB0 = J - r0*B0[t] - B0[t]*B1[t];
dB1 = p1*B0[t]*B1[t] - r1*B1[t] - B1[t]*B2[t];
dB2 = p2*B1[t]*B2[t] - r2*B2[t];

A = Flatten[Map[({B0[200], B1[200], B2[200]} /. # ) &,
sol = Table[NDSolve[{B0'[t] == dB0, B1'[t] == dB1, B2'[t] == dB2,
B0[0] == 1, B1[0] == 1, B2[0] == 1}, {B0, B1, B2},
{t, 0, 200}], {p1, 0.5, 1, 0.1}]], 1] // Transpose
Export["A.csv", A];
Do[Plot[Evaluate[{B0[t], B1[t], B2[t]} /. sol, {t, 0, 200}, PlotStyle -
>
{{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[0, 0, 1]}}]],
{p1, 0.5, 1, 0.1}]

Thanks!


Jens-Peer Kuska

unread,
Jun 14, 2007, 5:20:04 AM6/14/07
to
Hi,

you can make a function, that take the parameters of
your ode as arguments, something like


sol[omega_] :=
Module[{f, y, x},
f = y[x] /.
NDSolve[{y''[x] + omega^2*y[x] == 0, y[0] == 0, y'[0] == 1},
y[x], {x, 0, 2 Pi}][[1]];
Table[{x, f}, {x, 0., 2 Pi, 2 Pi/128}]
]

and

ListPlot[{sol[1], sol[2], sol[3]}]

will do that or you can define Rule[]s for the parameter

NDSolve[{y''[x] + omega^2*y[x] == 0, y[0] == 0, y'[0] == 1},
y[x], {x, 0, 2 Pi}] /. {{omega -> 1}, {omega -> 2}, {omega -> 3}}

Regards
Jens

dh

unread,
Jun 14, 2007, 6:06:46 AM6/14/07
to

Hi Solver

the variable p2 does not have a value. Presumably you wanted to write

p2=0.1 instead of p1=0.1.

If you are using version 6 you must additionally wrap the Plot statement

in a Print statement.

hope this helps, Daniel

rob

unread,
Jun 15, 2007, 4:25:28 AM6/15/07
to
Sir, I tried your code but the ListPlot does not work. I'm
using V5.1.
0 new messages