I really don't understand "The Shooting Method"

Question

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.

Answer 1

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$.

Answer 2

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) $$