Prof. Vladimir A. Dobrushkin
This tutorial contains many matlab scripts.
You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to
The idea of shooting method is to reduce the given boundary value problem to several initial value problems. Roughly speaking, we 'shoot' out trajectories in different directions until we find a trajectory that has the desired boundary value. We start with the Dirichlet boundary value problem for a linear differential equation of second order:
\[
x'' (t) = p(t)\,x' + q(t)\,x + r(t) \qquad \mbox{subject} \quad u(a) = \alpha, \quad u(b) = \beta
\]
over interval [a,b]. In this case, the solution to the boundary value problem is usually given by a linear combination of two functions, u(t) and v(t) that are solutions to the initial value problems.
\[
x (t) = u(t) + \frac{\beta - u(b)}{v(b)} \, v(t) ,
\]
where u(t) is a solution to the initial value problem
\[
u'' (t) = p(t)\,u' + q(t)\,u + r(t) \qquad \mbox{subject} \quad u(a) = \alpha, \quad u'(a) = 0;
\]
and v(t) is a solution to another initial value problem:
\[
v'' (t) = p(t)\,v' + q(t)\, v , \qquad \mbox{subject} \quad v(a) = 0, \quad v'(a) = 1.
\]
Theorem. Consider the Dirichlet boundary value problem for the linear second order differential equation
\[
y'' (x) = p(x)\,y' + q(x)\, y + r(x) , \qquad \mbox{subject} \quad y(a) = A, \quad y(b) = B,
\]
where A and B are some given real numbers. The coefficients of the differential equation, p(x), q(x), and r(x) are continuous on the given interval [a,b]. If there is exists a positive constant L so
\[
\begin{split}
q(x) >: 0 \quad \mbox{for all } x \in [a,b], \\
\left\vert p(x) \right\vert < L = \max_{a \le x \le b}\, \left\{ \, |p(x) | \, \right\} .
\end{split}
\]
Then the linear boundary value problem
\[
y'' (x) = p(x)\,y' + q(x)\, y + r(x) , \qquad \mbox{subject} \quad y(a) = A, \quad y(b) = B
\]
has a unique solution \( y = \phi (x) \) over \( a \le x \le b . \)
Example. Consider a harmonic oscillator subject to the Dirichlet boundary conditions
\[
y'' + y =0 , \qquad \mbox{subject} \quad y(0) = A, \quad y(\pi ) = B.
\]
This problem has infinite many solutions because if y(x) is its solution, then z(x) = y(x) + K sin (x) is also a solution for arbitrary constant K. Therefore, the conditions of the above theorem are violated. Indeed, since \( p(x) \equiv 0 , \) there does not exist a positive constant L (it is zero).
Example. Consider the Dirichlet boundary value problem
\[
\frac{{\text d}^2 w}{{\text d} r^2} + \frac{1}{r}\,\frac{{\text d} w}{{\text d} r} - \frac{w}{r^2} =0 , \qquad \mbox{subject} \quad w(3) = 6, \quad w'(8) = 2.
\]
Since the formula for the given boundary value problem is known to be
\[
w (r) = u(r) + \frac{2 - u(8)}{v(8)} \, v(r) ,
\]
we need to find u(r), the unique solution to the initial value problem
\[
\frac{{\text d}^2 u}{{\text d} r^2} + \frac{1}{r}\,\frac{{\text d} u}{{\text d} r} - \frac{u}{r^2} =0 , \qquad \mbox{subject} \quad u(3) = 6, \quad u'(3) = 0,
\]
and v(r), which is the solution to IVP:
\[
\frac{{\text d}^2 v}{{\text d} r^2} + \frac{1}{r}\,\frac{{\text d} v}{{\text d} r} - \frac{v}{r^2} =0 , \qquad \mbox{subject} \quad v(3) = 0, \quad v'(3) = 1.
\]
The given differential equation is the Euler one, so its general solution is
\[
u(r) = A\,r + \frac{B}{r} \qquad \mbox{and} \qquad v(r) = a\, r + \frac{b}{r}.
\]
To find the values of constants, we substitute their forms into the initial conditions:
\[
\begin{split}
u(3) = A\, 3 + \frac{B}{3} =1, \quad u' (3) = A - \frac{B}{9} =0 , \\
v(3) = a\,3 + \frac{b}{3} =0, \quad a- \frac{b}{9} =1 . \end{split}
\]
Therefore,
\[
u(r) = r + \frac{9}{r} \qquad \mbox{and} \qquad v(r) = \frac{1}{2}\, r - \frac{9}{2\,r}.
\]
Since u(8) = 73/8 and v(8) = 55/16, we find the required solution
\[
w (r) = r + \frac{9}{r} + \frac{2 - 73/8}{55/16} \left( \frac{1}{2}\, r - \frac{9}{2\,r} \right) = \frac{2}{55\,r} \left( 504 - r^2 \right) ,
\]
function L=linshoot(F1,F2,a,b,alpha,beta,m)
% Input: F1 and F2 are the systems of first-order equations
% representing the IVP's
% input as strings: 'F1', 'F2'
% a and b are the left and right end points
% alpha=x(a) and beta=x(b); boundary conditions
% m is the number of steps
% Output: L=[t',x'] where t' is the (m+1)x1 vector of abscissas
% and x is the (m+1)x1 vector of ordinates
% solve the system F1
za=[alpha,0];
% call for Runhe--Kutta 4 subroutine
[t,z]=rks4(F1,a,b,za,m);
u=z(:,1);
% solve the system F2
za=[0,1];
[t,z]=rks4(F2,a,b,za,m);
v=z(:,1);
% calculate the solution to the boundary value problem
x=u+(beta-u(m+1))*v/v(m+1);
L=[t',x];