I really don't understand "The Shooting Method"

I am attempting to learn from a textbook that has the following question:

The boundary-value problem $$y'' = 4(y-x), \qquad 0 \leq x \leq 1, \qquad y(0)=0, \, \, \, y(1)=2$$ has the solution $$y(x) = e^2(e^4 -1)^{-1} (e^{2x}-e^{-2x})+x$$. Use the linear shooting method to approximate the solution, and compare the results to the actual solution. Also, $$h = \frac{1}{2}$$.

This is not mentioned in the question itself, but by looking at the solution, I believe the solution wants us to give the numerical approximation of $$y(0.5)$$. Thank you.

So I'm mostly just focusing on the "approximate the solution" section. I've read so many textbooks but I just can't get my head around it, so I'm having to ask on here. I see that you make an initial approximation $$y'(a) = \lambda$$. And then you use a linear equation in the form:

$$\lambda = \frac{\beta - y_o(b)}{z(b)}$$ to determine $$\lambda$$.

I also think we can use:

$$y_\lambda(x) = y_0(x) + \lambda z(x)$$

I also know you get two equations. I have done it in this manner:

$$y'(x) = z(x)$$ $$z'(x) = 4y(x)-4x$$

And we then want to find $$y_0(x)$$ and $$z(x)$$ I think?. Using Euler's or Runge-Kutta or something? Let us use Runge-Kutta 4 for the sake of this post.

Any help is appreciated. I'm also stuck on the non-linear version, but I'm just going to try and get my head around the linear version for now... I also know this isn't the best written question, sorry. I'm just very lost.

Edit: As nobody has gotten the answer yet it seems, I thought I would post it. The textbook answer is: "0.82432432" as the approximation. Anybody know how to get that result? Thank you.

• For those, like me, that haven't an idea of what linear shooting method is : en.wikipedia.org/wiki/Shooting_method May 20 '19 at 13:27
• What is the question? Do you want someone to explain the shooting method to you and how it works? Or did you want help with your example? May 20 '19 at 14:08
• I have read many things explaining how the shooting method works, so I am good with that. I just do not understand how to actually work through an example, so I would like direct help with that, and what the final answer is. Hopefully this will provide a basis for me to go on and do other examples myself. May 20 '19 at 14:12

The value of $$y'(a)$$ depends on the method you use and the step size. However, because your system is linear, solving for $$y'(a)$$ is just a single step linear interpolation/extrapolation. Alternatively, you can keep track of $$y'(a)$$ during your iterations (since the DE is linear, every step of the approximate solution is going to look like $$y=y_0+y'(a)y_1$$) and solve for $$y'(a)$$ at the end.

First, convert your second-order equation to a first-order equation: $$\begin{pmatrix}y\\y'\end{pmatrix}' =\begin{pmatrix}y'\\4(y-x)\end{pmatrix}$$

For example, if you use Euler's method with step size $$h=1$$, using $$y(0)=0$$ and $$y_{approx}'(1)=\lambda$$ $$\begin{pmatrix} y_{approx}(x=1)\\y'_{approx}(x=1) \end{pmatrix}= \begin{pmatrix} 0\\\lambda \end{pmatrix}+ \begin{pmatrix} \lambda\\0 \end{pmatrix}$$ So $$\lambda_1=0$$ this scheme gives $$y_{approx}(1;\lambda_1)=0$$ and for $$\lambda_2=10$$ (say) we have $$y_{approx}(1;\lambda_2)=10$$. A linear interpolation gives the next candidate $$\lambda=\frac{y(1)-y_{approx}(1;\lambda_1)}{y_{approx}(1;\lambda_2)-y_{approx}(1;\lambda_1)}=2,$$ and of course you find $$y_{approx}(1;\lambda=2)=2$$.

As another example, again with Euler's method, this time $$h=\frac14$$. So \begin{align*} \begin{pmatrix}y_{approx}(1/4)\\y'_{approx}(1/4)\end{pmatrix} &=\begin{pmatrix}0\\\lambda\end{pmatrix} +\frac14\begin{pmatrix}\lambda\\0\end{pmatrix} =\begin{pmatrix}\frac14\lambda\\\lambda\end{pmatrix}\\ \begin{pmatrix}y_{approx}(2/4)\\y'_{approx}(2/4)\end{pmatrix} &= \begin{pmatrix}\frac14\lambda\\\lambda\end{pmatrix}+\frac14 \begin{pmatrix}\lambda\\\lambda-1\end{pmatrix} =\begin{pmatrix}\frac12\lambda\\\frac54\lambda-\frac14\end{pmatrix}\\ \begin{pmatrix}y_{approx}(3/4)\\y'_{approx}(3/4)\end{pmatrix} &= \begin{pmatrix}\frac12\lambda\\\frac54\lambda-\frac14\end{pmatrix}+\frac14 \begin{pmatrix}\frac54\lambda-\frac14\\2\lambda-2\end{pmatrix} =\begin{pmatrix}\frac{13}{16}\lambda-\frac1{16}\\\frac74\lambda-\frac34\end{pmatrix}\\ y_{approx}(1) &=\left(\frac{13}{16}\lambda-\frac1{16}\right)+\frac14\left(\frac74\lambda-\frac34\right)=\frac54\lambda-\frac14 \end{align*} So in the end you find $$\lambda=\frac95$$.

