Numerical methods for 1D ODEs

Consider the IVP:

, with the initial condition

We want to compute the numerical approximation for

, that is, to compute

for

using different numerical methods.

Euler's method

where h is the time step.

Notice, from the numerical results below, that the error is proportional to

We will approximate the value using N=10, 100, 1000, 10000 subdivisions

clear all;
f=@(t,x) x;
sol=exp(1); %we know the analytical solution
tini = 0;
xini = 1; %x(0)=1
tfin = 1;
N = [10,100,1000,10000];
t(1) = tini;
x(1) = xini;
for k = 1:4
    h = (tfin-tini)/N(k);
    npoints = N(k)+1;
    for i=2:npoints
        x(i) = x(i-1)+h*f(t(i-1),x(i-1));
        t(i)= t(i-1)+h;
    end
    fprintf('for N= %d, value x(1)= %14.12f error=%20.12e\n',N(k),x(end),abs(x(end)-exp(1)));
end
for N= 10, value x(1)= 2.593742460100 error=  1.245393683590e-01
for N= 100, value x(1)= 2.704813829422 error=  1.346799903752e-02
for N= 1000, value x(1)= 2.716923932236 error=  1.357896223150e-03
for N= 10000, value x(1)= 2.718145926825 error=  1.359016338189e-04

Mid point Method

Similar to the Euler method we can use the midpoint method:

Notice, from the numerical results below, that the error is proportional to

clear x t;
t(1) = tini;
x(1) = xini;
for k = 1:4
    h = (tfin-tini)/N(k);
    npoints = N(k)+1;
    for i=2:npoints
        x(i) = x(i-1)+h*f(t(i-1)+0.5*h,x(i-1)+0.5*h*f(t(i-1),x(i-1)));
        t(i) = t(i-1)+h;
    end
    fprintf('for N= %d, value x(1)= %14.12f error=%20.12e\n',N(k),x(end),abs(x(end)-exp(1)));
end
for N= 10, value x(1)= 2.714080846608 error=  4.200981850822e-03
for N= 100, value x(1)= 2.718236862560 error=  4.496589908864e-05
for N= 1000, value x(1)= 2.718281375752 error=  4.527072827720e-07
for N= 10000, value x(1)= 2.718281823929 error=  4.530157582394e-09

Implicit Euler Method

The implicit or backward Euler's method consists in taken the slope at the next point (already unknown).

This is an implicit equation on

. Because of that, unless any other information or simplification about the ODE is know, you have to solve a non-linear numerical equation at each time step.

An usual approach is to approximate using taylor expansion

does'nt depends on t,

. Therefore only the

is computed and then

and using

, then we obtain

that must be solved for

For the exemple we are dealing

and therefore

and

Thus, the implicit Euler method for this exemple can be written as

Exercise: Build the Matlab code for this example and compute the results for N=10, 100, 1000, 10000 subdivisions.

Runge-Kutta Methods

Runge-Kutta 2

with

where

Exercise: Build the Matlab code for this example and compute the results for N=10, 100, 1000, 10000 subdivisions.

Notice, from the numerical results below, that the error is proportional to

Runge-Kutta 4

with

where

Exercise: Build the Matlab code for this example and compute the results for N=10, 100, 1000, 10000 subdivisions.

Notice, from the numerical results below, that the error is proportional to

Home work:

Exercise1: Build the Euler Method as a Matlab function defined by

function [ t, X ] = eulerMethod ( f, timeRange, xini, N)

input:

f: function address for the differential equations

timeRange=[ tini, tfin ]: time interval 1x2 vector

xini: initial values nx1 vector (n the dimension of the problem)

N: number of subdivisions

output:

t: vector with all the time steps. t(1) = tini; t(end) =tfin;

X: trajectory with all the computed points, (N+1)xn Matrix

Exercise 2:

Do the same with the Mid Point Method as a Matlab function defined by

function [ t, X ] = midPointMethod ( f, timeRange, xini, N)

Exercise 3:

Do the same with the RK2 Method as a Matlab function defined by

function [ t, X ] = RK2Method ( f, timeRange, xini, N)

Exercise 4:

Do the same with the RK4 Method as a Matlab function defined by

function [ t, X ] = RK4Method ( f, timeRange, xini, N)

Exercise 5:

Apply all these methods to the solve the ODE

for

and

and compare with the exact solution

where