# Allow Mathematica to solve the differential equations with larger than 16 digits parameters

Consider the system of differential equations

test = ParametricNDSolve[{x'[t] ==
1.1*(1 - x[t])*a*y[t] - x[t]*a*z[t],
y'[t] == -1.1 (1 - x[t])*b*y[t] - c*y[t] + d*Exp[-t/0.02],
z'[t] == -x[t]*b*z[t] - c*z[t] + d*Exp[-t/0.02], y[0.02] == 0,
z[0.02] == 0, x[0.02] == 0.1}, {x, y, z}, {t, 0.02, 0.5}, {a, b, c,
d}]


If I set the parameters a,b,c,d to

av = 1;
bv = 2.75*10^6;
cv = 4*10^7;
dv = 10^16;


then the solutions are obtained. However, setting, for example, d = 10^37, and a = 10^-23, when calculating the solution

solutiontest1 = x[av, bv, cv, dv] /. test


ParametricNDSolve::ndsz: At t == 0.02, step size is effectively zero; singularity or stiff system suspected.

How to resolve this problem?

• @Rebel, usually, we do not accept answers in this site so quickly; this is so that other users can give other (and potentially more refined) answers. A 24-hour waiting period is common. – J. M. will be back soon Nov 14 '19 at 14:45

As the error says, the system is stiff, so you need to use a Method that is appropriate for stiff systems, either "BDF" or "StiffnessSwitching". From experience, I've found that it usually helps to also increase the WorkingPrecision to some high value and use exact numbers, eg 11/10 instead of 1.1 and so on.

For example, the following works:

a = 1;
b = 275/100*10^6;
c = 4*10^7;
d = 3*10^16;
p0 = 2/100;
test = NDSolve[{x'[t] == 11/10*(1 - x[t])*a*y[t] - x[t]*a*z[t],
y'[t] == -(11/10) (1 - x[t])*b*y[t] - c*y[t] + d*Exp[-t/p0],
z'[t] == -x[t]*b*z[t] - c*z[t] + d*Exp[-t/p0], y[p0] == 0,
z[p0] == 0, x[p0] == 1/10}, {x, y, z}, {t, p0, 1/2},
WorkingPrecision -> 25, Method -> "BDF", MaxSteps -> 10^6][];
LogLogPlot[{x[t], y[t], z[t]} /. test // Abs // Evaluate, {t, 0.02,0.5}]


Alternatively, you can also use "StiffnessSwitching", but I found it makes the solution a bit wiggly.

• The desired value of a parameter d is something of an order of 10^37, while of a parameter a is something like 10^-23. With this value, the error persists. – John Taylor Nov 14 '19 at 15:15
• You did not specify in the question what is the "desired" value of d. Please post all necessary information, otherwise any answer will be incomplete... – Rebel-Scum Nov 14 '19 at 15:16
• I have already added it. But initially, I mentioned values d > 10^16, which can be any, not necessarily close to 10^16, so the value 3*10^16 was your assumption for which the answer works. – John Taylor Nov 14 '19 at 15:18