Mathematica 9 is now available
Student Support Forum
-----
Student Support Forum: 'Bifurcation analysis of differential system' topicStudent Support Forum > General > "Bifurcation analysis of differential system"

Help | Reply To Topic
Author Comment/Response
Derek
06/28/08 12:05pm

Hello,

I've been trying to generate bifurcation diagrams for a system of differential equations but got stuck in the technicalities.
If anyone can show me how to do it simply then i'd be very grateful. I post my attempt below (from copying an example from a simpler system). if someone could show me where i'm going wrong then that would be great!! I would like to vary a parameter, for example v, and see how the dynamics change over time- is this do-able? I have seen it done for discrete single species systems but cannot find anything for something like this. Many thanks, Derek

system::::==================================
[{If[H[t] < 0, 0, a H[t] - ap P[t] + at T[t] - (H[t] b + bp P[t] + bt T[t])] == H'[t], If[P[t] < 0 || H[t] < 0, 0, P[t] ((L H[t])/(H0 + H[t]) - u) - (H[t] b + bp P[t] + bt T[t]) (P[t]/H[t]) + phi T[t]] == P'[t], If[T[t] < 0 || H[t] < 0, 0,
H[t] (theta - rho) - (H[t] b + bp P[t] + bt T[t]) (T[t]/H[t]) - v P[t]] == T'[t],

==============================================
attempt:::
a = 3; b = 1.31; u = .1; H0 = 10; ap = .3; L = 4; bp = 1.2;

at = .2; bt = 1.2; phi = .8; theta = .75; rho = .005; v = .213; tstop
= 50;
tmin = 0; tmax = 50; H0 = 5; P0 = 5; T0 = 5;

sol[v] = NDSolve[{If[H[t] < 0, 0,
a H[t] - ap P[t] + at T[t] - (H[t] b + bp P[t] + bt T[t])] ==
H'[t], If[P[t] < 0 || H[t] < 0, 0,
P[t] ((L H[t])/(H0 + H[t]) - u) - (H[t] b + bp P[t] +
bt T[t]) (P[t]/H[t]) + phi T[t]] == P'[t],
If[T[t] < 0 || H[t] < 0, 0,
H[t] (theta - rho) - (H[t] b + bp P[t] + bt T[t]) (T[t]/H[t]) -
v P[t]] == T'[t], H[tmin] == H0, P[tmin] == P0,
T[tmin] == T0}, {H[t], T[t], P[t]}, {t, tmin, tmax}];

Clear[bifurcationplot]

bifurcationplot[{vbegin_, vend_, vstep_}, timesteps_, takelast_] :=
Module[{H, iter, v, thelist = {}}, H[v_] := Evaluate[H[t] /. sol[v]]

iter[v_] := Table[sol[v], {t, 0, timesteps}];

For[v = vbegin, v < vend, v += vstep,
thelist =
Flatten[{thelist, ({v, #1} &) /@ Take[iter[v], -takelast]}, 1]];

ListPlot[thelist, PlotStyle -> PointSize[0.002],
PlotRange -> {Automatic, {0, 2000}}]]


bifurcationplot[{0, 10, 0.3}, 200, 20]

URL: ,
Help | Reply To Topic