'*************************************************************
'* Differential equation y"=f(x,y,y') by Stormer's method *
'* --------------------------------------------------------- *
'* SAMPLE RUN: *
'* Integrate y" = 8yy / (1 + 2x) from x=0 to x=1, *
'* with initial conditions: x(0)=0, y(0)=1 and y'(0)=-2 *
'* and compare with exact solution: y = 1 / (1 + 2x) *
'* *
'* Output file (stormer.lst): *
'* *
'* --------------------------------------------------------- *
'* Differential equation y"=f(x,y,y') by Stormer's method *
'* --------------------------------------------------------- *
'* X Y Y exact Error *
'* 0.000 1.000000 1.000000 0.0000000000 *
'* 0.010 0.980392 0.980392 0.0000000001 *
'* 0.020 0.961538 0.961538 0.0000000295 *
'* 0.030 0.943396 0.943396 0.0000000457 *
'* 0.040 0.925926 0.925926 0.0000000974 *
'* 0.050 0.909091 0.909091 0.0000001285 *
'* ... ... ... ... *
'* 0.950 0.344866 0.344828 0.0000381695 *
'* 0.960 0.342505 0.342466 0.0000388874 *
'* 0.970 0.340176 0.340136 0.0000396196 *
'* 0.980 0.337878 0.337838 0.0000403406 *
'* 0.990 0.335612 0.335570 0.0000410721 *
'* 1.000 0.333375 0.333333 0.0000418231 *
'* *
'* End of file. *
'* Basic Version By J-P Moreau *
'* (www.jpmoreau.fr) *
'*************************************************************
defint i-n
defdbl a-h,o-z
dim c(4),x(4),y(4),z(4)
h=0.01#
a1=1.08333333333333#
a2=-2#*(a1-1#)
a3=a1-1#
f$=" ##.#### ##.###### ##.###### ##.##########"
open "stormer.lst" for output as #2
cls
print #2, "-----------------------------------------------------------------"
print #2, " Differential equation y""=f(x,y,y') by Stormer's method"
print #2, "-----------------------------------------------------------------"
'initial conditions
x(1)=0# : y(1)=1# : z(1)=-2#
xx=x(1) : gosub 1200 : yex=Fx : er=0#
print #2," X Y Y exact Error"
print #2, using f$; x(1); y(1); yex; er
for k=1 to 2
'call Runge-Kutta for first 2 steps
gosub 2000
xx=x(k+1) : gosub 1200 : yex=Fx : er=abs(yex-y(k+1))
print #2, using f$; x(k+1); y(k+1); yex; er
next k
'main Stormer loop
10 'continue
for k=2 to 4
xx=x(5-k) : yy=y(5-k) : gosub 1100
c(k)=G
next k
y(4)=2#*y(3)-y(2)+h*h*(a1*c(2)+a2*c(3)+a3*c(4))
x(4)=x(3)+h : xx=x(4) : gosub 1200 : yex=Fx : er=abs(yex-y(4))
print #2, using f$; x(4); y(4); yex; er
for k=1 to 3
x(k)=x(k+1) : y(k)=y(k+1)
next k
if x(3)<1# then goto 10 'end x value = 1
print #2,
print #2, " End of file."
Close #2
print
print " Results in file stormer.lst"
print
END
1000 'real*8 function F(x,y,z)
F=zz
return
1100 'real*8 function G(x,y,z)
G=8#*yy*yy/(1#+2#*xx)
return
1200 'exact solution real*8 function Fx(x)
Fx=1#/(1#+2#*xx)
return
2000 'SUBROUTINE RK4D2(x,y,z,h,x1,y1,z1)
xx=x(k) : yy=y(k) : zz=z(k)
gosub 1000 : c1=F
gosub 1100 : d1=G
xx=x(k)+h/2# : yy=y(k)+h/2#*c1 : zz=z(k)+h/2#*d1
gosub 1000 : c2=F
gosub 1100 : d2=G
zz=z(k)+h/2#*d2 : gosub 1000 : c3=F
yy=y(k)+h/2#*c2 : gosub 1100 : d3=G
xx=x(k)+h
zz=z(k)+h*d3 : gosub 1000 : c4=F
yy=y(k)+h*c3 : gosub 1100 : d4=G
x(k+1)=x(k)+h
y(k+1)=y(k)+h*(c1+2#*c2+2#*c3+c4)/6#
z(k+1)=z(k)+h*(d1+2#*d2+2#*d3+d4)/6#
return
'End of file stormer.bas
EXPLANATION FILE OF PROGRAM STORMER
===================================
Solve Y"=f(x,y,y') with initial conditions by Stormer's method
--------------------------------------------------------------
The differential equation of order 2 can be replaced by a system of two
equations of order 1:
Given y"=f(x,y,y') with y(a) and y'(a), by calling u = y', the problem becomes
| u'=f(x,y,u)
| y'=u
| with y(a), u(a) given
We start from a particular form of the Taylor's formula, where the remainder
is under the form of an integral:
x+h
y(x+h) = y(x) + h y'(x) + Sum (x+h-t) y"(t) dt
x
In the same way
x-h
y(x-h) = y(x) - h y'(x) + Sum (x-h-t) y"(t) dt
x
By summing
x+h x-h
y(x+h) - 2y(x) + y(x-h) = Sum ... + Sum ...
x x
For the second integral, we change the variable: u = 2x - t
So x-h x+h
Sum (x-h-t) y"(t) dt = - Sum (x+h-u) y"(u) du
x x
x+h
du=-dt ==> = Sum (x+h-t) y"(2x-t) dt
x
Finally
x+h
y(x+h) - 2y(x) + y(x-h) = Sum (x+h-t)[y"(t)+y"(2x-t)] dt
x
We now use the interpolation polynomial recalling that y"(x)=f(x,y,y'):
x
n+1
y - y + y = Sum (x - t)[P(t)+P(2x - 1)] dt
n+1 n n-1 x n+1 n
n
x
2 n+1 0 1 2
= h Sum (x - t) [O0 Div (fn) + O1 Div (fn) +O2 Div (fn)] dt
x n+1
n
m 1 |(-s) (m)|
with O = (-1) Sum (1-s) |( ) + ( )| ds
m 0 |(m ) (s)|
(see file Adambash.txt).
This leads to Stormer's formulas:
Explicit: y - 2 y + y = (h²/12)[13f - 2 f + f ]
n+1 n n-1 n n-1 n-2
Implicit: y - 2 y + y = (h²/12)[f + 10 f + f ]
n+1 n n-1 n+1 n n-1
From [BIBLI 04]
---------------------------------------
End of file Stormer.txt
No comments:
Post a Comment