2.4.4 Euler's Method for Systems of Differential Equations

In the next example, we will illustrate Euler's method for first and second order ODEs. We first recall the basic idea for first order equations. Given an initial value problem of the form

$\displaystyle y'=f(x,y), \quad y(a)=c,$    

we want to find the approximate value of the solution at x=b for any given b with b>a.

Recall from the definition of the derivative that

$\displaystyle y'(x)\approx \frac{y(x+h)-y(x)}{h},$    

where h>0 is given and small. This and the DE together give $ f(x,y(x))\approx
\frac{y(x+h)-y(x)}{h}$ . Now solve for y(x+h):

$\displaystyle y(x+h)\approx y(x)+h f(x,y(x)).$    

If we call h f(x, y(x)) the ``correction term'' (for lack of anything better), call y(x) the ``old value of y'', and call y(x+h) the ``new value of y'', then this approximation can be re-expressed as

$\displaystyle y_{\text{new}} \approx y_{\text{old}}+h\cdot f(x,y_{\text{old}}).$    

If we break the interval from a to b into n steps, so that h=(b-a)/n, then we can record the information for this method in a table.

x y hf(x,y)
a c hf(a,c)
a+h c+hf(a,c) ...
a+2h ...  
...    
b=a+nh ??? ...

The goal is to fill out all the blanks of the table, one row at a time, until we reach the ??? entry, which is the Euler's method approximation for y(b).

The idea for systems of ODEs is similar.

Example: Numerically approximate z(t) at t=1 using 4 steps of Euler's method, where $ z''+tz'+z=0$ , $ z(0)=1$ , $ z'(0)=0$ .

One must reduce the 2nd order ODE down to a system of two first order DEs (using $ x=z$ , $ y=z'$ ) and apply Euler's method:

sage: t,x,y = PolynomialRing(RealField(10),3,"txy").gens()
sage: f = y; g = -x - y * t
sage: eulers_method_2x2(f,g, 0, 1, 0, 1/4, 1)
  t                x            h*f(t,x,y)                y       h*g(t,x,y)
  0                1               0.00000                0         -0.25000
1/4           1.0000             -0.062500         -0.25000         -0.23438
1/2          0.93750              -0.11719         -0.46875         -0.17578
3/4          0.82031              -0.15381         -0.61523        -0.089722
  1          0.66602              -0.16650         -0.66602          0.00000
Therefore, $ z(1)\approx 0.75$ .

We can also plot the points $ (x,y)$ to get an approximate picture of the curve. The function eulers_method_2x2_plot will do this; in order to use it, we need to define functions f and g which takes one argument with three coordinates: (t, x, y).

sage: f = lambda z: z[2]        # f(t,x,y) = y
sage: g = lambda z: -sin(z[1])  # g(t,x,y) = -sin(x)
sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)

At this point, P is storing two plots: P[0], the plot of x vs. t, and P[1], the plot of y vs. t. We can plot both of these as follows:

sage: show(P[0] + P[1])
(For more on plotting, see Section 2.5.)

See About this document... for information on suggesting changes.