# NDSolve for coupled differential equations with large coefficient

I have three coupled differential equations to solve. Those equations are

equation1=neta'[t] == -(7.6458*10^-23) neta[t] - (1.49912 neta[t])/t

equation2=rho'[t] == -(3.89898*10^36/t^3) + (3.97887*10^-23) neta[t]

equation3=npsi'[t] == (9 npsi[t])/(78125000000000000000000 \[Pi]) - (1.49912 npsi[t])/t


I have tried to solve these equations by using

soln = NDSolve[{equation1,equation2, equation3, neta[10^-10] == 0.0025,rho[10^-10]== 4*10^6 ,npsi[10^-10] == 0}, {neta,rho,npsi}, {t, 10^-10, 1}, AccuracyGoal -> Autometic, Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5,
"StiffnessTest" -> True}, MaxStepSize -> Autometic, PrecisionGoal -> Autometic]


It is showing the error

NDSolve::ndsz: At t == 1.*^-10, step size is effectively zero;singularity or stiff system suspected.


However, if I just skip the second equation and solve for the other two then it can be solved by the same code. Even if I skip the first part of equation2 then also it is solving. So I think the problem is in the first part of the equation2, a very large number is coming there and it fails to solve.

Can you please help in this regards?

Thank you.

Suppose I have a set of equations

eqn1=neta'[t] == -7.6458*10^-23 neta[t] - 3.566*10^-18 neta[t] Sqrt[rho[t]]

eqn2=rho'[t] == 3.97887*10^-23 neta[t] - 4.75467*10^-18 rho[t]^(3/2)

eqn3=npsi'[t] == 3.66693*10^-23 neta[t] - 3.566*10^-18 npsi[t] Sqrt[rho[t]]


Then if I do it in a similar way then

soln = NDSolve[{eqn1,eqn2, eqn3, neta[10^-10] == 0.0025,rho[10^-10]== 4*10^6 ,npsi[10^-10] == 0}, {neta,rho,npsi}, {t, 10^-10, 1}, AccuracyGoal -> Autometic, Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5, "StiffnessTest" -> True}, MaxStepSize -> Autometic, PrecisionGoal -> Autometic]


It shows same kind of error

NDSolve::ndtol: Tolerances requested by the AccuracyGoal and PrecisionGoal options could not be achieved at t == 1.*^-10.

• This system of equations has an exact solution. – Alex Trounev Feb 7 '19 at 15:38
• If you have equations where there are constants that whose ratios are 10^59, you have chosen a bad set of units. It seems to me that you should rescale t in this case by something like 10^23. – march Feb 7 '19 at 17:17

## 1 Answer

We write the system of equations in the form

equation1 = neta'[t] == a11*neta[t] + a12*neta[t]/t;

equation2 = rho'[t] == a21/t^3 + a22* neta[t];

equation3 = npsi'[t] == a31* npsi[t] + a32 npsi[t]/t;


The exact solution allows us to see the source of the problem.

soln = DSolve[{equation1, equation2, equation3}, {neta, rho, npsi}, t]

(*Out[]= {{neta -> Function[{t}, E^(a11 t + a12 Log[t]) C],
rho -> Function[{t}, -(a21/(2 t^2)) + C -
a22 t^(1 + a12) (-a11 t)^(-1 - a12) C Gamma[1 + a12, -a11 t]],
npsi -> Function[{t}, E^(a31 t + a32 Log[t]) C]}}*)

netaC = First[neta[t] /. soln]; rhoC =
First[rho[t] /. soln]; npsiC = First[npsi[t] /. soln];
eq1 = netaC /. t -> t0; eq2 = rhoC /. t -> t0; eq3 = npsiC /. t -> t0;
sol = Solve[
eq1 == neta0 && eq2 == rho0 && eq3 == npsi0, {C, C, C}];
netaT = netaC /. sol; rhoT = rhoC /. sol; npsiT = npsiC /. sol;

a11 = -(7.6458*10^-23); a12 = -1.49912; a21 = -3.89898*10^36; a22 =
3.97887*10^-23; a31 =
9./(78125000000000000000000 \[Pi]); a32 = -1.49912; neta0 = 0.0025; \
rho0 = 4.*10^6; npsi0 = 0.; t0 = 1.*10^-10;
{netaT, rhoT, npsiT}
{{(2.55117*10^-18 E^(7.6458*10^-33 - 7.6458*10^-23 t))/
t^1.49912}, {(-1.94949*10^56 + 2.44065*10^-51 I) + 1.94949*10^36/
t^2 - 9.28269*10^-52 Gamma[-0.49912, 7.6458*10^-23 t]}, {0.}}


To get a numerical solution, we need to increase WorkingPrecision -> 100and all coefficients and initial data set with this precision.

a11 = -7.6458100*10^-23; a12 = -1.49912100; a21 = \
-3.89898100*10^36; a22 = 3.97887100*10^-23; a31 =
SetPrecision[9/(78125000000000000000000 \[Pi]),
100]; a32 = -1.49912100; neta0 = 0.0025100; rho0 =
4.0100*10^6; npsi0 = 0; t0 = 1.0100*10^-10;
equation1 = neta'[t] == a11*neta[t] + a12*neta[t]/t;
equation2 = rho'[t] == a21/t^3 + a22* neta[t];
equation3 = npsi'[t] == a31* npsi[t] + a32 npsi[t]/t;
soln = NDSolve[{equation1, equation2, equation3, neta[t0] == neta0,
rho[t0] == rho0, npsi[t0] == npsi0}, {neta, rho, npsi}, {t, t0,
1.100}, WorkingPrecision -> 100, MaxSteps -> 10^6]

{LogLogPlot[neta[t] /. soln, {t, t0, 1}, PlotRange -> All,
PlotLabel -> "neta"],
LogLogPlot[Abs[rho[t]] /. soln, {t, t0, 1}, PlotRange -> All,
PlotLabel -> "rho"]}
`
• Thanks, Alex Trounev. – D.Nanda Feb 8 '19 at 4:11
• @D.Nanda If the answer solves your problem, you can click the checkmark sign to accept it :) – xzczd Feb 8 '19 at 5:08
• @Alex This solves the problem partially. I mean, if we deal with a system of equations where analytical solutions are not possible and I have to deal with such small numbers then it shows the same kind of error. – D.Nanda Feb 8 '19 at 5:47
• @D.Nanda See update. – Alex Trounev Feb 8 '19 at 12:26