0

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.

Penguin47
  • 57
  • 1
    For those, like me, that haven't an idea of what linear shooting method is : https://en.wikipedia.org/wiki/Shooting_method – Jean Marie 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? – Matthew Cassell 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. – Penguin47 May 20 '19 at 14:12

2 Answers2

2

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

user10354138
  • 33,239
  • Thank you so much for this. This is what I was looking for. During RK4, on $k_1 = hf(x_n, y_n)$, what function $f(x_n,y_n)$ are you using here to get $k_1 = h(\lambda, 0)$? – Penguin47 May 20 '19 at 16:26
  • I'm doing it on $f(x,(y,y'))=(y,y')'=(y',4(y-x))$. – user10354138 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! – Penguin47 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$. – user10354138 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. – user10354138 May 20 '19 at 17:04
  • When we are on our first rotation through RK4, shouldnt $x=0$ and $y=0$ always? As this is their initial values? Where do we find that they equal $\frac{1}{4}$? – Penguin47 May 20 '19 at 17:06
  • I'm not sure what you mean. For $k_1$ we evaluate at our initial value. For $k_2$ we start to move away and use an estimate at $x=\frac12h$, similarly $k_3,k_4$. – user10354138 May 20 '19 at 17:14
  • Thank you so much. It's just clicked. I get it! Thanks. – Penguin47 May 20 '19 at 17:52
  • How do you know $y_{approx}(1) = \frac{259}{144}\lambda - \frac{115}{144}$ basically. Other than that, I know how to do this. So thank you! – Penguin47 May 20 '19 at 18:57
  • Can you explain "A linear interpolation gives the next candidate..". Can you show explicitly what function you are using here? – Penguin47 May 20 '19 at 22:31
  • Please can you add in extra steps to show exactly how you got $y_{approx}(1) = \frac{5}{4} \lambda - \frac{1}{4}$. I do not see how you jump from the final steps of the euler/runge methods to suddenly obtaining functions for $y_{approx}(1)$? For example, I continued the euler method on, and I got $\lambda$ as my final $y{approx}(1)$ component? – Penguin47 May 20 '19 at 22:55
1

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

Lutz Lehmann
  • 126,666
  • 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. – Penguin47 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. – Penguin47 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. – Lutz Lehmann May 20 '19 at 14:59
  • I believe the textbook considers the "solution" to be all the values of $y(x_i)$. So for $h=\frac{1}{2}$, there is only one solution that we must find, which is $y(0.5)$. The actual answer to this is then: 0.8240271368 (can easily be found by plugging 0.5 into the equation given in the question). And the approximation is as stated in my last comment. We want to find the approximation. Thank you so much for your help so far. – Penguin47 May 20 '19 at 15:03
  • Then you should put it that way in the question: "As control value, give the numerical approximation of $y(0.5)$." or similar. Indeed that results in the value $0.8243243243243242$ while the value of the exact solution is $0.5+\frac{\sinh(1)}{\sinh(2)}=0.5+\frac{0.5}{\cosh(1)}=0.8240271368319427$. – Lutz Lehmann May 20 '19 at 15:10
  • I apologise. I wrote the question word-for-word how it was written in the textbook so I didn't think changing it would be sensible. I am now even more confused how the question can be interpreted to get different answers. Would you be able to show me how you get the solution $y(0.5) \approx 0.82432432$? Thanks. – Penguin47 May 20 '19 at 15:13
  • By performing a loop over the RK4 step with step size $h=0.5$. There is nothing mysterious about it. It should be done in vectorized form for the first-order system as that avoids easily committed errors in a component-wise implementation. // Are you expected to compute this manually? If not, what programming language / CAS system, ... did you try? – Lutz Lehmann May 20 '19 at 15:17
  • So I do a loop through RK4 for $y'(x) = z$ and for $z'(x) = 4y - 4x$ to find $y_0(x)$ and $z(x)$? What function do I use for $f(x_i, y_i)$ during RK4? Do I use $f(x_i,y_i) = 4y-4x$ and $f(x_i,y_i) = z$, respectively? Also does the final $y_{i+1} = ...$ for each give me the $y_0(x)$ and $z(x)$ I am after. (I'm sorry if you're banging your head against the wall due to my slowness here. I'm just struggling to grasp it.) I hope this makes sense. I know how to use RK4 generally, but not sure how to apply it here to get $y_0(x)$ and $z(x)$. Edit: I am doing this manually, yes. – Penguin47 May 20 '19 at 15:26
  • 1
    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)$. – Lutz Lehmann May 20 '19 at 15:33
  • How does this help me find $y_0(x)$ and $z(x)$? I'm starting to feel like a lost cause... – Penguin47 May 20 '19 at 15:52
  • 1
    You get $y_0(x)=u_1(x)$ and $z(x)=u_3(x)$, as per the construction of the state vector. – Lutz Lehmann May 20 '19 at 16:17
  • Brilliant. That makes sense. Thank you! I'm gonna have a go at it now. – Penguin47 May 20 '19 at 16:41