Computational Physics - Department of Physics

9.2 Shooting methods 285

We start with the Taylor expansion of the desired solution

y(x+h) =y(x)+hy(^1 )(x)+

h^2

2!

y(^2 )(x)+

h^3

3!

y(^3 )(x)+

h^4

4!

y(^4 )(x)+···

Herey(n)(x)is a shorthand notation for the nth derivativedny/dxn. Because the correspond-
ing Taylor expansion ofy(x−h)has odd powers ofhappearing with negative signs, all odd
powers cancel when we addy(x+h)andy(x−h)

y(x+h)+y(x−h) = 2 y(x)+h^2 y(^2 )(x)+

h^4

12

y(^4 )(x)+O(h^6 ).

We obtain

y(^2 )(x) =

y(x+h)+y(x−h)− 2 y(x)

h^2

−

h^2

12

y(^4 )(x)+O(h^6 ).

To eliminate the fourth-derivative term we apply the operator( 1 +h
2
12
d^2
dx^2 )on the differential
equation

y(^2 )(x)+

h^2

12

y(^4 )(x)+k^2 (x)y(x)+

h^2

12

d^2

dx^2

(

k^2 (x)y(x)

)

≈ 0.

In this expression they(^4 )terms cancel. To treat the generalxdependence ofk^2 (x)we
approximate the second derivative of(k^2 (x)y(x)by

d^2 (k^2 y(x))

dx^2

≈

(

k^2 (x+h)y(x+h)+k^2 (x)y(x)

)

+

(

k^2 (x−h)y(x−h)+k^2 (x)y(x)

)

h^2

.

We replace theny(x+h)with the shorthandyi+ 1 (and similarly for the other variables) and
obtain a final discretised algorithm for obtainingyi+ 1

yi+ 1 =

2

(

1 − 125 h^2 ki^2

)

yi−

(

1 + 121 h^2 k^2 i− 1

)

yi− 1

1 +h 122 k^2 i+ 1

+O(h^6 ),

wherexi=ih,ki=k(xi=ih)andyi=y(xi=ih)etc.
It is easy to add the termFisince we need only to take the second derivative. The final
algorithm reads then

yi+ 1 =

2

(

1 − 125 h^2 k^2 i

)

yi−

(

1 + 121 h^2 k^2 i− 1

)

yi− 1

1 +h 122 k^2 i+ 1

+

h^2

12

(Fi+ 1 +Fi− 1 − 2 Fi)+O(h^6 ).

Starting ati= 1 results in, using the boundary conditiony 0 = 0 ,

y 2 =

2

(

1 − 125 h^2 k 1 y 1

)

−

(

1 + 121 h^2 k^20 y 0

)

1 +h 122 k^22

+

h^2

12

(F 2 +F 0 − 2 F 1 )+O(h^6 ).

This equation carries a local truncation error proportional toh^6. This is an order better than
the fourth-order Runge-Kutta method which has a local errorproportional toh^5. The global
for the fourth-order Runge-Kutta is proportional toh^4 while Numerov’s method has an error
proportional toh^5. With few additional function evulations, we have achievedan increased
accuracy.
But even with an improved accuracy we end up with one unknown on the right hand side,
namelyy 1. The value ofy 1 can again be determined from the derivative aty 0 , or by a good
guess on its value. We need therefore an additional constraint on our set of equations before
we start. We could then add to the boundary conditions