# Shooting method implementation

I have a trouble of defining a Jacobian matrix for my problem. Basically, I have 4 differential equations to be solved.

\begin{aligned} \dot x_1(t)&=x_2(t)\\ \dot x_2(t)&=p_2(t)−\sqrt 2 x_1(t)e^{-αt}\\ \dot p_1(t)&=\sqrt 2p_2(t)e^{-αt}+x_1(t)\\ \dot p_2(t)&=−p_1(t) \end{aligned}

with initial and boundary values of: $$x_1(0)=1\\ p_2(0)=0\\p_1(1)=0\\p_2(1)=0$$

Following classic shooting method strategy, I suggest some values for $$x_2(0), p_1(0)$$, solve the system of diff equations using Dormand Prince method and obtain some values at the boundary.

Here I check, if the solution on my boundaries meets the criteria of $$p_1(1)=0,p_2(1)=0$$. If not, I have to suggest new values for the initial values of $$x_2(0), p_1(0)$$ to be iterated again. Since this is a two-parametric shooting, I can't use simple linear interpolation as in one-parametric shooting, so I have to use Newton's method for that.

Here is my question: How to implement Newton's method to iterate towards the right solution? I tried myself and it doesn't work. Can you at least help me with defining a Jacobian, for the given problem of iterative solution to $$x_2(0), p_1(0)$$?

• What literature have you read in this regard? The shooting method has been used for decades. It is well documented how to do what you want to do. Nov 20, 2018 at 14:01
• If you write your equation as $\dot{\mathbf{X}} = \mathbf{F}(\mathbf{X}, t)$, do you know how to compute the Jacobian of $\mathbf{F}$? Nov 20, 2018 at 14:37
• Wolfgang, in the literature they explain one-parametric shooting , so the next iterative value for shooting is determined by linear interpolation, which is straightforward. I have two-parametric shooting , so I need to go with Newton's method for the next iterative value. However, it seems that my Jacobian matrix consists of only two values( [$\dfrac{df}{dx_2}$, $\dfrac{df}{dp_1}$]) so I don't understand how to compute Jacobian itself. Nov 20, 2018 at 15:02
• nicoguaro, I didnt quite catch what you mean. Nov 20, 2018 at 15:15
• Please use the @ before the handle to notify people you commented. The Jacobian is the matrix with the partial derivatives of the function $\mathbf{F}$. Nov 20, 2018 at 17:49

Imagine you have the following system: $$\begin{equation} \left\{\begin{array}{c} \dot{y} = f(y,z,t)\\ \dot{z} = g(y,z,t) \end{array}\right. \end{equation}$$ being $$y$$ and $$z$$ vectors with appropiate BC $$y(0) = y_0$$ and $$z(T) = z_T$$ and $$f$$ and $$g$$ vector valued functions.

Since you don't have any initial condition for $$z$$, you have to solve the nonlinear equation: $$F_1(z_0) = \dot{z}(z_0,t)-g(y,z(z_0,t),t)=0$$ subjected to $$F_2(z_0) = z(z_0,T) - z_T$$

If you linearise last equations w.r.t. $$z_0$$ you obtain a system for $$\delta z_0$$ introducing a new variable: $$\zeta(t) =\frac{\partial z(z_0,t)}{\partial z_0}$$ Which is a matrix if $$z$$ has more than 1 component.

The equations are the following: $$F_1(z_0+\delta z_0)=F_1(z_0)+\frac{dF_1(z_0)}{dz_0}\delta z_0=\frac{dF_1(z_0)}{dz_0}\delta z_0=\left(\dot{\zeta}-\frac{\partial g}{\partial z}\zeta\right)\delta z_0=0$$ Since we solve $$F_1(z_0) = 0$$ exactly with a guessed value $$z_0$$.

We note that to solve the first linearised equation we need an initial value for $$\zeta$$, and according to the definition, in is exactly the identity matrix: $$\zeta(0) = \frac{\partial z(z_0,0)}{\partial z_0}=\frac{\partial z_0}{\partial z_0}=I$$ $$F_2(z_0+\delta z_0)=F_2(z_0)+\frac{dF_2(z_0)}{dz_0}\delta z_0=z(T,z_0)-z_T+\zeta(T)\delta z_0=0$$

Final equations

Therefore, you have to solve iteratively for $$z_0$$ the following equations: \begin{equation} \left\{\begin{aligned} \dot{y} &= f(y,z,t)\\ \dot{z} &= g(y,z,t)\\ \dot{\zeta} &= \frac{\partial g}{\partial z}\zeta \end{aligned}\right. \end{equation} With the initial conditions $$y(0) = y_0$$, $$z(0) = z_0^{n}$$, $$\zeta(0) = I$$, and the updating equations for $$z_0$$: $$z_0^{n+1} = z_0^n + \delta z_0\qquad \zeta(T)\delta z_0=z_T-z(T,z_0)$$

with $$z_0^0$$ a value close enough to the actual one.

• sorry, I didnt get this part: you linearise last equations w.r.t. z0 you obtain a system for δz0 introducing a new variable:. Can you explain it more detailed ? Anyway , thanks a lot. Nov 21, 2018 at 8:10
• You have a nonlinear system of equations. In order to solve it, you must go for an iterative method (for example newton's method) and for this reason you must linearise your system according to $z_0\rightarrow z_0+\delta z_0$ supposing $\delta z_0\ll 1$
– HBR
Nov 21, 2018 at 8:31
• Thanks a lot! I may come for more questions once I will implement this scheme today later : ) Nov 21, 2018 at 11:05
• Your solution is incomplete, as also $y=y(t; z_0)$ depends on $z_0$ and thus you need to include the derivatives of $f$ in the calculation of the Jacobian. Note also that the structure of the boundary conditions in the given problem prevents a separation of the variables in the form you used. Nov 21, 2018 at 22:38
• I am thinking of writing my Jacobian either like this $$\begin{bmatrix} \dfrac{dx_2(1)}{x_2(0)} & \dfrac{dx_2(1)}{p_1(0)} \\ \dfrac{dp_1(1)}{x_2(0)} & \dfrac{dp_1(1)}{p_1(0)} \end{bmatrix}$$ or like this $$\begin{bmatrix} \dfrac{dp_1(1)}{x_2(0)} & \dfrac{dp_1(1)}{p_1(0)} \\ \dfrac{dp_2(1)}{x_2(0)} & \dfrac{dp_2(1)}{p_1(0)} \end{bmatrix}$$ In the first case the functions are the product of the evaluation, second case describes functions of residuals. I think the first one is more correct. Nov 22, 2018 at 13:51