Thursday, October 12, 2017

Stormer's method

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