CODEWITHSUNDEEP
wolfram-mathematicamodelingmathematical-expressionsnonlinear-equationmanipulate
Raihanah Nazihah
Raihanah Nazihah

Reputation: 11

How fix error in Mathematica "NDSolve::ndsz: At t == ... step size is effectively zero; singularity or stiff system suspected" when using Manipulate?

I am trying to solve a nonlinear system of differential equations in Wolfram Mathematica using NDSolve and visualize the solution with Manipulate. However, I encounter the following error when running the code. Interestingly, the plot still appears in the output, but the warning message also shows up. Is this normal, or is there something wrong with my Mathematica code? I suspect the issue arises because my system of equations is quite complex and nonlinear, with fractions whose denominators contain variables. Is that the main cause of this error?

This is the code that I made:

 Module[{plt1, plt2, plt3, sol, S0 = SS0, Ih0 = IhIh0, Y0 = YY0}, sol =
   sol = 
    NDSolve[{S'[t] == 
       3.5*S[t]*(1 - (S[t] + Ih[t])/50) - 1*S[t]*Ih[t] - (
        0.34*(S[t])^2*Y[t])/(S[t] + 0.1*Ih[t]),
      Ih'[t] == 
       1*S[t]*Ih[t] - (2.2*(Ih[t])^2*Y[t])/(S[t] + 0.1*Ih[t]) - 
        0.02*Ih[t],
      Y'[t] == (0.9*0.34*(S[t])^2*Y[t])/(S[t] + 0.1*Ih[t]) - (
        0.02*2.2*(Ih[t])^2*Y[t])/(S[t] + 0.1*Ih[t]) - 0.022*Y[t],
      S[t /; t <= 0] == S0, Ih[t /; t <= 0] == Ih0, 
      Y[t /; t <= 0] == Y0}, {S[t], Ih[t], Y[t]}, {t, 0, 500}];
  plt1 = 
   ParametricPlot[{t, S[t]} /. sol, {t, 0, 500}, PlotRange -> All, 
    AspectRatio -> 1, PlotStyle -> {Red, Thick}, Axes -> {t, S}];
  plt2 = 
   ParametricPlot[{t, Ih[t]} /. sol, {t, 0, 500}, PlotRange -> All, 
    AspectRatio -> 1, PlotStyle -> {Green, Thick}, Axes -> {t, Ih}];
  plt3 = 
   ParametricPlot[{t, Y[t]} /. sol, {t, 0, 500}, PlotRange -> All, 
    AspectRatio -> 1, PlotStyle -> {Blue, Thick}, Axes -> {t, Y}];
  Show[plt1, plt2, plt3, ImageSize -> {300, 300}]],
 
 Style["Persamaan differensial :", Bold],
 Style["S'= rS(1-\!\(\*FractionBox[\(S + I\), \
\(K\)]\))-\[Lambda]SI-\!\(\*FractionBox[\(\*SuperscriptBox[\(pS\), \
\(2\)] Y\), \(S + \[Alpha]I\)]\)", Bold],
 Style["I'=\!\(\*FormBox[\(\[Lambda]SI - \
\*FractionBox[\(\*SuperscriptBox[\(cI\), \(2\)] Y\), \(S + \
\[Alpha]I\)] - \[Gamma]I\),
TraditionalForm]\) ", Bold],
 Style["Y'=\!\(\*FormBox[\(\*FractionBox[\(\*SuperscriptBox[\(\[Delta]\
pS\), \(2\)] Y\), \(S + \[Alpha]I\)] - \
\*FractionBox[\(\*SuperscriptBox[\(\[Eta]cI\), \(2\)] Y\), \(S + \
\[Alpha]I\)] - dY\),
TraditionalForm]\)", Bold],
 Delimiter,
 Style["parameters", Bold, 10],
 
 Delimiter,
 Style["initial conditions", Bold, 10],
 {{SS0, 4, "S0"}, 0, 10, .01, ImageSize -> Small, 
  Appearance -> "Labeled"},
 {{IhIh0, 1, "Ih0"}, 0, 10, .01, ImageSize -> Small, 
  Appearance -> "Labeled"},
 {{YY0, 6, "Y0"}, 0, 10, .01, ImageSize -> Small, 
  Appearance -> "Labeled"},
 ControlPlacement -> Left, SynchronousUpdating -> False]```

I got the plot and error message in the output. Is it normal or or is there something wrong with my Mathematica code? How can I resolve this error and still obtain a reliable solution? Are there specific methods or options in NDSolve designed to handle stiff systems? Should I focus on adjusting and finding the right values for S0, Ih0, and Y0 to prevent the error? I have attached the plot generated from Manipulate that illustrates the problem. I would greatly appreciate any explanations and suggestions to solve this issue. Thank you very much! 🙏

Upvotes: 1

Views: 37

Answers (1)

Bill
Bill

Reputation: 3967

Too long for a comment. If I do this

range=40;SS0=4;IhIh0=1;YY0=6;S0=SS0;Ih0=IhIh0;Y0=YY0;
sol=NDSolve[{
  S'[t]==3.5*S[t]*(1-(S[t]+Ih[t])/50)-1*S[t]*Ih[t]-(0.34*(S[t])^2*Y[t])/
    (S[t]+0.1*Ih[t]),
  Ih'[t]==1*S[t]*Ih[t]-(2.2*(Ih[t])^2*Y[t])/
    (S[t]+0.1*Ih[t])-
    0.02*Ih[t],
  Y'[t]==(0.9*0.34*(S[t])^2*Y[t])/
    (S[t]+0.1*Ih[t])-
    (0.02*2.2*(Ih[t])^2*Y[t])/
    (S[t]+0.1*Ih[t])-
    0.022*Y[t],
  S[t/;t<=0]==S0,Ih[t/;t<=0]==Ih0,Y[t/;t<=0]==Y0},{S[t],Ih[t],Y[t]},{t,0,range}];
ParametricPlot[{t,S[t]}/.sol,{t,0,range},PlotRange->All]
ParametricPlot[{t,Ih[t]}/.sol,{t,0,range},PlotRange->All]
ParametricPlot[{t,Y[t]}/.sol,{t,0,range},PlotRange->All]

then the three plots appear without error. Notice the plots show that S[t] and Ih[t] are very close to zero and getting much closer to zero and you have both of those in denominators. I suspect that means that the fraction finally blows up if I increase range beyond 40 and you want the range to go out to 500! Can you carefully go through your three equations and find a way to remove those "near zero" denominators?

In Mathematica any time you use a decimal point that means Mathematica will interpret that constant as a MachinePrecision value with only a small number of known good digits. You can change 0.02 to 2/100 and that will tell Mathematica that this is an exact rational value with infinite precision. But that is not enough. You can tell NDSolve to use WorkingPrecision->256 to try to use 256 digit precision in all calculations, but that will not change 0.02 to be a 256 digit precision constant. But this is still not enough if you cannot fix those denominators.

If I change stuff==morestuff+numerator/denominator to (stuff-morestuff)*denominator==numerator then I can increase the range a little before it fails. It looks like something is happening in your system around t==50 that may be making this fail.

You can also do a Google search for

Mathematica NDSolve "step size is effectively zero"

which should find posts other people have made when they have similar problems and answers to those might possibly help you.

Upvotes: 1

Related Questions