I need to solve a nonlinear ODE of the form

$$ \frac{d^2 y}{dx^2} + \frac{1}{x}\frac{dy}{dx}-\frac{1}{x^2}\sin(y)\cos(y)+\frac{2}{\alpha}\frac{\sin^2(y)}{x}-\sin(y)=0 $$

numerically, subject to the boundary conditions

$$ y(0)=\pi,\, y(\infty)=0 $$

I'd prefer to do this in python, but I'm having a lot of trouble zeroing in on only the decaying solutions. I tried a shooting method using a near-field asymptotic expansion to bypass the singularity at 0, and I know what the far-field asymptotic looks like up to some undetermined constant. The asymptotics are

$$ y \approx \pi-c_1 x $$ near zero, and

$$ y \approx \frac{c_2}{\sqrt{x}}e^{-x} $$

for large $x$.

Naturally it seems like $c_1$ would be the shooting parameter, but I'm stuck as to how to make judicious use of $c_2$ in some way. The solutions diverge from each other very quickly as $x$ gets even a little large.

The free parameter $\alpha$ is of order $1$. However, small variations in this parameter also trigger large changes in the convergence of solutions.

Any suggestions would be greatly appreciated!

Update: Huge thanks to Lutz' answer below -- I ended up using both a boundary value solver and single shooting depending on the regime that alpha is in. Single shooting is slower, but more reliable, while the boundary value solver is much faster, but more limited in its scope. Using good boundary conditions at the asymptotes ended up being the most important!

  • $\begingroup$ How do you get the square root in the far-field approximation? It is not from WKB, as the linearization $y''-y=0$ has constant coefficients. $\endgroup$ Commented Jun 25, 2023 at 12:22
  • $\begingroup$ If the first derivative term is retained then the solution satisfies a modified Bessel equation, from which the boundary conditions eliminate the growing solution and only keep the decaying one $\endgroup$
    – Ali Shakir
    Commented Jun 26, 2023 at 5:01
  • 1
    $\begingroup$ So you are using $x^2y''+xy'-(1+x^2)y=0$ as approximate equation for the far-field. // If you are sure of your asymptotic approximation, the constant is eliminated in the derivative of $\sqrt xe^xy(x)=c_2$. The resulting differential equation can be taken as boundary condition at some medium large upper bound. // I would use a BVP solver, not single-shooting. But if the upper boundary gets reduced to medium sizes, single shooting should still work well enough. $\endgroup$ Commented Jun 26, 2023 at 6:01
  • $\begingroup$ Ah, I see — so because the derivative of the rhs of the asymptotic expansion is 0 I just need to use that the lhs is approximately zero for some large-ish x? That makes a lot of sense, thank you. $\endgroup$
    – Ali Shakir
    Commented Jun 26, 2023 at 7:07

1 Answer 1


Eliminating the constants, the approximation close to $x=0$ is $y(x)=\pi+xy'(x)$ with $y'$ nearly constant. Or one could multiply the leading terms with $2x^2y'$ and integrate $$ 2x^2y'y''+2xy'^2-2\sin(y)\cos(y)y'\approx 0\\ (xy'(x))^2-\sin^2(y)\approx C $$ with $C=0$ for the given initial condition. So $xy'=\pm\sin(y)\approx \pm(\pi-y)$ and the (backwards) stable variant has the minus sign.

For the far field you decided to consider as leading terms $$ \sqrt{x}y''+\frac1{\sqrt x}y'-\sqrt xy\approx 0 $$ or with $u=\sqrt xy$ $$ u''-u=0 $$ The forwards stable solution satisfies the first-order equation $u'=-u$, which expands to $xy'+(x+\frac12)y=0$.

These conditions can now be inserted into any BVP solving method on some interval $[a,b]$ with $a>0$ small, $b>a$ medium sized. To distribute the singularity of the divisions by $x$ one can consider the second state variable $z(x)=xy'(x)$ so that $z'=xy''+y'=x(y''+\frac1xy')$, that is, $$ y'=\frac1xz\\ z'=\frac1x\sin(y)\cos(y)+x(\sin y-\frac2\alpha\sin^2(y)) $$

For a single shooting method thus $z_a=-\sin(y_a)=\sin(y_a-\pi)$ with a target $z_b+(b+\frac12)y_b=0$. Consider employing in the ODE solver a terminal event on a large circle, to capture divergence early.

For multiple shooting or collocation methods one uses the boundary conditions in equation form.

from numpy import pi, sin, cos, exp, logspace

alpha = 1.1

def ode(x,y): sy = sin(y[0]); return y[1]/x, x*(sy-2*sy*sy/alpha)+sy*cos(y[0])/x 

def bc(ya,yb): return (sin(ya[0])+ya[1]),(yb[1]+(b+0.5)*yb[0])

def solve(a,b):
    x = [a]
    while x[-1] < b: x.append(x[-1]*2)
    x = np.asarray(x); x[-1] = b; 

    y = [ pi*exp(-x), -pi*exp(-x)*x]

    res = solve_bvp(ode, bc, x, y, tol=1e-8)
    print(res.message, f" {len(res.x)} nodes")
    return res
res = solve(1e-4,15)

This works and gives as control output

The algorithm converged to the desired accuracy.  972 nodes

With this mediums small tolerance value the solution is also stable and compatible over combinations of $a\in$[1e-2,1e-4,1e-6], $b\in$[5,10,15,20].

  • $\begingroup$ so it turns out that my difficulties were partly caused by noting down my diff eq incorrectly. I missed a $\frac{1}{x}$ prefactor before the $\sin{y}^2$, which is now fixed in the original question. As far as I can tell, the far-field behavior should still be valid, but I'm trying to figure out how the near-field needs to be modified. Thanks for the detailed answer, it was immensely helpful! $\endgroup$
    – Ali Shakir
    Commented Jun 27, 2023 at 23:20

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.