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

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) &= -{1\over 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-{t\over 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}]

Mathematica graphics

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


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}]

Mathematica graphics


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}]

enter image description here