'******************************************************************** '* Differential equations with p variables of order 1 * '* by Runge-Kutta method of order 4 * '* ---------------------------------------------------------------- * '* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 by Marc * '* DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" * '* * '* Basic version 1.1 By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Example #1: integrate system of equations from x=0 to x=3: * '* y1' = y2 + y3 - 3*y1 * '* y2' = y1 + y3 - 3*y2 * '* y3' = y1 + y2 - 3*y3 * '* with initial conditions: y1(0)=1, y2(0)=2 and y3(0)=-1 * '* * '* * '* DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1 * '* of type yi' = f(y1,y2,...,yn), i=1..n * '* * '* number of variables: 3 * '* begin value x : 0 * '* end value x : 3 * '* y1 value at x0 : 1 * '* y2 value at x0 : 2 * '* y3 value at x0 : -1 * '* number of points : 7 * '* finesse : 30 * '* * '* X Y1 Y2 Y3 * '* -------------------------------------------------------- * '* 0.000000 1.000000 2.000000 -1.000000 * '* 0.500000 0.449466 0.584801 0.178795 * '* 1.000000 0.251358 0.269674 0.214727 * '* 1.500000 0.149580 0.152058 0.144622 * '* 2.000000 0.090335 0.090671 0.089664 * '* 2.500000 0.054738 0.054784 0.054648 * '* 3.000000 0.033193 0.033200 0.033181 * '* -------------------------------------------------------- * '* * '* Example #2: integrate system of equations from x=0 to PI*Sqrt(2) * '* y1' = y2 * '* y2' = -4*y1 - 3*y3 * '* y3' = y4 * '* y4' = -8*y1 - 2*y3 * '* with initial conditions: y1(0)=3, y2(0)=0, y3(0)=4, y4(0)=0 * '* * '* * '* DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1 * '* of type yi' = f(y1,y2,...,yn), i=1..n * '* * '* number of variables: 4 * '* begin value x : 0 * '* end value x : 4.442883 * '* y1 value at x0 : 3 * '* y2 value at x0 : 0 * '* y3 value at x0 : 4 * '* y4 value at x0 : 0 * '* number of points : 9 * '* finesse : 30 * '* * '* X Y1 Y2 Y3 Y4 * '* -------------------------------------------------------- * '* 0.000000 3.000000 0.000000 4.000000 0.000000 * '* 0.555360 0.000000 -8.485281 0.000000 -11.313708 * '* 1.110721 -3.000000 -0.000001 -4.000000 -0.000002 * '* 1.666081 -0.000001 8.485281 -0.000001 11.313708 * '* 2.221442 3.000000 0.000003 4.000000 0.000003 * '* 2.776802 0.000001 -8.485281 0.000002 -11.313708 * '* 3.332162 -3.000000 -0.000004 -4.000000 -0.000005 * '* 3.887523 -0.000002 8.485281 -0.000002 11.313708 * '* 4.442883 3.000000 0.000005 4.000000 0.000007 * '* -------------------------------------------------------- * '* * '* Release 1.1: Added example #2 (Sept. 2008). * '******************************************************************** DefInt I-N DefDbl A-H, O-Z 'ifi,ip : INTEGER Dim y(10), yi(10), t1(50), t2(50), t3(50), t4(50) Cls Print Print " DIFFERENTIAL EQUATIONS WITH P VARIABLES OF ORDER 1" Print " of type yi' = f(y1,y2,...,yn), i=1..n" Print PRINT " number of variables : "; : INPUT ip Print PRINT " begin value x : "; : INPUT xi PRINT " end value x : "; : INPUT xf For i = 0 To ip - 1 PRINT " y"; i + 1; " value at x0 : "; : INPUT yi(i) Next i PRINT " number of points : "; : INPUT m PRINT " finesse : "; : INPUT ifi 'call subroutine eqdifp GoSub 2000 End 'Example #1: y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3 '1000 'FUNCTION fp ' IF k = 0 THEN ' fp = y(1) + y(2) - 3# * y(0) ' ELSEIF k = 1 THEN ' fp = y(0) + y(2) - 3# * y(1) ' ELSEIF k = 2 THEN ' fp = y(0) + y(1) - 3# * y(2) ' ELSE ' fp = 0# ' END IF 'RETURN 'Example #2: y1'=y2, y2'=-4y1-3y3, y'3=y4, y'4=-8y1-2y3 1000 'FUNCTION fp If k = 0 Then fp = y(1) ElseIf k = 1 Then fp = -4# * y(0) - 3# * y(2) ElseIf k = 2 Then fp = y(3) ElseIf k = 3 Then fp = -8# * y(0) - 2# * y(2) Else fp = 0# End If Return '**************************************************************************** '* SOLVING DIFFERENTIAL SYSTEMS WITH P VARIABLES OF ORDER 1 * '* of type yi' = f(y1,y2,...,yn), i=1..n * '* ------------------------------------------------------------------------ * '* INPUTS: * '* m number of points to calculate * '* xi, xf begin, end values of variable x * '* yi table of begin values of functions at xi * '* ip number of independant variables * '* ifi finesse (number of intermediary points) * '* * '* OUTPUTS: * '* t1,t2 real vectors storing the results for first two functions, * '* y1 and y2. * '* ------------------------------------------------------------------------ * '* EXAMPLE: y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3 * '* Exact solution : y1 = 1/3 (exp(-4x) + 2 exp(-x)) * '* y2 = 1/3 (4exp(-4x) + 2 exp(-x)) * '* y3 = 1/3 (-5exp(-4x)+ 2 exp(-x)) * '**************************************************************************** 2000 'subroutine Eqdifp Dim ta(10), tb(10), tc(10), td(10), z(10) h = (xf - xi) / ifi / (m - 1) ip = ip - 1 t1(1) = yi(0) t2(1) = yi(1) t3(1) = yi(2) t4(1) = yi(3) For k = 0 To ip y(k) = yi(k): z(k) = yi(k) Next k For i = 1 To m ni = (i - 1) * ifi - 1 For j = 1 To ifi X = xi + h * (ni + j) For k = 0 To ip y(k) = z(k) Next k For k = 0 To ip GoSub 1000 ta(k) = h * fp Next k For k = 0 To ip y(k) = z(k) + ta(k) / 2# Next k X = X + h / 2# For k = 0 To ip GoSub 1000 tb(k) = h * fp Next k For k = 0 To ip y(k) = z(k) + tb(k) / 2# Next k For k = 0 To ip GoSub 1000 tc(k) = h * fp Next k For k = 0 To ip y(k) = z(k) + tc(k) Next k X = X + h / 2# For k = 0 To ip GoSub 1000 td(k) = h * fp Next k For k = 0 To ip z(k) = z(k) + (ta(k) + 2# * tb(k) + 2# * tc(k) + td(k)) / 6# Next k Next j t1(i + 1) = z(0) t2(i + 1) = z(1) t3(i + 1) = z(2) t4(i + 1) = z(3) Next i 'Display(t1,t2,t3,t4,m,p,xi,xf) GoSub 3000 Return 3000 'subroutine Display h = (xf - xi) / (m - 1) X = xi - h Print Print " X"; For i = 1 To ip + 1 Print " Y"; i; Next i Print Print "--------------------------------------------------------" For i = 1 To m X = X + h If ip = 1 Then Print USING; " ##.###### ##.###### ##.######"; X; t1(i); t2(i) ElseIf ip = 2 Then Print USING; " ##.###### ##.###### ##.###### ##.######"; X; t1(i); t2(i); t3(i) ElseIf ip > 2 Then Print USING; " ##.###### ##.###### ##.###### ##.###### ###.######"; X; t1(i); t2(i); t3(i); t4(i) End If Next i Print "--------------------------------------------------------" Return 'End of file teqdifp.bas
Wednesday, September 20, 2017
Runge-Kutta method of order 4
Monday, September 18, 2017
Adams-Moulton Prediction-Correction Method
'******************************************************************** '* Solve Y' = F(X,Y) with initial conditions using the Adams- * '* Moulton Prediction-Correction Method * '* ---------------------------------------------------------------- * '* SAMPLE RUN * '* (Integrate Y' = -Y + X/(1+X)^2 from X=0 to X=1 with initial * '* condition Y(0) = 1 ) * '* * '* X Y Y True |Y-Y True| * '* --------------------------------------------- * '* 0.000000 1.000000 1.000000 0.000000 * '* 0.050000 0.952381 0.952381 0.000000 * '* 0.100000 0.909091 0.909091 0.000000 * '* 0.150000 0.869569 0.869565 0.000004 1 * '* 0.200000 0.833340 0.833333 0.000006 1 * '* 0.250000 0.800008 0.800000 0.000009 1 * '* 0.300000 0.769241 0.769231 0.000010 1 * '* 0.350000 0.740752 0.740741 0.000011 1 * '* 0.400000 0.714298 0.714286 0.000012 1 * '* 0.450000 0.689668 0.689655 0.000012 1 * '* 0.500000 0.666679 0.666667 0.000013 1 * '* 0.550000 0.645174 0.645161 0.000013 1 * '* 0.600000 0.625013 0.625000 0.000013 1 * '* 0.650000 0.606073 0.606061 0.000013 1 * '* 0.700000 0.588248 0.588235 0.000013 1 * '* 0.750000 0.571441 0.571429 0.000013 1 * '* 0.800000 0.555568 0.555556 0.000012 1 * '* 0.850000 0.540553 0.540541 0.000012 1 * '* 0.900000 0.526327 0.526316 0.000012 1 * '* 0.950000 0.512832 0.512821 0.000011 1 * '* 1.000000 0.500011 0.500000 0.000011 1 * '* ---------------------------------------------------------------- * '* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en * '* Basic et en Pascal By Claude Nowakowski, Edition du * '* P.S.I., 1984". * '* * '* Quick Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '******************************************************************** 'See explanation file adambash.txt '--------------------------------- DefInt I-N DefDbl A-H, O-Z Option Base 0 'index begins from zero Dim X(3), Y(3) H = 0.05 'integration step X(0) = 0#: Y(0) = 1# 'Initial conditions EC = 0.000001 'Precision F$ = " ##.###### ##.###### ##.###### ##.###### ##" Cls Print " X Y Y True |Y-Y True| " Print " ---------------------------------------------" Print USING; F$; X(0); Y(0); Y(0); X(0) 'Start with Runge-Kutta For K = 0 To 1 XX = X(K): YY = Y(K): GoSub 1000: C1 = C XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C1: GoSub 1000: C2 = C YY = Y(K) + H / 2# * C2: GoSub 1000: C3 = C X(K + 1) = X(K) + H XX = X(K + 1): YY = Y(K) + H * C3: GoSub 1000: C4 = C Y(K + 1) = Y(K) + H * (C1 + 2 * C2 + 2 * C3 + C4) / 6# XX = X(K + 1): GoSub 2000: ER = Abs(F - Y(K + 1)) Print USING; F$; X(K + 1); Y(K + 1); F; ER Next K 100 K = 2 XX = X(K): YY = Y(K): GoSub 1000: C1 = C XX = X(K - 1): YY = Y(K - 1): GoSub 1000: C2 = C XX = X(K - 2): YY = Y(K - 2): GoSub 1000: C3 = C X(K + 1) = X(K) + H YP = Y(K) + H / 12# * (23 * C1 - 16 * C2 + 5 * C3) L = 0 200 XX = X(K + 1): YY = YP: GoSub 1000: C1 = C XX = X(K): YY = Y(K): GoSub 1000: C2 = C XX = X(K - 1): YY = Y(K - 1): GoSub 1000: C3 = C YC = Y(K) + H / 12# * (5 * C1 + 8 * C2 - C3) 'PRINT YC If Abs(YP - YC) > EC Then YP = YC: L = L + 1: GoTo 200 End If Y(K + 1) = YC: XX = X(K + 1): GoSub 2000: ER = Abs(F - Y(K + 1)) Print USING; F$; X(K + 1); Y(K + 1); F; ER; L For K = 0 To 2 X(K) = X(K + 1): Y(K) = Y(K + 1) Next K If X(2) < 1# Then GoTo 100 End If End 1000 C = -YY + XX / ((1# + XX) ^ 2) Return 2000 F = 1# / (1# + XX) Return
Monday, September 11, 2017
Adams-Bashforth Method
'******************************************************************** '* Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the * '* Adams-Bashforth Method * '* ---------------------------------------------------------------- * '* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en * '* Basic et en Pascal By Claude Nowakowski, Edition du * '* P.S.I., 1984" [4]. * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* (Solve Y' = -Y + X/((1+X)*(1+X)) with Y(0)=1 for X0=0 up to * '* X1=1.0, exact solution is Y = 1 / (1+X). * '* * '* X Y Y exact Error * '* ----------------------------------------------- * '* 0.000000 1.000000 1.000000 0 * '* 0.000000 1.000000 1.000000 0.000000 * '* 0.050000 0.952381 0.952381 0.000000 * '* 0.100000 0.909091 0.909091 0.000000 * '* 0.150000 0.869525 0.869565 0.000040 * '* 0.200000 0.833265 0.833333 0.000068 * '* 0.250000 0.799910 0.800000 0.000090 * '* 0.300000 0.769125 0.769231 0.000106 * '* 0.350000 0.740623 0.740741 0.000117 * '* 0.400000 0.714160 0.714286 0.000125 * '* 0.450000 0.689525 0.689655 0.000131 * '* 0.500000 0.666533 0.666667 0.000134 * '* 0.550000 0.645026 0.645161 0.000135 * '* 0.600000 0.624865 0.625000 0.000135 * '* 0.650000 0.605926 0.606061 0.000134 * '* 0.700000 0.588103 0.588235 0.000133 * '* 0.750000 0.571298 0.571429 0.000131 * '* 0.800000 0.555428 0.555556 0.000128 * '* 0.850000 0.540416 0.540541 0.000125 * '* 0.900000 0.526194 0.526316 0.000121 * '* 0.950000 0.512703 0.512821 0.000118 * '* 1.000000 0.499886 0.500000 0.000114 * '* ----------------------------------------------- * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '******************************************************************** DefDbl A-H, O-Z DefInt I-N Option Base 0 H = 0.05 'integration step Dim B(4), X(4), Y(4) 'Initial conditions X(0) = 0# 'starting X X1 = 1# 'ending X Y(0) = 1# 'initial Y 'write header Cls Print Print " X Y Y exact Error " Print " ------------------------------------------------" 'write initial line F$ = "#####.######" Print USING; F$; X(0); Print USING; F$; Y(0); Print USING; F$; Y(0); Print " 0" 'use Runge-Kutta to start For K = 0 To 1 GoSub 1000 'call RK4 XX = X(K + 1): GoSub 600: ER = Abs(FX - Y(K + 1)) Print USING; F$; X(K + 1); Print USING; F$; Y(K + 1); Print USING; F$; FX; Print USING; F$; ER Next K 'main integration loop While X(2) < X1 For I = 1 To 3 XX = X(3 - I): YY = Y(3 - I): GoSub 500 B(I) = F Next I X(3) = X(2) + H Y(3) = Y(2) + H * (23# * B(1) - 16# * B(2) + 5# * B(3)) / 12# XX = X(3): GoSub 600: ER = Abs(Y(3) - FX) Print USING; F$; X(3); Print USING; F$; Y(3); Print USING; F$; FX; Print USING; F$; ER For K = 0 To 2 X(K) = X(K + 1): Y(K) = Y(K + 1) Next K Wend Print " ------------------------------------------------" End 'User defined function Y'=F(XX,YY) 500 F = -YY + XX / ((1# + XX) * (1# + XX)) Return 'Exact solution Y=FX(XX) 600 FX = 1# / (1# + XX) Return 'Runge-Kutta method to calculate first points only 1000 XX = X(K): YY = Y(K): GoSub 500: C1 = F XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C1: GoSub 500: C2 = F XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C2: GoSub 500: C3 = F XX = X(K) + H: YY = Y(K) + H * C3: GoSub 500: C4 = F X(K + 1) = X(K) + H Y(K + 1) = Y(K) + H * (C1 + 2 * C2 + 2 * C3 + C4) / 6# Return 'end of file adambash.bas
EXPLANATION FILE OF PROGRAM ADAMBASH
====================================
Ordinary Differential Equations Y' = F(x,y)
-------------------------------------------
Linked Steps Method
-------------------
We have seen that Euler, Runge-Kutta... methods can be put under the general
form:
y = y + h phi(x , y , h)
n+1 n n n
and each point y of the solution is only determined from previous point, y
n+1 n
Calculations are made for y1, y2,...,yn-1, yn, yn+1... These methods are cal-
led "with separate steps": each step is independant from the previous one. To
obtain a given accuracy, each step has to be divided into intermediate steps.
Let us take an example:
For the Runge-Kutta method of order 4, we use the formulas:
k1 = h f(xn, yn)
k2 = h f(xn+h/2, yn+k1/2)
k3 = h f(xn+h/2, yn+k2/2)
k4 = h f(xn+h, yn+k3)
y = y + (1/6)(k1+2k2+2k3+k4)
n+1 n
Another method consists in using the previous calculations to improve speed;
the y point is evaluated from y , y and y points.
n+1 n n-1 n-2
The general formula is:
y = a y + a y + ... + a y
n+1 0 n 1 n-1 k n-k
+ h (b f + b f + b f + ... + b f
-1 n+1 0 n 1 n-1 k n-k
This can be written as:
k k
y = Sum a y + h Sum b f
n+1 j=0 j n-j j=-1 j n-j
These methods are called "with linked steps" and f(x,y) evaluations are only
made at points x0, x1, x2,..., xn.
If b = 0, the process is explicit; y is directly obtained by applying the
-1 n+1
formula.
If b <> 0, the process is implicit: we must solve an equation of the form
-1
y = phi(y ) to obtain y .
n+1 n+1 n+1
For the three first points, y1, y2, y3, we have no previous points to calcu-
late them: these methods cannot start by themselves like "with separate steps"
methods.
The truncation error Tn is estimated by using the Taylor formula for y(x )
n+j
(jh)² (jh)^p (p)
y(x ) = y(x ) + jhy'(x ) + ----- y"(x ) + ... + ------ y (x )
n+j n n 2! n p! n
+ h^p eps(h)
and also for f(x , y(x )):
n+j n+j
f(x , y(x )) = y'(x +jh)
n+j n+j n p-1
h (jh) (p)
= y'(x ) + j y"(x ) + ... + -------- y (x ) + h^p eps(h)
n n (p-1)! n
Example: let us estimate the truncation error for the implicit formula:
y = y + (h/2) (f + f )
n+1 n n+1 n
T = y(x ) - y
n+1 n+1 n+1
3
y(x ) = y(x ) + hy'(x ) + (h^2/2) y"(x ) + (h^3/3) y"'(x ) + h eps(h)
n+1 n n n n
2 3
(h/2) f = (h/2) y'(x + h) = (h/2) y'(x ) + (h /2) y"(x ) + (h /4)
n+1 n n n
3
y"'(x ) + h eps(h)
n
(h/2) f = (h/2) y'(x )
n n
Hence y(x ) - y = y(x ) - y + (h/2) (f + f ) =
n+1 n+1 n+1 n n+1 n
3 3
= (h /12) y"'(x ) + h eps(h)
n
Here the method is of order two.
The methods with separate steps can be integrated by the formula
xn+2h
y(x ) - y(x ) = Sum f(t,y(t)) dt
n+1 n-M xn
and by applying the Simpson's formula
b b - a a + h
Sum (f(t) dt = ----- {f(x) + 4f(-----) + f(b)}
a 6 2
Hence y - y = (h/3)[f + 4f + f ]
n+2 n n n+1 n+2
Here we have an implicit process.
In a more general way, knowing y , y , y , we can calculate f , f ,
n n-1 n-2 n n-1
f and approximate y' = f(x,y) by an interpolation polynomial at points
n-2
x , x , x :
n n-1 n-2 x
n+1
y = y + Sum P(x) dx
n+1 n-M x
n-M
where y is an approximation of y(x ) and y is an approximation of
n-M n-M n+1
y(x ).
n+1
In the case of an implicit method, y(x ) can be approximated by the formula:
n+1
x
p n+1
y = y + Sum P(x) dx
n+1 n-M x
n-M
p p
This allows evaluating f = f(x , y ) and we can resume the interpolation
n+1 n+1 n+1
step with a new polynomial P*(x). y is calculated by the correction formula
n+1
x
c n+1
y = y + Sum P*(x) dx
n+1 n-M x
n-M
As the points are equally spaced and indices are decreasing, xn, xn-1, xn-2...
we can use the Newton formula with back diferences:
Div (f ) = f - f
n n n-1
2
Div (f ) = f - 2f +f
n n n-1 n-2
--------------------------
k |k| |k| k
Div (f ) = f - | | f | | f +...+ (-1) f
n n |1| n-1 |2| n-2 n-k
Hence x 2
p n+1 Div(fn) Div (fn)
y = y + Sum [f + ------- (x-x ) + -------- (x-x )(x-x ) + ...
n+1 n-M x n h n 2h^2 n n-1
n-M
k
Div (fn)
+ -------- (x-x )(x-x )...(x-x )] dx
k! h^k n n-1 n-k+1
Let us put u = (x-x )/h, then du = dx/dh and
n 2
p 1 Div (fn)
y = y + Sum [f + Div(f )u + -------- u(u+1) + ...
n+1 n-M -M n n 2
k
Div (fn)
+ -------- u(u+1)(u+2)...(u+k-1) h du
k!
After integration:
p _ _ _ 2 _ k
y = y + h [P f + P Div(f ) + P Div (f ) + ... + P Div (f )]
n+1 n-M 0 n 1 n 2 n k n
_ 1 1
where P = -- Sum u(u+1)(u+2)...(u+j-1) du
j j! -M
j |j| |j| j
and Div (f ) = f - | | f | | f +...+ (-1) f
n n |1| n-1 |2| n-2 n-j
|j| j!
with | | = --------- (Newton's coefficients)
|m| m! (j-m)!
Fianally p
y = y
n+1 n-M + h [P f + P f + ... + P f ]
0 n 1 n-1 k n-k
Example: M=3, k=2
x 2
p n+1 Div(fn) Div (fn)
y = y + Sum [f + ------- (x-x ) + -------- (x-x )(x-x )] dx
n+1 n-3 x n n n 2h² n n-1
n-3
Let us put u = (x-xn)/h, then
2
p 1 Div (fn)
y = y + Sum [f + Div(f ) u + -------- u (u+1)] h du
n+1 n-3 -3 n n 2
After integration:
p 2
y = y + h [4f - 4 Div(f ) + (8/3)Div (f )
n+1 n-3 n n n
This leads to the Milne's formula (explicit process):
p
y = y + (4h/3) [2f - f + 2f ]
n+1 n-3 n n-1 n-2
We do in a similar way for the implicit process. For example, with M=1 and
k=2, the Milne's corrector formula is:
c
y = y + (h/3) [f + 4f + f ]
n+1 n-1 n+1 n n-1
This corresponds to the numerical integration of
x
n+1
Sum f(x) dx by the Simpson's method.
x
n-1
For the Adams-Moulton's implicit formulas, we have:
p h
k=1 y = y + - [3f - f ]
n+1 n 2 n n-1
c h
y = y + - [f + f ]
n+1 n 2 n+1 n
p h
k=2 y = y + -- [23f - 16f + 5f ]
n+1 n 12 n n-1 n-2
c h
y = y + -- [5f + 8f - f ]
n+1 n 12 n+1 n n-1
p h
k=3 y = y + -- [55f - 59f + 17f - 9f ]
n+1 n 24 n n-1 n-2 n-3
c h
y = y + -- [9f + 19f - 5f +f ]
n+1 n 24 n+1 n n-1 n-2
p h
k=4 y = y + --- [1901f - 2984f + 2616f - 1274f + 251f ]
n+1 n 720 n n-1 n-2 n-3 n-4
c h
y = y + --- [251f + 646f - 264f +106f - 19f ]
n+1 n 720 n+1 n n-1 n-2 n-3
In these formulas, the main part of what is left aside in the integration
corresponds to the truncation error.
3
c h "'
Example: for y = y + (h/2) (f + f ) - -- y (ksi)
n+1 n n+1 n 12
with x <= ksi <= x
n n+1
Another way to obtain an explicit or implicit process is using a recursive
formula:
y = a y + a y + ... + a y
n+1 0 n 1 n-1 k n-k
+ h [b f + b f + ... + b y ]
-1 n+1 0 n k n-k
For y, the polynomial of greatest degree is:
y(x) = 1, y'(x) = 0
y(x) = x, y'(x) = 1
y(x) = x², y'(x) = 2x
---------------------
m m-1
y(x) = x , y'(x) = mx
So
1 = a0 + a1 + ... + ak
x = (x-h)a + (x-2h)a + ... + ha + b +b + b + ... +b
0 1 k-1 -1 0 1 k
x² = (x-h)²a + (x-2h)²a + ... + ha
0 1 k-1
+ 2 [xb +(x-h)b + ... + hb
-1 0 k-1
--------------------------------------------------
l l l
x = (x-h) a + (x-2h) a + ... + ha
0 1 k-1
l-1 l-1
+ l [x b +(x-h) b + ... + hb
-1 0 k-1
Example: Nystroem's explicit formula:
3
y = y + h Sum b f
n+1 n-1 j=0 j n-j
3 2
y = x , y' = 3x
3 3
==> y = x , y = (x-2h)
n+1 n-1
2 2 2
f = y' = 3(x-h) , f = 3(x-2h) , f = 3(x-3h)
n n n-1 n-2
3 2 2 2 2
==> x = (x-2h) + 3h [b (x-h) + b (x-2h) + b (x-3h)
0 1 2
By developing and identifying:
b0 + b1 + b2 = 2
b0 + 2b1 +3b2 = 2
b0 + 4b1 + 9b2 = 8/3
We find: b0 = 7, b1 = -2, b2 = 1.
So we have the Adams's formulas:
Explicit, order 2: y = y + h/2 (3f - f )
n+1 n n n-1
order 3: y = y + h/12 (23f - 16f + 5f )
n+1 n n n-1 n-2
Implicit, order 2: y = y + h/2 (f + f )
n+1 n n+1 n
and Nystroem's formulas:
Explicit, order 2: y = y + 2h f
n+1 n-1 n
order 3: y = y + h/3 (7f - 2f + f )
n+1 n+1 n n-1 n-2
Implicit, order 3: y = y + h/12 (f + 4f -f )
n+1 n-1 n+1 n n-1
From [BIBLI 04].
------------------------------------------------
End of file adambash.txt
Thursday, September 7, 2017
Euler-Romberg Method
'******************************************************************** '* Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the * '* Euler-Romberg Method * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* (Solve Y' = X*X*Y with Y(0)=1 for X0=0 up to Xn=1.1, exact * '* solution is Y = Exp(X^3/3) ). * '* * '* X Y Y Error Number of * '* estimated exact subdivisions * '* ---------------------------------------------------------- * '* 0.1 1.000333 1.000333 0.00000000 4 * '* 0.2 1.002670 1.002670 0.00000001 4 * '* 0.3 1.009041 1.009041 0.00000006 4 * '* 0.4 1.021562 1.021562 0.00000014 4 * '* 0.5 1.042547 1.042547 0.00000027 4 * '* 0.6 1.074654 1.074655 0.00000086 4 * '* 0.7 1.121125 1.121126 0.00000107 4 * '* 0.8 1.186094 1.186095 0.00000126 4 * '* 0.9 1.275067 1.275069 0.00000133 4 * '* 1.0 1.395611 1.395612 0.00000114 4 * '* 1.1 1.558410 1.558412 0.00000047 4 * '* ---------------------------------------------------------- * '* * '* Ref.: "Methodes de calcul numerique By Claude Nowakowski, Tome 2 * '* PSI Editions, France, 1981" [BIBLI 04]. * '* * '******************************************************************** 'Program EulerRomberg DefInt I-N DefDbl A-H, O-Z Option Base 0 NMAX = 100 H = 0.1 'initial integration step ER = 0.000001 'desired precision LA = 10 'maximum number of subdivisions NC = 10 'number of calculations NC = (Xn+1 - X1)/H Dim X(NMAX), Y(NMAX), T(20) 'Initial conditions X(0) = 0#: Y(0) = 1# 'write header Cls Print Print " X Y Y Error Number of " Print " estimated exact subdivisions " Print "---------------------------------------------------------" 'main integration loop For N = 0 To NC XC = X(N): YC = Y(N) XX = XC: YY = YC: GoSub 1000 T(1) = Y(N) + H * F L = 1: LM = 2: ET = 1# While L < LA And ET >= ER XC = X(N): YC = Y(N) For J = 1 To LM XC = XC + H / LM XX = XC: YY = YC: GoSub 1000 YC = YC + H / LM * F Next J T(L + 1) = YC: M = 1: K = L: MM = 2: ET = 1# If K > 1 Then While ET >= ER And K > 1 T(K) = (MM * T(K + 1) - T(K)) / (MM - 1) ET = Abs(T(K) - T(K - 1)) M = M + 1: K = K - 1: MM = MM * 2 Wend End If If K = 1 Then L = L + 1: LM = LM * 2 End If Wend X(N + 1) = X(N) + H: Y(N + 1) = T(K) XX = X(N + 1): GoSub 2000: YEX = FX EF = Abs(YEX - Y(N + 1)) Print USING; "###.#"; X(N + 1); Print USING; "#####.######"; Y(N + 1); Print USING; "#####.######"; YEX; Print USING; "#####.########"; EF; Print " "; L Next N Print "---------------------------------------------------------" End 'of main program 'Y' = F(X,Y) 1000 'Function F(XX,YY) F = XX * XX * YY Return 'Exact solution FX(XX) 2000 'Function FX(XX) FX = Exp(XX ^ 3 / 3) Return 'end of file eulromb.bas
Wednesday, September 6, 2017
Explanation File for Euler-Romberg Method
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<Lmax, increment L then goto step a) else display a
message "No convergence!"
g) if n < N, increment n, go to 2. else end program.
Note: the program EULROMB is made such as the solution Yk is
stored for all the points (this is not always necessary).
[From BIBLI 04].
SOURCE WEBSITE: http://jean-pierre.moreau.pagesperso-orange.fr/eqdiff.html
Saturday, September 2, 2017
VB6: Waking a Sleeping Giant
Don't Get Me Started - VB6: Waking a Sleeping Giant:
https://msdn.microsoft.com/en-us/magazine/mt632280
Subscribe to:
Posts (Atom)