MHB Backward Euler Method for 2x2 systems

[mathmari]

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) ...

 

[I like Serena]

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]

Backward Euler is:
tⁿ⁺¹=tⁿ+h
y1ⁿ⁺¹=y1ⁿ+h*f(tⁿ⁺¹,y1ⁿ⁺¹,y2ⁿ⁺¹) (1)
y2ⁿ⁺¹=y2ⁿ+h*g(tⁿ⁺¹,y1ⁿ⁺¹,y2ⁿ⁺¹) (2)
or not??

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

The backward euler for a 2x2 system is similar? this function will be used at the equation (1) only for y1ⁿ⁺¹, or for both y1ⁿ⁺¹ and y2ⁿ⁺²?

 

[I like Serena]

Backward Euler is:
tⁿ⁺¹=tⁿ+h
y1ⁿ⁺¹=y1ⁿ+h*f(tⁿ⁺¹,y1ⁿ⁺¹,y2ⁿ⁺¹) (1)
y2ⁿ⁺¹=y2ⁿ+h*g(tⁿ⁺¹,y1ⁿ⁺¹,y2ⁿ⁺¹) (2)
or not??

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

The backward euler for a 2x2 system is similar? this function will be used at the equation (1) only for y1ⁿ⁺¹, or for both y1ⁿ⁺¹ and y2ⁿ⁺²?

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]

This is the code in Matlab that I wrote :

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:

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?

 

[I like Serena]

    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:

       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]

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?

 

[I like Serena]

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]

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?

 

[I like Serena]

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]

Ok... Thank you very much! :D

 

[mathmari]

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 ?

 

[I like Serena]

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]

So the code will be:

...
        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) ??

 

[I like Serena]

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]

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:

 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()

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);

 

[I like Serena]

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]

Ok! Thank you very much! :D