'************************************************************* '* 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