## Differential Equations

### © 2003, Ingemar Bjerle

#### ingemar.bjerle@telia.com

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:

Application Window

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`