Lab 9: Euler's method for systems

All the Second Project options rely on using a numerical method to approximate solutions of a system of 1st order differential equations.

In Lab 3, you used Euler's method on a single 1st order DE. Here we extend it to work on systems, i.e. a collection of more than one DE to be solved simultaneously.

Here is Euler's method applied to the Emily and Jacob system. You will notice that the Euler step looks exactly as in Lab 3. The difference is that X and F are now vectors rather than scalars (numbers).

1. Euler's method for a system

Copy and run the following code:

from resources306 import *

def euler_step(F,t,X,h):
         newX = X + h*F(t,X)  # this vector arithmetic is done elementwise
         newt = t + h
         return newt,newX

def FEJ(t,X):
         # This function computes and returns the rates of change of the dependent variables
         # X' for the Emily and Jacob system
         # the independent variable t may or may not be needed (here not)
         E,J = X  # unpack elements of the dependent variable vector X
         Eprime = J
         Jprime = -E
         return np.array([Eprime,Jprime])

t0     = 0
tfinal = 8
nsteps = 20
h = (tfinal-t0)/nsteps
d = 2            # the dimension of system

data = np.zeros((nsteps+1,1+d))  # we will store t and X in the rows of this array

t = t0
X = 0.,1.      # Here is where the initial values of the dependent variables are specified
data[0] = t,*X

for i in range(nsteps):
         t,X = euler_step(FEJ,t,X,h)
         data[i+1] = t,*X

print(data)

Observe that the t-values are in the first column of the "data" array, and the E and J values are in the 2nd and 3rd columns respectively.

2. Plotting the results

Make some plots of the results:

E and J vs. t

t,E,J = data.T
plt.plot(t,E,label='E')
plt.plot(t,J,label='J')
plt.xlabel('t')
plt.legend();

The phase plane picture, E vs J

plt.subplot(111,aspect=1) # force the scales on horizontal and vertical axes to be the same
plt.plot(E,J)
plt.xlabel('E'); plt.ylabel('J',rotation=0);

3. Choosing an appropriate step size

The step size used above is much too large to obtain an accurate approximation of the exact solution. You can verify this by running the code again with a larger number of steps (hence smaller step size). Try 10 times as many steps:

t0     = 0
tfinal = 8
nsteps = 200
h = (tfinal-t0)/nsteps
d = 2 # dimension of system

data = np.zeros((nsteps+1,1+d))  # we will store t and X in the rows of this array

t = t0
X = 0.,1.       #  initial values of the dependent variables
data[0] = t,*X

for i in range(nsteps):
         t,X = euler_step(FEJ,t,X,h)
         data[i+1] = t,*X

t,E,J = data.T
plt.figure()
plt.plot(t,E,label='E')
plt.plot(t,J,label='J')
plt.xlabel('t')
plt.legend();

plt.figure()
plt.subplot(111,aspect=1) # force the scales on horizontal and vertical axes to be the same
plt.plot(E,J)
plt.xlabel('E'); plt.ylabel('J',rotation=0);

You will see the approximate solution curve in the phase plane is almost a closed loop. In this system, we happen to know the exact solution: it is E(t) = sin(t), J(t) = cos(t), which is exactly a circular cycle in the E,J plane.

To get an accurate numerical approximation, we will need to

4. Enough steps for good accuracy

See how many steps you need to take with Euler's method to get a loop that looks pretty much closed in the phase plane plot. You will find it is quite a large number!

5. A better method than Euler's

Look up the Runge-Kutta method in the textbook, p73, Exercise 1.7.6, and implement it.

You will write a function to replace "euler_step" that has this form:

def rk_step(F,t,X,h):

        ...  # implement the RK formulas here


        newX = X + h* ...
        newt = t + h
        return newt,newX

Note that the textbook's xi is our t, and their yi is our X. Otherwise you can more or less just transcribe the formulas on p73.

6. Test the Runge-Kutta method

See how much more accurate it is than Euler's method on the E,J system.

I recommend you use the RK method you have just written for your Second Project.

7. Application to Van der Pol system

If you have time, it might be interesting to apply your RK method to the Van der Pol oscillator in which we saw a "limit cycle" on Day 20:

x′ = y

y′ =  − x + (1 − x2)y

with an initial condition near, but not exactly at, (0,0). Compare with the pictures in the textbook on p288, where tfinal = 30.