# Forward Euler Method for ODE system

• stvoffutt
The above is true but misses the point. Forward Euler is known to be unstable for this problem over the interval [0,5] (hence my use of "overestimate" and "underestimate"). The question is what you learned from your experiment, and I wager it is nothing except that you probably had a bug in your code. I'm not saying that you did, but you should have learned something from your experiment, and if you didn't, then you need to repeat the experiment until you do. If you are confident that your code is bug-free, I suggest that you repeat your experiment with different step sizes and different initial values. You should be able to find examples where both forward and backward Euler blow up and examples where
stvoffutt

## Homework Statement

Solve the following system for $$0<t<5$$
$$u^\prime = u-e^{-2t} v, u(0) = 1$$
$$v^\prime = u+3v, v(0) = -2$$
using Forward Euler method and implement the numerical scheme into a MATLAB code.

## Homework Equations

Forward Euler : $$\vec x^(\prime)_{n+1} = \vec F(t,\vec x)$$
$$\frac{\vec x_{n+1} - \vec x_n}{k} = \vec F(t,\vec x)$$
$$\vec x_{n+1} = \vec x_n + k\vec F(t,\vec x)$$

## The Attempt at a Solution

This seems simple enough but I am having trouble within MATLAB code. First,
$$\vec F(t,\vec x) = \left(\begin{array}{cc} u-e^{-2t} v\\ u+3v \end{array} \right)$$

So, I implemented this into Matlab and got that both u(t) and v(t) diverge to +inf and -inf respectively. Below is my Matlab code.

Code:
function result = hw_2_2(a,b,k)
T = b-a;
N = T/k;
t = [a:k:T];
y_forward = zeros(2,N+1);
y_forward(:,1) = [1;-2];
% y_backward(:,1) = [1;-2];

for n = 1:N

y_forward(:,n+1) = fun_FE(y_forward(:,n),t(n),k);
% y_backward(:,n+1) = fun_BE(y_backward(:,n),t(n+1),k);

end
u = y_forward(1,:);
v = y_forward(2,:);
figure;
plot(t,u,'-r')
title('forward euler')
figure';
plot(t,v,'-b')

% hold off;
% hold on;

% title('backward euler')
% w= y_backward(1,:)
% z = y_backward(2,:)
% plot(t,w,'-r')
% plot(t,z,'-b')

function solFE = fun_FE(u0,t,k)
A = [1 -exp(-2*t);1 3];
solFE = u0 + k*A*u0;

function solBE = fun_BE(u0,t,k)
solBE = inv(1-k*[1,-exp(-2*t);1,3])*u0;

Eventually I will need to use this code to do the same thing for Backward Euler but I wanted to make sure that this code was OK for Forward Euler method first before I proceeded. The function call that I used is
Code:
hw_2_2(0,5,0.005)
Again, I get solution curves that diverge to + - inf.

Strictly speaking, the following is incorrect:
Code:
T = b-a;
N = T/k;
t = [a:k:T];
That happens to work with a=0, which is what you have.

I can't see what you are doing wrong, but this shouldn't blow up. At what time step is your script blowing up? Identifying this point of divergence might help.

Try using the Matlab debugger. Step through the first few iterations to make sure you are doing things correctly, then set a break point so you can advance to just before the point at which things go awry.

stvoffutt said:

## Homework Statement

Solve the following system for $$0<t<5$$
$$u^\prime = u-e^{-2t} v, u(0) = 1$$
$$v^\prime = u+3v, v(0) = -2$$
using Forward Euler method and implement the numerical scheme into a MATLAB code.

## Homework Equations

Forward Euler : $$\vec x^(\prime)_{n+1} = \vec F(t,\vec x)$$
$$\frac{\vec x_{n+1} - \vec x_n}{k} = \vec F(t,\vec x)$$
$$\vec x_{n+1} = \vec x_n + k\vec F(t,\vec x)$$

## The Attempt at a Solution

This seems simple enough but I am having trouble within MATLAB code. First,
$$\vec F(t,\vec x) = \left(\begin{array}{cc} u-e^{-2t} v\\ u+3v \end{array} \right)$$

So, I implemented this into Matlab and got that both u(t) and v(t) diverge to +inf and -inf respectively. Below is my Matlab code.

