Welcome to our community

Be a part of something great, join today!

Backward Euler Method for 2x2 systems

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
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) ... :confused:
 

Klaas van Aarsen

MHB Seeker
Staff member
Mar 5, 2012
8,780
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) ... :confused:
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
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047
Ok... Thank you very much!!!!! :D
 

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047
Ok! Thank you very much!!!!!!!! :eek: :D