Finding decaying solution using shooting method

Question

Consider the following simple model equation: $$ \frac{d^2y}{dx^2}=y-y^2,$$ which satisfies the decaying boundary condition $$ \lim_{x\to \infty}y(x)=0,\quad \lim_{x\to \infty }y'(x)=0.$$ The solution to the above boundary value problem is given by $$ y(x)=\frac{3}{2}\text{sech}^2(\frac{x}{2}).$$ In particular, $y(0)=3/2$ and $y'(0)=0$.

Now I want to use the shooting method to find a good numerical approximation to the above solution. First, observe that the solution is an even function. Thus we have $y'(0)=0$ and we only need to guess $y_0$ with $y(0)=y_0$. Let $y(x,y_0)$ denote the function that solves the initial value problem $$ \frac{d^2y}{dx^2}=y-y^2,\quad y(0)=y_0,\quad y'(0)=0.$$ My question is that for the decaying boundary condition, what is my target condition? Is it $$ F(y_0):=|y(L,y_0)|+|y'(L,y_0)|=0$$ for sufficiently large $L$?

The following figure shows the $F(y_0)$ versus $y_0$ curve plotted by Mathematica when $L=10$. It seems that there is a singularity at $y_0=3/2$ such that the value of $F(y_0)$ is extremely large when $y_0$ passes through $3/2$. So I do not think that above target condition is a good one. $F(y_0)$ versus $y_0$ curve when $L=10$" />

Answer 1

You are moving towards a saddle point in phase space. The only way to reach it is via one of the branches of the stable manifold. As in the linearization at $y=0$ the quadratic term falls away, one is left with $y''=y$ which has as solution $y(x)=Ce^{-x}+De^x$. For the stable solution you need $D=0$. Every error, method related or floating point, will give a non-zero value to $D$ which leads to a solution rapidly growing away from the $x$-axis.

Thus the condition one should enforce at $x=L$ is $D=0$ or $$y(L)+y'(L)=0$$ in this first approximation. The equation can thus be solved as a boundary value problem with this condition and continued as $y(L)e^{L-x}$ for $x>L$.

Further, one can integrate the equation once, after multiplying with $2y'$, $$ y'^2=y^2-\frac23y^3+C. $$ For the desired solution this gives $C=0$ and thus either $y(0)=0$ or $y(0)=\frac32$. Along the solution to the saddle point, at every point one has thus the identity $$ y'(x)+y(x)\sqrt{1-\frac23y(x)}=0. $$ This gives the exact second boundary condition without using a perturbation or other approximation method.

Using the substitution $u=\sqrt{1-\frac23y(x)}\iff y=\frac32(u^2-1)$, $y'=3uu'$, this leads quickly to the exact solution.

Using these pieces to produce an improved shooting graph gives

f = lambda y,t: [y[1],y[0]*(1-y[0])]
v = np.linspace(0,1.6,201)
def shoot(v,L):
    y = odeint(f,[v,0],[0,L],atol=1e-12,rtol=1e-12,hmax=0.1)[-1]
    return y[1]+y[0]*abs(1-2/3*y[0])**0.5
for L in np.linspace(4,7,7): 
    plt.plot(v,[shoot(vv,L) for vv in v])
plt.grid(); plt.show()

As one can see, with increasing $L$ the sensitivity of the target value from the initial value grows rapidly, meaning the solutions reach the saddle point and escape along the unstable manifold. The velocity of the escape then gets increased by the non-linear term for larger $|y|$.