---

In formulating DE it is essential to understand clearly the meaning of a derivative dx/dt. The symbol dx/dt merely states "the rate of change of x with respect to t". If x is related to t, then dx/dt is the slope of the curve at any point t. A often used shorthand for dx/dt is x´.

Example: x´=v. This DE is used very often by car drivers. Assume you are out driving at a velocity of v=50 mph. What distance have you been driving in 1 hour? Just by multiplying 50 with one you get the distance x=50 miles. You have perhaps without knowing it used the integrated form of x´=v, that is x=v*t. By integrating x´=v you get the formula. By rewriting x´=v to dx=v*dt and by integrating both sides the result is:

Now if the velocity varies, v is a function of time and now you must let the car computer take over the calculation. With constant velocity you can set dt=1h but now if the velocity is read every second by the computer dt=1s and v is a function of t and v is in this case written as v(t). The integrated form will now look like:

If v(t) increases linear from 0 to 50 mph in one hour, v(t)=t*50/3600

In a for-next loop the calculation is done

    x=0: dt=1:t=0
    for i =1 to 3600
    vt=t*50/3600
    x=x+vt/3600*dt
    t=t+dt
    next i
    print x	;"    ";t'	x=25 mph

x´= v is a first order DE and always by handling DE on a computer they must be of first order. What to do if the DE is of higher order x´´ or higher? A second order DE can always be reduced to two first order equations and a third order to three equations. The first step in preparing a DE for numerical solving is to reduce the order.

Preparing a DE for numerical solution:

Derivation of a function y

    a=2, b=3
    y  = x^2+a*x+b
    y´ = 2*x+a
    y´´= 2  ( original DE)

Integration of y´´:

    y´=2*a
    y=x^2+ax+b

The DE y´´=2 will now be converted to two first order equations.

The first substitution that always should be made is y´= u. This is the first equation that will be used.

The next step is to get the derivative of the substitution that is y´´=u´.

The second equation can be found by inserting the derivative of the substitution (y´´=u´) in the original DE

giving u´ =2

We now have two first order DE:

    y´ =u
    u´ =2

with initial boundary conditions:

    x = 0
    y = b   (according to: y  = x^2+a*x+b with x=0 inserted)
    y´= a   (according to: y´ = 2*x+a with x=0 inserted)

Integration methods:

The most simple integration method is the EULER. In this method from the slope at the beginning of the interval dx, dy and du are calculated from dy=u*dx and du=2*dx. The method requires a small interval size to give accurate results.

A very popular and accurate method is the RUNGE KUTTA. In this method 4 estimates of dy and du are made: one at the beginning two in the middle and one at the end of the interval. These four values are weighed together.

Both methods are given in the program below.

Finally a LB3 program with a parser is enclosed. With that program it is possible to solve any initial boundary value problem with up to five first order DE. The program is using Runge Kutta. Be aware that not all DE give a stable solution. However in most cases it is possible to find a x-region with a stable solution.

   'PROGRAM WITH EULER AND RUNGE KUTTA SUBROUTINES
    gosub [INI]
    PRINT" X       RUNGE     EULER"
    for i=1 to n
     gosub [RUNGE]
     print X;"    ";Y;"    ";
     gosub [EULER]
     print Y1
    next i
    PRINT" X       RUNGE     EULER"
    PRINT "STEPS= ";n
    input r$
    end

[INI]
    a=2:b=3
    Xin=0:  Xend=10
    Yin=b:  Uin=a
    Y=Yin:  U=Uin:  X=Xin
    Y1=Y:   U1=U:  X1=X
    n=100
    DX=(Xend-Xin)/n
    RETURN

[EQ]       '  SUBROUTINE EQUATIONS
   REM  DIFF EQ:  D2Y/DX2=2
   REM  WRITTEN AS TWO FIRST ORDER:  DY/DX=U:     DU/DX=2
   REM  BOUNDARY CONDITIONS: Xin0: Xend=10: Yin=b: Uin=a
        'A=DY/DX; B=DU/DX
       '*********************
        A=U
        B=2
       '*********************
       RETURN




[RUNGE]      ' SUBROUTINE RUNGE
  rem  VARIABLES IN SUBROUTINE: INDEPENTENT: X  DEPENDENT:  Y, U
  GOSUB [EQ]
      		          K1=A*DX:  L1=B*DX
       X=X+DX/2:     Y=Y+K1/2:  U=U+L1/2
  GOSUB [EQ] :        K2=A*DX:  L2=B*DX
                     Y=Y-K1/2:  U=U-L1/2
                     Y=Y+K2/2:  U=U+L2/2
  GOSUB [EQ] :        K3=A*DX:  L3=B*DX
                     Y=Y-K2/2:  U=U-L2/2
        X=X+DX/2:      Y=Y+K3:  U=U+L3
  GOSUB [EQ]:         K4=A*DX:  L4=B*DX
                       Y=Y-K3:  U=U-L3
      Y=Y+K1/6+K2/3+K3/3+K4/6:  U=U+L1/6+L2/3+L3/3+L4/6
       return

[EQ1]       '  SUBROUTINE EQUATIONS
  REM  DIFF EQ:  D2Y/DX2=2
  REM  WRITTEN AS TWO FIRST ORDER:  DY/DX=U1:     DU1/DX=2
  REM  BOUNDARY CONDITIONS: Xin=0: Xend=10: Yin=b: Uin=a
        'A=DY/DX:   B=DU1/DX
       '*********************
        A=U1
        B=2
       '*********************
       RETURN

[EULER]      ' SUBROUTINE EULER
   rem  VARIABLES IN SUBROUTINE: INDEPENTENT: X1  DEPENDENT:  Y1, U1
        GOSUB [EQ1]
        K1=A*DX:  L1=B*DX
        Y1=Y1+K1: U1=U1+L1
        X1=X1+DX
       Return

Enclosed program: diff_eq-parser.bas

datafile: diff_eq2.dat - case 8

---