Numerically solving an ODE whose right-hand side involves an integral

Question

In a 1955 paper on surfaces of constant negative Gaussian curvature, Marc-Henri Amsler defines a function $w(x)$ by giving the ordinary differential equation
$$w'(x)=frac{1}{x}int_0^x sin(w(t))dt,$$
with the initial condition $w(0)=pi/2$.  So the slope $w'(x)$ is the average value of $sin(w(x))$ over the interval $[0..x]$.  In particular, we have $w'(0)=sin(pi/2)=1$.

I asked Mathematica to compute $w(x)$ for me numerically by using NDSolve:

NDSolve[{w'[x] == NIntegrate[Sin[w[t]], {t, 0, x}]/x, w[0] == Pi/2}, w[x],
             {x, -2, 2}]

but Mathematica (version 10.0.2) responds with a slew of errors, the first batch of which repeatedly report that

NIntegrate::nlim: t = x is not a valid limit of integration. >>

Is this ODE outside of Mathematica's region of expertise, or I am doing something stupid in how I have phrased my Mathematica request?

Michael E2 · Accepted Answer

Being a mathematician, I resist fudging by cutting off the singularity by some small eps = 10^-12.  But if you're an engineer, you should be satisfied @Nasser's answer, which likely yields single-precision accuracy, which is pretty good after all.

A little analysis shows
$$eqalign{
w''(x) &= -{1over x^2} int_0^x sin w(t) ; dt + {sin w(x) over x} cr
 &=  int_0^x {sin w(x) - sin w(t) over x^2} ; dt cr
 &= int_0^x {(x-t) over x^2}{sin w(x) - sin w(t) over x-t} ; dt cr
 &buildrel rm{MVT} over = int_0^x {(x-t) over x^2},cos w(xi) , w'(xi),(x-t)  ; dt cr
 &= int_0^x {left(1-{tover x}right)^2},cos w(xi) , w'(xi)  ; dt cr
}$$
where $xi = xi_t$ is between $x$ and $t$ and depends on $t$ for a fixed $x$. Of the three factors in the integrand as $x rightarrow 0$, the first lies between $0$ and $1$, the cosine approaches zero since $w(xi) rightarrow pi/2$, and $w'(xi)$ approaches $1$.  Therefore $w''(x) rightarrow 0$ as $x rightarrow 0$.
Thus we should define our ODE so that $w''(0) = 0$.  Piecewise is the tool for that.

sol = NDSolve[{w''[x] == Piecewise[{{(Sin[w[x]] - w'[x])/x, x != 0}}],
     w[0] == Pi/2, w'[0] == 1}, w, {x, 0, 4 Pi}];
Plot[w[x] /. sol, {x, 0, 4 Pi}, AxesOrigin -> {0, 0}]

It differs from Nasser's answer by less than 3 * 10^-8 over the interval.

Nasser · Answer

There is a singularity at x=0.
First, differentiate and use Leibniz rule, this gives the non-linear ODE
$$
   x w''(x)+ w'(x) = sin(w(x))
$$
But NDSolve can't solve this due to singularity at $x=0$
eq = x w''[x] + w'[x] == Sin[w[x]];
NDSolve[{eq, w[0] == Pi/2, w'[0] == 1}, w, {x, 0, 1}]
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`.

So using trick to bypass $x=0$ (maybe a better trick exists, but the idea is to avoid zero)
sol = With[{eps = 10^-12},
  NDSolve[{w''[x] == 
     If[x < eps, -w'[x]/eps + Sin[w[x]]/eps, -w'[x]/x + Sin[w[x]]/x], 
    w[0] == Pi/2, w'[0] == 1}, w, {x, 0, 4 Pi}]
  ];
Plot[w[x] /. sol, {x, 0, 4 Pi}, AxesOrigin -> {0, 0}]

Carl Woll · Answer

The equivalent ODE (as given by @Nasser) can be solved by NDSolveValue when using the {"EquationSimplification"->"Residual"} Method:

sol = NDSolveValue[
    {
    x w''[x] + w'[x] == Sin[w[x]],
    w[0] == Pi/2,
    w'[0] == 1
    },
    w,
    {x, 0, 12},
    Method->{"EquationSimplification"->"Residual"}
];

Plot[sol[x], {x, 0, 12}, AxesOrigin -> {0,0}]

Numerically solving an ODE whose right-hand side involves an integral

3 Answers

Add your own answers!

Ask a Question