With RK4 and $$h=\frac12$$, you get \begin{align*} k_1&=(\frac\lambda2, 0)\\ k_2&=(\frac\lambda2, \frac\lambda2 - \frac12)\\ k_3&=(\frac58\lambda - \frac18, \frac\lambda2 - \frac12)\\ k_4&=(\frac34\lambda - \frac14, \frac54\lambda - \frac54) \end{align*} So $$(y_{approx}(h),y'_{approx}(h))=\left(\frac7{12}\lambda - \frac1{12},\frac{37}{24}\lambda - \frac{13}{24}\right)$$ and a further RK4 step yields $$y_{approx}(1)=\frac{259}{144}\lambda - \frac{115}{144}.$$ So $$\lambda=\frac{403}{259}\approx1.556$$, $$y_{approx}(h)=\frac{61}{74}\approx 0.8243$$.

• Thank you so much for this. This is what I was looking for. During RK4, on $k_1 = h*f(x_n, y_n)$, what function $f(x_n,y_n)$ are you using here to get $k_1 = h*(\lambda, 0)$? May 20 '19 at 16:26
• I'm doing it on $f(x,(y,y'))=(y,y')'=(y',4(y-x))$. May 20 '19 at 16:27
• Ah, I see! So during $k_2 = (\frac{\lambda}{2}, k_1 - h)$, how does this simplify down to what you have? The $k_1$ being of the form $( , )$ throws me off. And on $k_3$, why isn't the LHS of the output $\frac{\lambda}{2}$ like the others, if the output for the LHS is always $h*y'$. I'm really close to understanding it... thank you. These are my last questions I think! May 20 '19 at 16:38
• For $k_2$, we evaluate $f(x=\frac14, y=\frac14\lambda,y'=\lambda)$ and get the second coordinate ($y''$) to be $4(y-x)=4(\frac14\lambda-\frac14)=\lambda-1$. Multiplying by $h$ gives the $\frac12\lambda-\frac12$. May 20 '19 at 16:49
• For $k_3$, we are evaluating $f$ at $x=\frac14$ and $(y,y')=(0,\lambda)+\frac12 k_2=(\frac14\lambda,\frac54\lambda-\frac14)$. So $f=(\frac54\lambda-\frac14, \frac14\lambda-\frac14)$ and hence $k_3=hf$ is what I wrote above. May 20 '19 at 17:04

Linear case

In principle you could also a second solution $$y_1$$ with $$y_1'(0)=1$$. Then we know that any other solution is an affine combination, $$y=y_0+λ(y_1-y_0)$$.

What your text does is to use $$z=y_1-y_0$$ from the beginning. The differential equation for $$z$$ is thus $$z''=y_1''-y_0''=4(y_1-x)-4(y_0-x)=4(y_1-y_0)=4z$$ with initial condition $$z'(0)=y_1'(0)-y_0'(0)=1$$.

In the numerical solution you can use any solver and solve the equations for $$y_0$$ and $$z$$ either separately or as a system. The quality of the solution will of course depend on the quality of the solver. If you use a method with adaptive step size, then it is better to solve both equations simultaneously.

Non-linear case

What you do exactly depends on what root finding method you want to employ. If $$\phi(t; t_0,x_0)$$ denotes the flow/global solution/general IVP solution with initial condition $$\phi(t_0; t_0,x_0)=x_0$$, then $$λ\mapsto f(λ)=\phi_1(b; a,(0,λ))-β$$ is a scalar function. You may now apply methods like bisection, secant method, Brent, Newton etc. like for any other scalar function.

You get the complication that the numerical solver might introduce additional noise in the evaluation of $$f$$. The whole interval where the evaluation of $$f$$ gives values below that noise level could potentially be the location of the root.

To get at least the search direction somewhat correct it would be helpful to integrate all the occurring solutions with the same time subdivision, which again means either a fixed-step integration with the same step size or a simultaneous integration with variable step size. For instance for the secant method you need two solutions $$y_k,y_{k+1}$$ for initial values $$y_k'(0)=λ_k$$. To the general purpose ODE integrator you would then pass the combined vector, here of dimension 4, and compute the derivatives $$(\dot y_k,\dot y_{k+1})=(f(t,y_k),\,f(t,y_{k+1}))$$ for both components separately. For the Newton method, you would need a derivative in the linear approximation $$y_k+ε v_k$$. While you could use difference quotients, you can also solve for $$v_k$$ directly as $$\dot y+ε\dot v=f(t,y+εv)=f(t,y)+ε\partial_yf(t,y)v+O(ε^2) \\\implies (\dot y_k,\dot v_k) =(f(t,y_k),\,\partial_yf(t,y_k)v_k)$$

• Thank you. This is as far as I got in my own working. Say we have to use Runge-Kutta 4 to solve for $y_0(x)$ and $z(x)$. What would the final answer be, where $h=\frac{1}{2}$? The textbook has the final answer (using RK4 I think) but I am unsure how to get there. May 20 '19 at 14:17
• The final answer in the text-book is 0.82432432. I have edited this into the post. Would you mind showing me how I can get to that final solution? Thank you. May 20 '19 at 14:43
• No, I can't. The derivative $y'(0)$ for the exact solution $y(x)=x+\frac{\sinh(2x)}{\sinh(2)}$ has the value $1.5514411295435664$, and the numerically computed value for $h=0.5$ is $1.555984555984556$. I do not see where your value could occur. May 20 '19 at 14:59
• No, you do one loop for the system with state vector $u=(u_1,u_2,u_3,u_4)=(y,y',z,z')$ and derivatives function $f(x,u)=(u_2, 4(u_1-x), u_4, 4u_3)$. May 20 '19 at 15:33
• You get $y_0(x)=u_1(x)$ and $z(x)=u_3(x)$, as per the construction of the state vector. May 20 '19 at 16:17