NDSolve inside NDSolve in Mathematica

Question

I need to use NDSolve which in turn uses the solution from another ODE as function in terms of output from another NDSolve.

If I use the exact solution from the first differential equation inside the NDSolve, it's OK. But when I use the same solution in the form of function (which uses InterpolatingFunction) it does not work.

I believe, it's got to do with the structure of NDSolve output. Could anyone please enlighten me on this. Will be of great help!

The code is:

---

feq = 2 V alpha fip F''[fi] - (V^2 - (V^2 + sigma - 2 fi) (F'[fi])^2 + (F'[fi])^4

Frange[lo_, hi_] := Module[{fii, sol}, sol = NDSolve[{(feq == 0 /.fi -> fii), F[0] == 0}, F, {fii, lo, hi}]]

eqpois = fi''[x] == ne[x] - F[fi[x]]/.sol

NDSolve[{eqpois, fi'[0] == 0, fi[0] == 0}, fi, {x,0,1}]

Here in order to find F[phi], I need to solve the 1st diff eq that is feq, which is solved by NDSolve inside the function Frange[lo,hi]. The solution is then used inside the second equation eqpois, which has to be solved using NDSolve again. The problem comes up in the second NDSolve, which does not produce the result. If I use the analytical solution of F[phi] in eqopis, then there is no problem.

---

Example Problem

I have done a little experiment with this. Let's take an example of coupled ODEs

1st eqn : dg/dx = 2f(g) with initial condition g(0) = 1

The function f(y) is a solution from another ODE, say,

2nd eqn : df/dy = 2y with IC f(0) = 0

The solution of the 2nd ODE is f(y) = y^2 which when put into the the 1st ODE becomes

dg/dx = 2 g^2 and the final solution is g(x) = 1/(1-2x)

The issue:

When I use DSolve, it finds the answer correctly

In[39]:= s = DSolve[{f'[y] == 2 y, f[0] == 0}, f, y]

Out[39]= {{f -> Function[{y}, y^2]}}

In[40]:= ss = DSolve[{g'[x] == 2 (f[g[x]]/.First@s), g[0] == 1}, g, x]

Out[40]= {{g -> Function[{y}, 1/(1 - 2 x)]}}

The problem comes when I use NDSolve

In[41]:= s = NDSolve[{f'[y] == 2 y, f[0] == 0}, f, {y, 1, 5}]

Out[41]= {{f -> InterpolatingFunction[{{1., 5.}}, <>]}}

In[42]:= ss1 = NDSolve[{g'[x] == 2 (Evaluate[f[g[x]]/.First@s1]), g[0] == 1}, g, {x, 1, 2}]

Out[42]= {}

The erros are:

During evaluation of In[41]:= InterpolatingFunction::dmval: Input value {2.01726} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

During evaluation of In[41]:= InterpolatingFunction::dmval: Input value {2.01726} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

During evaluation of In[41]:= InterpolatingFunction::dmval: Input value {2.04914} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

During evaluation of In[41]:= General::stop: Further output of InterpolatingFunction::dmval will be suppressed during this calculation. >>

During evaluation of In[41]:= NDSolve::ndsz: At y == 0.16666654771477857, step size is effectively zero; singularity or stiff system suspected. >>

During evaluation of In[41]:= NDSolve::ndsz: At y == 0.16666654771477857, step size is effectively zero; singularity or stiff system suspected. >>

Any help in this regard will be highly appreciated!

--- Madhurjya

Answer 1

I got your simple example to work with a little mod ..

 f0 = First@First@DSolve[{f'[y] == 2 y, f[0] == 0}, f, y]
 g0 = g /. 
     First@First@DSolve[{g'[x] == 2 (f[g[x]] /. f0), g[0] == 1}, g, x]
 fn = f /. First@First@NDSolve[{f'[y] == 2 y, f[0] == 0}, f, {y, 0, 10}]
 gn = g /. 
      First@First@
         NDSolve[{g'[x] == 2 (fn[g[x]]), g[0] == 1}, g, {x, 0, 9/20}]
 GraphicsRow[{
  Plot[{g0@x, gn@x}, {x, 0, 9/20}, 
       PlotStyle -> {{Thick, Black}, {Thin, Red, Dashed}}],
  Plot[{f@x /. f0, fn@x}, {x, 0, 2}, 
       PlotStyle -> {{Thick, Black}, {Thin, Red, Dashed}}]}]

note we need to ensure the y range in the first NDSolve is sufficient to cover the expected range of g from the second. That is where all those interpolation range errors come from.