Mathematica Asked on May 5, 2021
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?
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.
Correct answer by Michael E2 on May 5, 2021
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}]
Answered by Nasser on May 5, 2021
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}]
Answered by Carl Woll on May 5, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP