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}]