# Backward Euler Method for 2x2 systems

#### mathmari

##### Well-known member
MHB Site Helper
Hi!!!
I want to write a code in Matlab for the Backward Euler Method for 2x2 systems, using the fixed point iteration to find the yn+1.
y1n+1=y1n+h*f(tn+1,y1n+1,y2n+1) (1)
y2n+1=y2n+h*g(tn+1,y1n+1,y2n+1) (2)
Could you tell how I use the fixed point iteration??
At (1) the fixed point iteration will calculate y1n+1, y2n+1 will be calculated at (2) but it is already used in the equation (1) ...

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Hi!!!
I want to write a code in Matlab for the Backward Euler Method for 2x2 systems, using the fixed point iteration to find the yn+1.
y1n+1=y1n+h*f(tn+1,y1n+1,y2n+1) (1)
y2n+1=y2n+h*g(tn+1,y1n+1,y2n+1) (2)
Could you tell how I use the fixed point iteration??
At (1) the fixed point iteration will calculate y1n+1, y2n+1 will be calculated at (2) but it is already used in the equation (1) ...
Hi!

I think that should be:
\left\{\begin{aligned} t_{n+1}&=t_n+h \\ y_{1,n+1}&=y_{1n}+h\ f(t_n,y_{1n},y_{2n}) \\ y_{2,n+1}&=y_{2n}+h\ g(t_n,y_{1n},y_{2n}) \\ \end{aligned}\right.
where $h$ is the time step,
$f(t_n,y_1,y_2)$ is the derivative with respect to time of $y_1$,
and $g(t_n,y_1,y_2)$ is the derivative with respect to time of $y_2$.

There are more advanced methods, like the method of Heun, and notably the method of Runge Kutta.
The latter is usually preferred.
Either way, Euler's method is based on the last iteration and not on the iteration that is coming up.

#### mathmari

##### Well-known member
MHB Site Helper
Backward Euler is:
tn+1=tn+h
y1n+1=y1n+h*f(tn+1,y1n+1,y2n+1) (1)
y2n+1=y2n+h*g(tn+1,y1n+1,y2n+1) (2)
or not??

The backward euler for y'=F(t,y) is
yn+1=yn+h*f(tn+1,stage(yn,h,tn))
where stage is a function that appreciates the fixed point to calculate the yn+1.

The backward euler for a 2x2 system is similar??? this function will be used at the equation (1) only for y1n+1, or for both y1n+1 and y2n+2???

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Backward Euler is:
tn+1=tn+h
y1n+1=y1n+h*f(tn+1,y1n+1,y2n+1) (1)
y2n+1=y2n+h*g(tn+1,y1n+1,y2n+1) (2)
or not??

The backward euler for y'=F(t,y) is
yn+1=yn+h*f(tn+1,stage(yn,h,tn))
where stage is a function that appreciates the fixed point to calculate the yn+1.

The backward euler for a 2x2 system is similar??? this function will be used at the equation (1) only for y1n+1, or for both y1n+1 and y2n+2???
Ah, my mistake. I only just now read the article about backward Euler.

Anyway, the article says, that you can (sometimes) use the process:
$$y_{k+1}^{[0]} = y_k, \quad y_{k+1}^{[i+1]} = y_k + h f(t_{k+1}, y_{k+1}^{})$$
to find the next iteration $y_{k+1}$ by running i from 0 to a high enough value.

Generalizing to your problem, that becomes:
\left\{\begin{aligned} y_{1,k+1}^{[0]} &= y_{1k}, &\quad y_{1,k+1}^{[i+1]} &= y_{1k} + h f(t_{k+1}, y_{1,k+1}^{},y_{2,k+1}^{}) \\ y_{2,k+1}^{[0]} &= y_{2k}, &\quad y_{2,k+1}^{[i+1]} &= y_{2k} + h f(t_{k+1}, y_{1,k+1}^{},y_{2,k+1}^{}) \end{aligned}\right.

#### mathmari

##### Well-known member
MHB Site Helper
This is the code in Matlab that I wrote :

Code:
function [t,y1,y2]=BackwardEulerSystem(y1_0,y2_0,t0,T,N)
h=(T-t0)/N;
t=zeros(1,N+1);
y1=zeros(1,N+1);
y2=zeros(1,N+1);
y1(1)=y1_0;
y2(1)=y2_0;
t(1)=t0;
for n=1:N
[Y1,Y2]=stage(y1(n),y2(n),y1(n+1),y2(n+1),h,t(n));
t(n+1)=t(n)+h;
y1(n+1)=y1(n)+h*S(t(n+1),Y1,y2(n+1));
y2(n+1)=y2(n)+h*G(t(n+1),y1(n+1),Y2);
end
t=t';
y1=y1';
y2=y2';
where the function stage is:
Code:
function [Y1,Y2]=stage(y1n,y2n,y1N,y2N,h,tn)
M=5;
TOL=1e-5;
x1(1)=y1n+h*S(tn,y1N,y2N);
x2(1)=y2n+h*G(tn,y1N,y2N);
for m=1:M
x1(m+1)=y1n+h*S(tn+h,x1(m),x2(m));
x2(m+1)=y2n+h*G(tn+h,x1(m),x2(m));
if (abs(x1(m+1)-x1(m))<TOL || abs(x2(m+1)-x2(m))<TOL)
Y1=x1(m+1);
Y2=x2(m+1);
return
end
x1_m=x1(m+1);
x2_m=x2(m+1);
end
Y1=x1(M);
Y2=x2(M);
for example, S=-y2, G=y1...

At the function BackwardEulerSystem, when I call the function stage(y1(n),y2(n),y1(n+1),y2(n+1),h,t(n)),
is the y1(n+1) and y2(n+1) right or do I have to replace them with something else???

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Code:
    for n=1:N
[Y1,Y2]=stage(y1(n),y2(n),y1(n+1),y2(n+1),h,t(n));
t(n+1)=t(n)+h;
y1(n+1)=y1(n)+h*S(t(n+1),Y1,y2(n+1));
y2(n+1)=y2(n)+h*G(t(n+1),y1(n+1),Y2);
end
Zooming in on this part for now.

You shouldn't use y1(n+1) when it hasn't been assigned a value yet.
So your stage() function shouldn't have y1(n+1),y2(n+1) as input parameters.
And your next iteration should be like:
Code:
       y1(n+1)=y1(n)+h*S(t(n+1),Y1,Y2);
y2(n+1)=y2(n)+h*G(t(n+1),Y1,Y2);

#### mathmari

##### Well-known member
MHB Site Helper
So, do you mean that I have not to use y1(n+1) and y2(n+1) at all at the function stage(), or could I replace them with some other parameters???

#### Klaas van Aarsen

##### MHB Seeker
Staff member
So, do you mean that I have not to use y1(n+1) and y2(n+1) at all at the function stage(), or could I replace them with some other parameters???
You should not use them at all.
They are supposed to be calculated by the stage() function.
Its output [Y1,Y2] is the actual y1(n+1) and y2(n+1).

#### mathmari

##### Well-known member
MHB Site Helper
So, I have also change the function stage(), since I used as arguments the y1N,y2N (that corresponds to y1(n+1), y2(n+1)), right??? But how can I do that???

#### Klaas van Aarsen

##### MHB Seeker
Staff member
So, I have also change the function stage(), since I used as arguments the y1N,y2N (that corresponds to y1(n+1), y2(n+1)), right??? But how can I do that???
The stage() function should implement the following algorithm:
\left\{\begin{aligned} y_{1,k+1}^{[0]} &= y_{1k}, &\quad y_{1,k+1}^{[i+1]} &= y_{1k} + h f(t_{k+1},\ y_{1,k+1}^{},\ y_{2,k+1}^{}) \\ y_{2,k+1}^{[0]} &= y_{2k}, &\quad y_{2,k+1}^{[i+1]} &= y_{2k} + h f(t_{k+1},\ y_{1,k+1}^{},\ y_{2,k+1}^{}) \end{aligned}\right.

The input is $y_{1k}$ and $y_{2k}$, which are in your case y1(n) and y2(n).
Then, after a number of iterations you'll get $y_{1,k+1}^{}$ and $y_{2,k+1}^{}$, which are actually y1(n+1) and y2(n+1).
Those 2 are the output of the algorithm and should be returned as Y1 and Y2.

#### mathmari

##### Well-known member
MHB Site Helper
Ok... Thank you very much!!!!!

#### mathmari

##### Well-known member
MHB Site Helper
Now I'm studying again the code and I'm facing difficulties... Could you explain me why at the function stage(), the "if" loop has to stop either if |(x1(m+1)-x1(m)|<TOL or |x2(m+1)-x2(m)|<TOL and not if both of them are <TOL ????

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Now I'm studying again the code and I'm facing difficulties... Could you explain me why at the function stage(), the "if" loop has to stop either if |(x1(m+1)-x1(m)|<TOL or |x2(m+1)-x2(m)|<TOL and not if both of them are <TOL ????
Both of them need to be below TOL.
So that's a mistake in the code.

#### mathmari

##### Well-known member
MHB Site Helper
So the code will be:
Code:
......
if (abs(x1(m+1)-x1(m))<TOL && abs(x2(m+1)-x2(m))<TOL)
Y1=x1(m+1);
Y2=x2(m+1);
return
end
......
???
or do have I to use two "if" loops or something else, because at the above "if", if abs(x1(m+1)-x1(m))<TOL but abs(x2(m+1)-x2(m))<TOL is not true, Y1 does not get the value x1(m+1) ??

Last edited:

#### Klaas van Aarsen

##### MHB Seeker
Staff member
It is correct as you have it now.
If one of the two is not within TOL, the surrounding loop must continue until both are.
When they finally are both within TOL, Y1 and Y2 will both get their final value.

#### mathmari

##### Well-known member
MHB Site Helper
Great..
I compared my results with the results using ode45 and they are not the same...
I do not know where my mistake is...
Either I did something wrong at the following "for" loop:
Code:
 for n=1:N
[Y1,Y2]=stage(y1(n),y2(n),h,t(n));
t(n+1)=t(n)+h;
y1(n+1)=y1(n)+h*S(t(n+1),Y1,Y2);
y2(n+1)=y2(n)+h*G(t(n+1),Y1,Y2);
end
or at the following loop of the function stage()
Code:
for m=1:M
x1(m+1)=y1n+h*S(tn+h,x1(m),x2(m));
x2(m+1)=y2n+h*G(tn+h,x1(m),x2(m));
if (abs(x1(m+1)-x1(m))<TOL && abs(x2(m+1)-x2(m))<TOL)
Y1=x1(m+1);
Y2=x2(m+1);
return
end
x1(m)=x1(m+1);
x2(m)=x2(m+1);
end
Y1=x1(M+1);
Y2=x2(M+1);

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Well, both fragments look okay to me...

So you may need to "debug" your code.

You can do so by printing intermittent results with both methods and see if they are similar.
And what usually works pretty well, is make a graph of the points that are found.
Often you can see what's going wrong then.

#### mathmari

##### Well-known member
MHB Site Helper
Ok! Thank you very much!!!!!!!!