Re: NDSolve, Do loop, and Plot

*To*: mathgroup at smc.vnet.net*Subject*: [mg77439] Re: NDSolve, Do loop, and Plot*From*: DBK <boydkramer at gmail.com>*Date*: Fri, 8 Jun 2007 05:42:06 -0400 (EDT)*References*: <f48eqh$df9$1@smc.vnet.net>

On Jun 7, 4:15 am, chuck009 <dmili... at comcast.com> wrote: > Your code worked fine on my machine (5.2). Also, to get it to work on a list of w values, I just defined your code as a pure function and then used the Map construct applied to a list of w values (the {2,3,4} list at the end): > > Map[mysol[#] = NDSolve[{x'[ > t] == -y[t] - x[t]^2, y'[t] == # x[t] - y[t]^3, x[0] == y[0] == 1}, {x, > y}, {t, 20}]; > Plot[{Evaluate[x[t] /. mysol[#]], > Evaluate[y[t] /. mysol[#]]}, {t, 0, 20}] &, {2, 3, 4}] Thanks all for the many alternative ways of doing this and for your patience. I am still learning. I've got three follow up questions. 1. Each code that was sent plots x[t] and y[t] on one graph for one value of a changing parameter w (w has three values). This gives me a total of three graphs. How do I plot x[t] three times on one graph for each value of the parameter w and the same for y[t]. This would give me 2 graphs in total, each with 3 lines. 2. The code I originally submitted was a simple example meant to help me understand various approaches to my actual problem with consists of 6 differential equations. I've tried each approach with the more complicated system of equations and am not having much success. Code below using one suggestion. What am I doing wrong? 3. I also would like to evaluate the stability of the Jacobian matrix using numerical methods. Any thoughts on the best way to approach this? ClearAll["Global`*"]; $TextStyle={FontFamily->"Arial",FontSize->12,FontSlant->"Italic"}; time0; timeend=time-100; herbstart=2100; piscstart=1400; algaestart=0.3; coralstart=0.3; herbharveststart=0.5; pischarveststart=0.5; aha=0.06349206; ahp=0.00001020; ahm=0.00001905; ahh=0; aph=0.0000068; app=0.007; apm=0.00004286; ra=0.28; ka=0.80; aah=0.00003333; aac=0.80; rc=0.15; kc=0.70; aca=0.60; pop=50; lag=1.0; effortP=2.0; effortH=2.0; PriceP=2.0; PriceH=2.0; CostP=1; CostH=1; slope=7.0; ha=0.30; harvestH=effortH ahm h[t]; harvestP=effortP apm p[t]; profitH=harvestH(PriceH-CostH/(h[t]+1)); profitP=harvestP (PriceP-CostP/(p[t]+1)); profitMEAN=sH[t] profitH+sP[t] profitP; equations={a'[t]==(ra a[t] (ka-a[t]-aac c[t]))/ka-aah h[t] \ a[t],c'[t]==(rc c[t] (kc-c[t]-(aca \ a[t]^slope)/(a[t]^slope+ha^slope)))/kc,h'[t]==aha h[t] \ a[t]-ahh h[t]-ahp h[t] p[t]-pop sH[t] ahm effortH \ h[t],p'[t]==aph p[t] h[t]-app p[t]-pop sP[t] apm effortP \ p[t],sH'[t]==(lag sH[t] (profitH-profitMEAN))/profitMEAN,sP'[t]==(lag sP[t] (profitP-profitMEAN))/profitMEAN}; mysol[lag_?NumericQ]:=NDSolve[{equations,sH[0]==herbharveststart,sP[\ 0]==pischarveststart,h[0]==herbstart,p[0]==piscstart,a[0]==algaestart, \ c[0]==coralstart},{h,p,sH,sP,a,c},{t,0,time}][[1]] Do[Print[Plot[Evaluate[{h[t],p[t],sH[t],sP[t],a[t],c[t]}/.mysol[lag]], \ {t,0,200}]],{lag,1,3}]