Euler's Method (First Order Runge-Kutta)

• Intro
• First Order
• Second
• Fourth
• Printable

Contents

Problem Statement

Consider the situation in which the solution, y(t), to a differential equation

is to be approximated by computer starting from some known initial condition, y(t₀)=y₀ (note that the tick mark denotes differentiation). The following text develops an intuitive technique for doing so, and then presents several examples. This technique is known as "Euler's Method" or "First Order Runge-Kutta".

Euler's Method (Intuitive)

A First Order Linear Differential Equation with No Input

Consider the following case: we wish to use a computer to approximate the solution of the differential equation

or

with the initial condition set as y(0)=3. For this case the exact solution can be determined to be (y(t)=3e^-2t, t≥0) and is shown below. Since we know the exact solution in this case we will be able to use it to check the accuracy of our approximate solution.

There are several ways to develop an approximate solution, we will do so using the Taylor Series for y(t) expanded about t=0 (in general we expand around t=t₀).

We now restrict our solution to a short time step, h, after t=0 and truncate the Taylor series after the first derivative

We call the value of the approximation y*(h), and we call the derivative y'(0)=k_1.

This is shown on the graph below for h=0.2

There are several important details to observe:

• The graph on the right is simply a detail of the shaded area of the graph on the left.
• The first order 4aylor Approximation is just a straight line starting at the initial value with a slope of -6 (i.e., k₁=y'(0)=-2y(0)=-6).
• The approximation at t=h=0.2 is just the initial value plus the slope multiplied by the time step, h; y*(h)=y*(0.2)=y(0)+k₁y'(0)=1.8 or about 10% error.
• Obviously a smaller h would result in a better approximation, and a larger h in a worse approximation.

To find the value of the approximation after the next time step, y*(2h), we simply repeat the process using our approximation, y*(h) to estimate the derivative at time h (we don't know y(h) exactly, so we can only estimate the derivative - we call this estimate k_1).

In general we move forward one step in time from t₀ to t₀+h

%% Example 1
% Solve y'(t)=-2y(t) with y0=3
y0 = 3;                  % Initial Condition
h = 0.2;% Time step
t = 0:h:2;               % t goes from 0 to 2 seconds.
yexact = 3*exp(-2*t)     % Exact solution (in general we won't know this)
ystar = zeros(size(t));  % Preallocate array (good coding practice)

ystar(1) = y0;           % Initial condition gives solution at t=0.
for i=1:(length(t)-1)
    k1 = -2*ystar(i);  % Previous approx for y gives approx for derivative
    ystar(i+1) = ystar(i) + k1*h; % Approximate solution for next value of y
end
plot(t,yexact,t,ystar);
legend('Exact','Approximate');

A First Order Linear Differential Equation with Input

Adding an input creates no particular difficulty to our approximate solution (but makes the exact solution significantly harder to find). Consider an input of cos(4t).

so

Note: the exact solution in this case is $y(t) = 2.9{e^{ - 2{\kern 1pt} t}} + 0.1\cos (4t) + 0.2\sin (4t)$.

%% Example 2
% Solve y'(t)=-2y(t)+cos(4t) with y0=3
y0 = 3;                  % Initial Condition
h = 0.2;% Time step
t = 0:h:2;               % t goes from 0 to 2 seconds.
% Exact solution, hard to find in this case (in general we won't have it)
yexact = 0.1*cos(4*t)+0.2*sin(4*t)+2.9*exp(-2*t);
ystar = zeros(size(t));  % Preallocate array (good coding practice)

ystar(1) = y0;           % Initial condition gives solution at t=0.
for i=1:(length(t)-1)
    k1 = -2*ystar(i)+cos(4*t(i));  % Previous approx for y gives approx for deriv
    ystar(i+1) = ystar(i) + k1*h; % Approximate solution for next value of y
end
plot(t,yexact,t,ystar);
legend('Exact','Approximate');

A First Order Non-Linear Differential Equation

Non-linear differential equations can be very difficulty to solve analytically, but pose no particular problems for our approximate method. Consider the differential equation given by

the solution is (described here)

To obtain the approximate solution, we write an expression for y'(t) and use it to find k₁.

%% Example 3
% Solve y'(t)=y(t)(1-2t) with y0=1
y0 = 1;                  % Initial Condition
h=0.2;                   % Time step
t = 0:h:2;               % t goes from 0 to 2 seconds.
% Exact solution, hard to find in this case (in general we won't have it)
yexact = exp(t-t.^2);
ystar = zeros(size(t));  % Preallocate array (good coding practice)

ystar(1) = y0;           % Initial condition gives solution at t=0.
for i=1:(length(t)-1)
  k1 = ystar(i)*(1-2*t(i));   % Approx for y gives approx for deriv
  ystar(i+1) = ystar(i) + k1*h; % Approximate solution at next value of y
end
plot(t,yexact,t,ystar);
legend('Exact','Approximate');

A Higher Order Linear Differential Equation

Though the techniques introduced here are only applicable to first order differential equations, the technique can be use on higher order differential equations if we reframe the problem as a first order matrix differential equation. Consider the 3^rd order equation (with initial conditions)

If we introduce new variables, q₁, q₂, and q₃, we can write three first order differential equations

We can further simplify by writing this as a single first order matrix differential equation.

A more in depth discussion of the procedure of going from n^th order differential equation to first order nxn matrix differential equation is here.

We now proceed as before except the variable being integrated and the quantity k₁ are both vectors.

Note: the exact solution to this problem is $y(t) = {1 \over 4}+{{\rm{e}}^{ - t}}{\mkern 1mu} \left( {\cos (t) - {5 \over 2}\sin (t)} \right) - {5 \over 4}{{\rm{e}}^{ - 2{\kern 1pt} t}}$, for t≥0.

%% Example 4
% Solve y'''(t)+4y''(t)+6y'(t)+4y(t)=gamma(t),  y''(0)=0, y'(0)=-1, y(0)=0
q0 = [0; -1; 0];         % Initial Condition (vector)
h = 0.2;% Time step
t = 0:h:5;               % t goes from 0 to 2 seconds.
A = [0 1 0; 0 0 1; -4 -6 -4];  % A Matrix
B = [0; 0; 1];                 % B Matrix

yexact = 1/4 + exp(-t).*(cos(t) - 5*sin(t)/2) - 5*exp(-2*t)/4;  % Exact soln
qstar = zeros(3,length(t));  % Preallocate array (good coding practice)

qstar(:,1) = q0;             % Initial condition gives solution at t=0.
for i=1:(length(t)-1)
  k1 = A*qstar(:,i)+B*1;            % Approx for y gives approx for deriv
  qstar(:,i+1) = qstar(:,i) + k1*h; % Approximate solution at next value of q
end
plot(t,yexact,t,qstar(1,:));        % ystar = first row of qstar
legend('Exact','Approximate');
hold off

Euler's Method (The Math)

The math for this method, the first order Runge-Kutta (or Euler's Method) is fairly simple to understand, and has been discussed before. If we write the differential equation as

and write the approximation to the derivative as

We expand y(t) around t₀ assuming a time step h

and drop all terms after the linear term. Because all of the dropped terms are multiplied by h² or greater, we say that the algorithm is accurate to order h² locally, or O(h²) (if h is small the other terms that are multiplied by h³, h⁴... which will be even smaller, and can be dropped as well).

This gives us our approximate solution at the next time step.

$$y^*({t_0} + h) = y^*({t_0}) + {k_1}h$$

Sinc the number of steps over the whole interval is proportional to 1/h (or O(h_^-1)) we might expect the overall accuracy to be the O(h²)·O(h_^-1)=O(h). A rigourous analysis proves that this is true. Such an analysis can be found in references about numerical methods such as the book Applied Numerical Methods, by Carnahan, Luther and Wilkes. Because the error is O(h), where h is raised to the power 1, this is called a first order algorithm.

Moving on

With a little more work we can develop an algorithm that is accurate to higher order than O(h). The next page describes just such a method.

References