Friday, September 22, 2017

The Runge-Kutta Method

'********************************************************************
'*             Differential equations of order 1                    *
'*             by Runge-Kutta method of order 4                     *
'*                                                                  *
'*                             Basic version by J-P Moreau, Paris   *
'*                                     (www.jpmoreau.fr)            *
'* ---------------------------------------------------------------- *
'* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc  *
'*             DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
'*             [BIBLI 03].                                          *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'*                                                                  *
'* Example: integrate y'=4x*(y+sqrt(y))/1+x^2 from x=0 to x=1       *
'*                                                                  *
'*        DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1          *
'*                     of type y' = f(x,y)                          *
'*                                                                  *
'*   begin value x   : ? 0                                          *
'*   end value x     : ? 1                                          *
'*   y value at x0   : ? 1                                          *
'*   number of points: ? 11                                         *
'*   finesse         : ? 10                                         *
'*                                                                  *
'*       X            Y                                             *
'* -------------------------                                        *
'*   0.000000     1.000000                                          *
'*   0.100000     1.040400                                          *
'*   0.200000     1.166400                                          *
'*   0.300000     1.392400                                          *
'*   0.400000     1.742400                                          *
'*   0.500000     2.250000                                          *
'*   0.600000     2.958400                                          *
'*   0.700000     3.920400                                          *
'*   0.800000     5.198400                                          *
'*   0.900000     6.864400                                          *
'*   1.000000     9.000000                                          *
'*                                                                  *
'********************************************************************
DefInt I-N
DefDbl A-H, O-Z

'ifi:INTEGER

  Cls
  Print
  Print "    DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1"
  Print "             of type y' = f(x,y)"
  Print
  print "    begin value x   : "; : input xi
  print "    end value x     : "; : input xf
  print "    y value at x0   : "; : input yi
  print "    number of points: "; : input m
  print "    finesse         : "; : input ifi

  Dim t(m + 1)

  'equadiff1(fp,t,xi,xf,yi,m,ifi)
  GoSub 2000

End


1000 'Example: y'=4x(y+rac(y))/(1+x²)
    fp = 4# * xx * (yy + Sqr(yy)) / (1# + xx * xx)
Return

'***************************************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER 1        *
'*                       of type y' = f(x,y)                               *
'* ----------------------------------------------------------------------- *
'*  INPUTS:                                                                *
'*    fp        Given equation to integrate (see test program)             *
'*    xi, xf    Begin, end values of variable x                            *
'*    yi        Begin Value of y at x=xi                                   *
'*    m         Number of points to calculate                              *
'*    fi        finesse (number of intermediate points)                    *
'*                                                                         *
'*  OUTPUTS:                                                               *
'*    t       : real vector storing m results for function y               *
'***************************************************************************
2000 'Subroutine  Equadiff1
  If ifi < 1 Then Return
  h = (xf - xi) / ifi / (m - 1)
  Y = yi
  t(1) = yi

  For i = 1 To m
    ni = (i - 1) * ifi - 1
    For j = 1 To ifi
      X = xi + h * (ni + j)
      xx = X: yy = Y
      GoSub 1000 'calculate current fp(x,y)
      a = h * fp
      xx = X + h / 2#: yy = Y + a / 2#
      GoSub 1000 'calculate fp(x+h/2,y+a/2)
      B = h * fp
      yy = Y + B / 2#
      GoSub 1000 'calculate fp(x+h/2,y+b/2)
      c = h * fp
      X = X + h
      xx = X: yy = Y + c
      GoSub 1000 'calculate fp(x,y+c)
      d = h * fp
      Y = Y + (a + B + B + c + c + d) / 6
    Next j
    t(i + 1) = Y
  Next i
  'Affiche(t,m,xi,xf)
  GoSub 3000
Return

3000 'Subroutine Affiche
  h = (xf - xi) / (m - 1)
  X = xi - h
  Cls
  Print
  Print "      X         Y     "
  Print "----------------------"
  For i = 1 To m
    X = X + h
    Print using; " ##.######  ##.######"; X; t(i)
  Next i
Return

'End of file teqdif1.bas



EXPLANATION FILE OF PROGRAM TEQUDIF1
====================================

  
  The Runge-Kutta Method
  ----------------------


  We present here the Runge-Kutta method of order 4 to integrate an ODE

  of order 1:       Y' = F(X, Y)

  The development of Y around x  coincidates with its Taylor development
                               n
  of order 4:

       y    = y  + h y' + (h^2/2) y" + (h^3/6) y"' + (h^4/24) y""
        n+1    n      n            n            n              n

  Other orders may also be used for the Runge-Kutta method, for instance

  the Euler's method is a Runge-Kutta of order one.

  Order 4 is often used because it is a good compromise between speed and
  
  accuracy.

  From a current point (x , y ), the next point is determined by:
                         n   n

   K  = h f(x , y )
         0        n   n

 k  = h f(x + h/2, y + K /2)
         1        n        n   0

 k  = h f(x + h/2, y + K /2)
         2        n        n   1

 k  = h f(x + h, y + K )
         3        n      n   2

  and   y   = y  + (1/6)(K + 2 K + K )
         n+1   n          0     1   3

  It can be proved that the error at each step is about h^5.

  The only drawback is the time of calculation: at each step, four function

  evaluations are necessary.

  From [BIBLI 03].

-------------------------------------------------------------
End of file rkutta.txt.

No comments:

Post a Comment