(* * predprey_bif.mma - predator prey bifurcation analysis * * RMM, 8 Jul 06 *) <1.6, d->0.56, b->0.6, k->125, a->3.2, c->50} bifpar = {r->1.6, d->0.56, b->0.6, k->125} App = FullSimplify[ Jac[Fpp[H, L], {H, L}] /. eqpts[[1]] ]; eigpp = FullSimplify[ Eigenvalues[App] ]; (* Create a function to compute the real part of the largest eigenvalue *) maxeig[aval_, cval_] := Max @@ Thread[Re[eigpp /. {a->aval, c->cval} /. bifpar]]; (* Compute equilibrium points/stability as a function of a *) eqval1 = {}; eqval2 = {}; eqval1 = {}; stab1 = {}; stab2 = {}; stab1 = {}; avec = {}; mxval1 = {}; mnval1 = {}; cmin={}; cmax={}; For[aval = 1.35; ind=1, aval < 10, aval += 0.05; ind++, (* Compute the equilibrium points for these parameter values *) AppendTo[eqval1, H /. eqpts[[1]] /. a->aval /. nompar]; AppendTo[eqval2, H /. eqpts[[2]] /. a->aval /. nompar]; AppendTo[eqval3, H /. eqpts[[3]] /. a->aval /. nompar]; (* Determine the stability of the equilibrium point *) AppendTo[stab1, Max @@ Thread[Re[eigpp /. a->aval /. nompar]]]; (* Simulate the system and determine excursions *) L0 = L /. eqpts[[1]] /. a->aval /. nompar; H0 = H /. eqpts[[1]] /. a->aval /. nompar; soln = NDSolve[ { H'[t] == Fpp[H[t], L[t]][[1]], L'[t] == Fpp[H[t], L[t]][[2]], H[0] == H0 + 5, L[0] == L0 + 5 } /. a->aval /. nompar, {H, L}, {t, 0, 250} ]; tvec = Array[#&, 500]/10 + 200; AppendTo[mxval1, Max @@ (MapThread[H, {tvec}] /. soln)]; AppendTo[mnval1, Min @@ (MapThread[H, {tvec}] /. soln)]; (* Figure out some stability boundaries *) AppendTo[cmin, c /. FindRoot[maxeig[aval, c] == 0, {c, 10, 10, 110}]]; AppendTo[cmax, c /. FindRoot[maxeig[aval, c] == 0, {c, 150, 50, 200}]]; AppendTo[avec, aval]; ]; (* Generate a bifurcation diagram *) eqplot = ListPlot[Transpose[{avec, eqval1}], PlotJoined->True, PlotStyle->Thickness[0.01]]; mxplot = ListPlot[Transpose[{avec, mxval1}], PlotJoined->True, PlotStyle->Dashing[{0.025, 0.025}]]; mnplot = ListPlot[Transpose[{avec, mnval1}], PlotJoined->True, PlotStyle->Dashing[{0.025, 0.025}]]; Show[{eqplot, mxplot, mnplot}]; (* Generate simulations corresponding to a couple of different parameters *) L0 = L /. eqpts[[1]] /. a->1 /. nompar; H0 = H /. eqpts[[1]] /. a->1 /. nompar; soln2 = NDSolve[ { H'[t] == Fpp[H[t], L[t]][[1]], L'[t] == Fpp[H[t], L[t]][[2]], H[0] == H0 + 0.5, L[0] == L0 - 0.5 } /. a->1 /. nompar, {H, L}, {t, 0, 100} ]; par2 = Plot[{H[t] /. soln2, L[t] /. soln2}, {t, 0, 100}]; L0 = L /. eqpts[[1]] /. a->3 /. nompar; H0 = H /. eqpts[[1]] /. a->3 /. nompar; soln3 = NDSolve[ { H'[t] == Fpp[H[t], L[t]][[1]], L'[t] == Fpp[H[t], L[t]][[2]], H[0] == H0 + 0.5, L[0] == L0 + -0.5 } /. a->3 /. nompar, {H, L}, {t, 0, 100} ]; par3 = Plot[{H[t] /. soln3, L[t] /. soln3}, {t, 0, 100}]; L0 = L /. eqpts[[1]] /. a->7.5 /. nompar; H0 = H /. eqpts[[1]] /. a->7.5 /. nompar; soln1 = NDSolve[ { H'[t] == Fpp[H[t], L[t]][[1]], L'[t] == Fpp[H[t], L[t]][[2]], H[0] == H0 + 0.5, L[0] == L0 - 0.5 } /. a->7.5 /. nompar, {H, L}, {t, 0, 500} ]; par1 = Plot[{H[t] /. soln1, L[t] /. soln1}, {t, 0, 500}]; (* Plot stable versus unstable region *) ContourPlot[maxeig[x, y], {x, 1.35, 10}, {y, 0, 200}, Contours->{0}, PlotPoints->100]; (* Create a file with the relevant data as well *) Export["predprey.dat", Transpose[{avec, eqval1, stab1, mnval1, mxval1, cmin, cmax}]];