Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: NDSolve: avoiding negative values?

• To: mathgroup at smc.vnet.net
• Subject: [mg45191] Re: [mg45162] NDSolve: avoiding negative values?
• From: Selwyn Hollis <sh2.7183 at misspelled.erthlink.net>
• Date: Sat, 20 Dec 2003 05:55:56 -0500 (EST)
• References: <200312191157.GAA29369@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Hi Maanu,

The negative values are coming as a result of numerical error
in a solution that's very close to zero.

Here's an improvement. It uses WorkingPrecision->20 and
MaxStepSize->.1 within NDSolve. Also note that I changed real
constants like 10.7 to exact forms like 107/10.

sol := NDSolve[{
V'[t] == U[t]*(1-V[t]/Kv[t]) - (A*V[t]*H[t])/(V[t]+α*M[t]+B),
M'[t] == u[t]*M[t]*(1-M[t]/Km) - (α*A*M[t]*H[t])/(V[t]+α*M[t]+B),
H'[t] == R*H[t]*((A*(V[t] + α*M[t]))/(V[t] + α*M[t] + B) - G[t]),
V[0] == 1000, M[0] == 1000, H[0] == 2}, {V, M, H}, {t, 0, 100},
MaxSteps->Infinity, WorkingPrecision->20, MaxStepSize->.1]

A = 15;
α = 1/2;
B = 70;
Km = 1000;
R = 107/(10*A);

m[t_] := t - 12*IntegerPart[t/12];
U[t_] := 0 /; m[t] <= 10;
U[t_] := 10000/2 /; m[t] > 10;
u[t_] := 0 /; m[t] <= 10;
u[t_] := 12/2 /; m[t] > 10;
G[t_] := (63*A)/100 /; m[t] <= 10;
G[t_] := (44*A)/100 /; m[t] > 10;
Kv[t_] := 100 /; m[t] <= 10;
Kv[t_] := 1000 /; m[t] > 10;

sol1 = sol; Plot[{V[t]/.sol1, M[t]/.sol1, H[t]/.sol1}, {t,0,100},
    PlotRange -> All, PlotPoints -> 200,
    PlotStyle -> {Hue[0.3], Hue[0.6], Hue[0.]}];

-----
Selwyn Hollis
http://www.math.armstrong.edu/faculty/hollis


On Dec 19, 2003, at 6:57 AM, Maano Aunapuu wrote:

> Hi,
>
> I'm a biologist who tries to use Mathematica to model population
> fluctuations. When trying a set of differential equations that are
> supposed to work properly I get the problems with keeping the system
> persistent or avoiding equations taking negative values (which in a
> biological population is absurdity).
> ------------------------------------
> (*Defining the system*)
> sol := NDSolve[{V'[t] == U[t] (1 - V[t]/Kv[t]) - (A V[t] H[t])/(V[t]
> +\[Alpha] M[t] + B),
>       M'[t] == u[t] M[t] (1 - M[t]/Km) - (\[Alpha] A M[t] H[t])/(V[t]
> +\[Alpha] M[t] + B),
>       H'[t] == R H[t] (A (V[t] + \[Alpha] M[t])/(V[t] + \[Alpha] M[t]
> +B) - G[t]),
>       V[0] == 1000,
>       M[0] == 1000,
>       H[0] == 2},
>     {V, M, H}, {t, 0, 10000}, MaxSteps -> 10000]
>
> (*Defining the constant parameters*)
> A = 15;
> \[Alpha]  = 0.5;
> B = 70;
> Km = 1000;
> R = 10.7/A;
>
> (*As I want to simulate seasons with winter during 0=<t<10 and with
> summer 10<t<12
> I define a parameter m that will have only values 0 - 12*)
> m[t_] := (t - 12*IntegerPart[t/12]);
>
> (*Now I define the parameters whose values depend on t (via m) *)
> U[t_] := 0 /; m[t] \[LessEqual] 10;
> U[t_] := 10000/2 /; m[t] > 10;
> u[t_] := 0 /; m[t] \[LessEqual] 10;
> u[t_] := 12/2 /; m[t] > 10;
> G[t_] := 0.63*A /; m[t] \[LessEqual] 10;
> G[t_] := 0.44*A /; m[t] > 10;
> Kv[t_] := 100 /; m[t] \[LessEqual] 10;
> Kv[t_] := 1000 /; m[t] > 10;
>
> (*And now I wish to see how the system behaves, so I plot it)
> sol1 = sol; Plot[{V[t] /. sol1, M[ t] /. sol1, H[t] /. sol1}, {t, 0,
> 19}, PlotRange -> All, PlotPoints -> 200, PlotStyle -> {Hue[0.3],
> Hue[0.6], Hue[0.0]}];
> -----------------------------
> And now the problems:
>
> 1)
> I get the error message:
> step size is effectively zero; singularity or stiff system suspected.
>
> What does it mean exactly?
> I guess the reason for the message is due H[t] taking negative values,
> and therefore the two other variables V[t]&M[t] shooting into infinity.
>
> As this is a biological system, the negative values are absurd. The
> lowest realistic value should be zero. Is it possible somehow to tell 
> to
> Mathematica it? Or to set some boundaries? So this is my main question
> -how to prohibit the variables reaching negative values?
>
> When I play around with StartingStepSize and decrease it, I get rid of
> the error message. But the results depend and differ very much 
> depending
> on its value. If  StartingStepSize-> 1/20 then H[t]->0, and
> V[t]=M[t]=1000 whereas in case of StartingStepSize-> 1/40 H[t]->0,
> V[t]=0.00004, M[t]=0.22. Why do these results differ so much?
>
> 2)
> If I change the initial conditions to MaxSteps -> 10000,
> StartingStepSize -> 1/20, Method -> {"FixedStep", Method ->
> "LinearlyImplicitModifiedMidpoint"} then I get the system to be
> persistent, even if all the variables now and then reach negative
> values. And I guess therefore I can not relay on the results.
>
> thanks a lot for any help
> Maano
> =¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=
> Maano Aunapuu
> MSc, PhD student
>
> Animal Ecology
> Dept. of Ecology and Environmental Science
> Umeå University
> SE -901 87 Umeå, SWEDEN
> e-mail: maano.aunapuu at eg.umu.se
> http://www.eg.umu.se/sve/personal/maano_aunapuu.htm
> http://www.norrmejerier.se/plupp/
> =¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=
>
>

• References:
  • NDSolve: avoiding negative values?
    • From: Maano Aunapuu <maano.aunapuu@eg.umu.se>