Student Support Forum: 'Bifurcation analysis of differential system' topicStudent Support Forum > General > Archives > "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