Mathematical Methods for Physics and Engineering : A Comprehensive Guide

27.6 DIFFERENTIAL EQUATIONS

We assume that this can be simulated by a form

yi+1=yi+α 1 hfi+α 2 hf(xi+β 1 h, yi+β 2 hfi), (27.76)

which in effect uses a weighted mean of the value ofdy/dxatxiand its value at

some point yet to be determined. The object is to choose values ofα 1 ,α 2 ,β 1 and

β 2 such that (27.76) coincides with (27.75) up to the coefficient ofh^2.

Expanding the functionfin the last term of (27.76) in a Taylor series of its

own, we obtain

f(xi+β 1 h, yi+β 2 hfi)=f(xi,yi)+β 1 h

∂fi

∂x

+β 2 hfi

∂fi

∂y

+O(h^2 ).

Putting this result into (27.76) and rearranging in powers ofh, we obtain

yi+1=yi+(α 1 +α 2 )hfi+α 2 h^2

(

β 1

∂fi

∂x

+β 2 fi

∂fi

∂y

)

. (27.77)

Comparing this with (27.75) shows that there is, in fact, some freedom remaining

in the choice of theα’s andβ’s. In terms of an arbitraryα 1 (=1),

α 2 =1−α 1 ,β 1 =β 2 =

1

2(1−α 1 )

.

One possible choice isα 1 =0.5, givingα 2 =0.5,β 1 =β 2 = 1. In this case the

procedure (equation (27.76)) can be summarised by

yi+1=yi+^12 (a 1 +a 2 ), (27.78)

where

a 1 =hf(xi,yi),

a 2 =hf(xi+h, yi+a 1 ).

Similar schemes giving higher-order accuracy inhcan be devised. Two such

schemes, given without derivation, are as follows.

(i) To orderh^3 ,

yi+1=yi+^16 (b 1 +4b 2 +b 3 ), (27.79)

where

b 1 =hf(xi,yi),

b 2 =hf(xi+^12 h, yi+^12 b 1 ),

b 3 =hf(xi+h, yi+2b 2 −b 1 ).