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. – Wolfgang Bangerth Nov 20 '18 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}$? – nicoguaro Nov 20 '18 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. – Farid Hasanov Nov 20 '18 at 15:02
• nicoguaro, I didnt quite catch what you mean. – Farid Hasanov Nov 20 '18 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}$. – nicoguaro Nov 20 '18 at 17:49

Imagine you have the following system: $$$$\left\{\begin{array}{c} \dot{y} = f(y,z,t)\\ \dot{z} = g(y,z,t) \end{array}\right.$$$$ 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: \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. 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. – Farid Hasanov Nov 21 '18 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 '18 at 8:31
• Thanks a lot! I may come for more questions once I will implement this scheme today later : ) – Farid Hasanov Nov 21 '18 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. – LutzL Nov 21 '18 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. – Farid Hasanov Nov 22 '18 at 13:51