Shooting method implementation

Question

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

Answer 1

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.