Explanation File for Euler-Romberg Method ========================================= We want to solve the Cauchy's problem: | y' = f(x,y) | | with y(x0) = y0 using the Euler-Romberg method. For that, we calculate approximate values for y1, y2,...yn, the exact solution is y(x1), y(x2),...y(xn), using an iterative algorithm based on: * repeatedly apply the Euler method in the same interval with an integration step halved after each iteration. * linearly extrapolate to obtain a better approximation. This procedure allows getting the same precision all along the computation. Let be (xn, yn) and a starting step h0, the approximation yn+1 is obtained by building up a table in the same way than in the integration Romberg method. 0 We start with estimating an initial value of yn+1, noted y0 by Euler method, using step h0: 0 y0 = yn + h0 f(xn,yn) then the new step is h0/2 and we do the same computations for all the points of the sub-interval: 0 h0 y0 1 0 L=1 h0/2 y0 y1 2 2 1 0 L=2 h0/2 y0 y1 y2 ---------------------------------- k k k-1 1 0 L=k h0/2 y0 y ... y y 1 k-1 k 0 The first element y0 of the column 0 is an approximation of y(x ) with the step h0=x - x . The other elements are the n+1 n-1 n successive approximations of y(x ) given by the Euler's formula n+1 for h0/2, h0/4... i.e.: L=1 h1=h0/2 x0=xn y0=yn f0=f(x0,y0) x1=xn+h1 y1=y0+h1f0 f1=f(x1,y1) 1 x2=xn+2h1 y2=y1+h1f1 -> Y0 = Y2 ------- k The elements Y of the column m (m=1,2...) are obtained by using m the linear extrapolation formula: k k+1 k Ym-1[h-h ] - Ym-1[h-h ] k+m k Y (h) = ------------------------- m h - h k k+m and when h -> 0: k k+m k+1 k k Y [0-h0/2 ] - Y [0-h0/2 ] m-1 m-1 Y = ------------------------------- m k k+m h0/2 - h0/2 Dividing by h0/2, we obtain: m k+1 k k 2 Ym-1 - Ym-1 Y = ------------- m m 2 - 1 k Observing how the table elements Ym are calculated, we can see that it is not necessary to keep in memory a table with two dimen- k-1 k-2 0 k sions, we first calculate Y1, then Y2 ... until Yk. The table Ym is then replaced by the one-dimension table, Tk: k iteration # step Ym Tk --------------------------------------------------- 0 0 h Y0 T1 1 0 1 h/2 Y0 Y1 T2 T1 2 1 0 2 h/4 Y0 Y1 Y2 T3 T2 T1 ---- ---- --------- --------- k k 0 k h/2 Y0 ... Yk Tk ... T1 --------------------------------------------------- To fully understand the mechanism, we recommend the reader to 0 1 0 2 manually calculate a few elements: Y1, Y1 and Y2, Y1, etc. Or else: m 2 T - T K=1 k T = ------------- (k = L...1) k m 2 - 1 To sum up, we have in the program EULROMB the following steps: 1. Define the function f(x,y), the starting values, X(1), Y(1), the starting step H, the maximum error ER, the max. number of subdivisions LA and the number of calculations NC. NC = (X - X ) / H n+1 1 2. Approximate Yn+1 by Euler's method: 0 0 Y0 = Yn + h f(xn,yn) and T1 = Y0 3. Calculate the Tk elements of the extrapolation table; we initialize the iteration with L=1. a) XC = X ; YC = Y n n approximate Y by Euler's method successively dividing n+1 the step by 2 ; T is the value of Y given by L+1 n+1 Euler. b) Initialize the column count m=1; k=L is the T index. k c) calculate Ym or rather T using: k m 2 T - T k+1 k T = ------------- k m 2 - 1 d) if k<1 goto subroutine f. e) test convergence: if |Tk - Tk+1| < ER, the iteration is finished, otherwise increment m, decrement k, go to step c) f) if L