Backward Euler Method for 2x2 systems

In summary, the conversation discusses backward Euler method for 2x2 systems in Matlab using fixed point iteration. The method involves calculating y1n+1 and y2n+1 at equation (1) and (2) respectively, but y2n+1 is already used in equation (1). The process can be generalized for other problems and the code in Matlab involves a function called stage that calculates the next iteration. It is important to not use y1(n+1) and y2(n+1) in the stage function, as they should be calculated by it.
  • #1
mathmari
Gold Member
MHB
5,049
7
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:
 
Mathematics news on Phys.org
  • #2
mathmari said:
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.
 
  • #3
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?
 
  • #4
mathmari said:
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.$$
 
  • #5
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?
 
  • #6
mathmari said:
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);
 
  • #7
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?
 
  • #8
mathmari said:
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).
 
  • #9
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?
 
  • #10
mathmari said:
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.
 
  • #11
Ok... Thank you very much! :D
 
  • #12
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 ?
 
  • #13
mathmari said:
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.
 
  • #14
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 by a moderator:
  • #15
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.
 
  • #16
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);
 
  • #17
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.
 
  • #18
Ok! Thank you very much! :eek: :D
 

FAQ: Backward Euler Method for 2x2 systems

1. What is the Backward Euler Method for 2x2 systems?

The Backward Euler Method is a numerical method used to approximate the solutions of a system of two differential equations. It is an implicit method, meaning that it uses the value of the dependent variable at the next time step to approximate the value at the current time step.

2. How does the Backward Euler Method work?

The Backward Euler Method uses the backward difference formula to approximate the derivatives of the dependent variables. This approximation is then used in the original differential equations to calculate the values of the dependent variables at the next time step. This process is repeated until the desired number of time steps is reached.

3. When is the Backward Euler Method useful?

The Backward Euler Method is useful when solving stiff systems of differential equations, where the step size needs to be very small to ensure accuracy. It is also useful when the system is nonlinear, as it can handle larger changes in the dependent variables compared to other methods.

4. What are the advantages of using the Backward Euler Method?

One advantage of the Backward Euler Method is that it is unconditionally stable, meaning that it can handle large step sizes without causing numerical instability. It is also relatively easy to implement and can handle a wide range of systems of differential equations.

5. What are the limitations of the Backward Euler Method?

One limitation of the Backward Euler Method is that it is a first-order method, meaning that the error in the approximations can be relatively large compared to higher-order methods. It also requires the solution of nonlinear equations at each time step, which can be computationally expensive.

Similar threads

Back
Top