From Euler forward method in to Euler backward.

Question

I have a system of linear differential equations

$$y'=Ay+b,$$

where

A = [0 0 0 1;0 0 1 0;0 0 1 0;0 0 0 -1];
b = [0; 0; 1; 0];
y0 = [1; 0; 0; -2];
tinit = 0;
h = 0.1;
tFinal = 2;

Solving it with Euler forward, worked nicely by using $$y_{k+1}=y_k+hf(t_{k},y_k),$$

the functioncode is:

function yOut = eulerforward(F,t0,h,tFinal,y0)
    Y = y0;
    yOut = Y;
    for t = t0 : h : tFinal
            s = F(t,Y);
            Y = Y + h*s;
        yOut = [yOut Y];
    end
end

and finally

F = @(t,Y) A*Y + b;
K1 = eulerforward(F,tinit,h,tFinal,y0);

Now, Euler backward is formulated by

$$y_{k+1}=y_k+hf(t_{k+1},y_{k+1}).$$

How can I modify my Euler forward code so that it uses the Euler backward method?

Is $s$ the only parameter that will change?

Answer 1

using backward Euler you need to find the zero of $F(\mathbf{y})=\mathbf{y} - \mathbf{y_n} - kf(\mathbf{y})$, where $y_n$ is the solution at time $n$. Of course this zero is the next value of the solution. Usually you need to solve a non-linear system, but if the problem is linear, then you have just to solve a linear system at every iteration (for example using $\verb|\|$ in Octave).

In this case $f(y)=Ay+b$, so the numerical scheme is:

$y_{n+1}=y_n + k(Ay_{n+1} + b)$, so we have $(I-kA)\cdot y_{n+1}=y_n + kb$, with the initial value $y_0$ given.

If you want to use a function called "backward Euler" you will need to use Newton's method.

Here's a runable MatLab code for your specific problem.

y0=[1; 0; 0; -2];
tinit = 0;
h = 0.1;  
tFinal = 2;
ts=tFinal/h; %timesteps


%matrix
A = [0 0 0 1;0 0 1 0;0 0 1 0;0 0 0 -1];
b = [0; 0; 1; 0];

y=NaN(4,ts+1);
b = [0; 0; 1; 0];
y(:,1)=y0;
%loop
for n=1:ts
    y(:,n+1)=(eye(4)-h*A)\(y(:,n)+h*b);   
end
t=linspace(tinit,tFinal,ts+1);
plot(t,y(1,:),'b-o',t,y(2,:),'r-o',t,y(3,:),'g-o',t,y(4,:),'k-o')
legend('first','second','third','fourth equation')