Services & ResourcesWolfram ForumsMathGroup Archive

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

NDSolve: avoiding negative values?

• To: mathgroup at smc.vnet.net
• Subject: [mg45162] NDSolve: avoiding negative values?
• From: Maano Aunapuu <maano.aunapuu at eg.umu.se>
• Date: Fri, 19 Dec 2003 06:57:29 -0500 (EST)
• Reply-to: maano.aunapuu at eg.umu.se
• Sender: owner-wri-mathgroup at wolfram.com

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/
=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=¤=

• Follow-Ups:
  • Re: NDSolve: avoiding negative values?
    • From: Selwyn Hollis <sh2.7183@misspelled.erthlink.net>