0
$\begingroup$

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.
$\endgroup$
  • $\begingroup$ This system of equations has an exact solution. $\endgroup$ – Alex Trounev Feb 7 at 15:38
  • $\begingroup$ 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. $\endgroup$ – march Feb 7 at 17:17
2
$\begingroup$

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[1]], 
  rho -> Function[{t}, -(a21/(2 t^2)) + C[2] - 
     a22 t^(1 + a12) (-a11 t)^(-1 - a12) C[1] Gamma[1 + a12, -a11 t]],
   npsi -> Function[{t}, E^(a31 t + a32 Log[t]) C[3]]}}*)

 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[1], C[2], C[3]}];
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.6458`100*10^-23; a12 = -1.49912`100; a21 = \
-3.89898`100*10^36; a22 = 3.97887`100*10^-23; a31 = 
 SetPrecision[9/(78125000000000000000000 \[Pi]), 
  100]; a32 = -1.49912`100; neta0 = 0.0025`100; rho0 = 
 4.0`100*10^6; npsi0 = 0; t0 = 1.0`100*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"]}

fig1

$\endgroup$
  • $\begingroup$ Thanks, Alex Trounev. $\endgroup$ – D.Nanda Feb 8 at 4:11
  • $\begingroup$ @D.Nanda If the answer solves your problem, you can click the checkmark sign to accept it :) $\endgroup$ – xzczd Feb 8 at 5:08
  • $\begingroup$ @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. $\endgroup$ – D.Nanda Feb 8 at 5:47
  • $\begingroup$ @D.Nanda See update. $\endgroup$ – Alex Trounev Feb 8 at 12:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.