Code:
function result = hw_2_2(a,b,k)
T = b-a;
N = T/k;
t = [a:k:T];
y_forward = zeros(2,N+1);
y_forward(:,1) = [1;-2];
% y_backward(:,1) = [1;-2];

for n = 1:N

y_forward(:,n+1) = fun_FE(y_forward(:,n),t(n),k);
% y_backward(:,n+1) = fun_BE(y_backward(:,n),t(n+1),k);

end
u = y_forward(1,:);
v = y_forward(2,:);
figure;
plot(t,u,'-r')
title('forward euler')
figure';
plot(t,v,'-b')

% hold off;
% hold on;

% title('backward euler')
% w= y_backward(1,:)
% z = y_backward(2,:)
% plot(t,w,'-r')
% plot(t,z,'-b')

function solFE = fun_FE(u0,t,k)
A = [1 -exp(-2*t);1 3];
solFE = u0 + k*A*u0;

function solBE = fun_BE(u0,t,k)
solBE = inv(1-k*[1,-exp(-2*t);1,3])*u0;

Eventually I will need to use this code to do the same thing for Backward Euler but I wanted to make sure that this code was OK for Forward Euler method first before I proceeded. The function call that I used is
Code:
hw_2_2(0,5,0.005)
Again, I get solution curves that diverge to + - inf.

I do not have access to Matlab, so I cannot comment directly on your code, etc. However, the forward Euler method is known to be unstable, so could easily give you trouble. The exact solution grows very large very quickly, so an unstable method might well diverge.

When I solve the system in closed-form using Maple 14, I get solutions in terms of Bessel Y and J functions. As a matter of possible interest, the numerical values of u(t) and v(t) at t = 5 are
$$u(5) = 1051.1258341055402991 \\ v(5) =-3564922.0735228273290$$
So, v starts at -2 at t = 0 and ends at about -3,000,000 by t = 5.

Ray Vickson said:
I do not have access to Matlab, so I cannot comment directly on your code, etc. However, the forward Euler method is known to be unstable, so could easily give you trouble. The exact solution grows very large very quickly, so an unstable method might well diverge.
Euler's method, done properly, doesn't blow up (i.e., go to infinity), at least not over the interval [0,5]. In this case, Euler's method underestimates in an absolute sense.

D H said:
Euler's method, done properly, doesn't blow up (i.e., go to infinity), at least not over the interval [0,5]. In this case, Euler's method underestimates in an absolute sense.

Agreed---if the arithmetic is exact and free of roundoff errors (which is never). Even if the solution does not go to ##\pm \infty## the computed values may be essentially useless---again because forward Euler + roundoff errors are a bad combination. Theoretical underestimates can, possibly, turn into actual overestimates when roundoff errors are thrown into the mix. It really is hard to say what will happen when one uses theoretically unstable methods, like forward Euler.

The definition of stability is a bit artificial: It's how a technique fares against ##\dot x = -x##. Forward Euler has problems with exponential decay if the step size is too large while backward Euler is stable. That a technique is stable against an exponential decay doesn't necessarily mean it will fare well against some other problem. This problem is exemplary. It's backward Euler that has a much greater tendency to blow up (go to infinity in finite time) than does forward Euler.

Both techniques give bounded values for this problem (advancing state from t=0 to t=5) and using 1000 steps. The result isn't particularly accurate from either technique. You need to use better scheme than either forward Euler or backward Euler on this problem.

Thank you everyone for your replies. I suppose our instructor just gave this exercise to us to illustrate the instability of forward/backward Euler method for this particular problem. I tested my program on another system which gave me the correct answer so I know it's not my code. What conditions caused this problem to be unstable? The eigenvalues of the system?

I implemented both backward and forward Euler for this problem in other languages. I did not see your problem of solutions tending to infinity with forward Euler over the interval [0,5]. You didn't answer my question: Where in [0,5] does the solution first go to infinity? I suspect it might be the very last step. You did something that is always a bit suspect, which is to compute the number of steps as range/step_size. It's much better to compute the step size as range/number_steps. The problem is that 5.0/0.005 is not 1000 thanks to the double precision representation used by default in Matlab.

