I am solving the initial value problem of two coupled ODEs numerically. The equations are

eqs = {(2.74315901125834` - 2.771065822236567`/x^2) y1[x] + 2.743156268102072`*^-6 y2[x] + 
    2.7989379082285724`*^-18 y1''[x] == 0, 
  1.` y1[x] + 1.0000000000634328` y2[x] + 1.0203348386583513`*^-12 y2''[x] == 0, 
  y1[0.5] == 0, y2[0.5] == 1, y1'[0.5] == 0, y2'[0.5] == 0}

and I want to solve it with NDSolve

sol = NDSolve[eqs, {y1, y2}, {x, 0.5, 1.5}]

However, the results returned with errors

NDSolve::ndsz: At x == 0.5000000000000003`, step size is effectively zero; singularity or stiff system suspected.

and the result is quite large, which is wrong.

I think the problem is from the tiny number, with the orders of $10^{-12}$ and $10^{-18}$.

I looked for solutions on the Internet and learned that using "ImplicitRungeKutta" method may make sense. So I tried

sol = NDSolve[eqs, {y1, y2}, {x, 0.5, 1.5}, 
  Method -> {"FixedStep", 
    Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 10}}, 
  StartingStepSize -> 1/100]

and I tried different settings of StartingStepSize and DifferenceOrder, but they all failed. For large StartingStepSize such as $10^{-3}$, it returns

NDSolve::ndcf: Repeated convergence test failure at x == 1.004`; unable to continue.

the same with large DifferenceOrder. If I choose "DifferenceOrder" -> 10" and "StartingStepSize -> 1/100", there are no errors, but the result is also quite large, which is in the order of $10^{102}$.

I wonder whether there are some methods to deal with tiny numbers in solving ODEs, or maybe I should use Python or Matlab and write the code using Runge-Kutta method which may perform better.

  • 3
    $\begingroup$ you can try adding Method -> {"StiffnessSwitching"} and see what happens. But the correct way to handle this is to scale/normalize your ode itself. $\endgroup$
    – Nasser
    Commented May 8 at 4:25
  • $\begingroup$ @Nasser Thank you for your comment. I have tried using Method -> {"StiffnessSwitching"}, but two hours have passed and the results haven't been calculated. I wonder if I can scale these ODEs because the tiny coefficient is not an overall coefficient. $\endgroup$
    – Link
    Commented May 8 at 9:19
  • $\begingroup$ @Link This is a linear system. Have you any chance to estimate solution around initial conditions? $\endgroup$ Commented May 8 at 14:20
  • 1
    $\begingroup$ Eliminating the boundary conditions and replacing x^2 by c allows DSolve to give a very rough approximation involving terms like E^(9.98524*10^7 x). So, you should expect enormous answers from NDSolve. $\endgroup$
    – bbgodfrey
    Commented May 8 at 17:34


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.