I calculated the num_steps the way I did because I usually need to change the step_size for one problem. For example I would solve a problem with step size k = 1/4,1/8,1/16,... etc. I ran the program and Matlab does not give INF in any spaces within the solution vectors but the v vector is scaled by 10^6 so the answers are very large which is always a "red flag" in my mind.

So in general it is always best use num_steps as an input and calculate the step_size? I'm still learning Matlab and there are a lot of small things that get over looked in my math classes (such as double precision).

stvoffutt said:
I ran the program and Matlab does not give INF in any spaces within the solution vectors but the v vector is scaled by 10^6 so the answers are very large which is always a "red flag" in my mind.
It's not a red flag in this case. As Ray Vickson noted earlier, the exact solution for v(5) is about -3.56 million. These integrals grow very large in short time, particularly so for v.

So in general it is always best use num_steps as an input and calculate the step_size? I'm still learning Matlab and there are a lot of small things that get over looked in my math classes (such as double precision).
Yep. There are lots of surprising, nasty outcomes when you use doubles. Consider the following:
Code:
N = 49;
stepsize = 1.0/N;
should_be_exactly_one = N*stepsize;
disp (1.0 - should_be_exactly_one);
should_be_exactly_one = sum(step_size*ones(N));
disp(should_be_exactly_one - 1.0);
Note: I have not tested the above with Matlab. I have tested the equivalent in a number of other languages and they all say that ##49 \frac 1 {49} < 1## (by about 1.1e-16) but that ##\sum_{n=1}^{49} \frac 1 {49} > 1## (by about 6.7e-16).

stvoffutt said:
So in general it is always best use num_steps as an input and calculate the step_size? I'm still learning Matlab and there are a lot of small things that get over looked in my math classes (such as double precision).

The general take-home messages (for any programming language) are

1. Integers have exact values, and floating point variables do not, in general.
2. Errors never cancel out when you hope they should.

So, it's better to set num_steps = 2, 4, 8, etc, calculate the step size directly from (b-a)/num_steps. Then do exactly the correct number of steps, rather than testing t against b.

For the same reason, if you are doing a large number of equal sized steps, it is better to calculate the time at the end of the i'th step by t = i*(b-a)/num_steps, rather than repeatedly doing t = t + step_size and letting the error in t accumulate.

If you use the array operations in Matlab, it's reasonable to assume they do the arithmetic the "best" way to minimize this type of error - but if it doubt, and the documentation doesn't answer the question, you can always do a test yourself, similar to D.H's post.

## 1. What is the Forward Euler Method for ODE system?

The Forward Euler Method is a numerical method used to solve ordinary differential equations (ODEs). It is an explicit method, meaning that the solution at a given time step depends only on the solution at the previous time step. It is a first-order method, which means that the error in the solution is proportional to the size of the time step.

## 2. How does the Forward Euler Method work?

The Forward Euler Method works by approximating the solution to an ODE at discrete time steps. It starts with an initial value for the solution and then uses the derivative of the solution to calculate the solution at the next time step. It repeats this process for all time steps until the desired solution is reached.

## 3. What are the advantages of using the Forward Euler Method?

The Forward Euler Method is relatively simple to implement and computationally efficient. It is also useful for solving linear ODEs or for obtaining a rough estimate of the solution to a nonlinear ODE. Additionally, it can handle stiff systems of equations and can be used for adaptive time stepping.

## 4. What are the limitations of the Forward Euler Method?

One of the main limitations of the Forward Euler Method is that it is a first-order method, which means that the error in the solution can accumulate over time. This can lead to inaccurate results, particularly for nonlinear ODEs. It is also not suitable for solving ODEs with discontinuities or for long-term simulations.

## 5. How does the Forward Euler Method compare to other numerical methods for ODEs?

The Forward Euler Method is a relatively simple and straightforward method, but it is not always the most accurate or efficient. Other methods, such as the Runge-Kutta method or the Adams-Bashforth method, may provide more accurate results for certain types of ODEs. It is important to consider the specific characteristics of the ODE system when choosing a numerical method for solving it.

Replies
7
Views
1K
Replies
24
Views
3K
Replies
2
Views
2K
Replies
6
Views
1K
Replies
1
Views
1K
Replies
6
Views
1K
Replies
4
Views
3K
Replies
1
Views
3K
Replies
1
Views
2K
Replies
1
Views
4K