Forward Euler Method for ODE system

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

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

hw_2_2(0,5,0.005)

Again, I get solution curves that diverge to + - inf.

 

[D H]

Strictly speaking, the following is incorrect:

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.

 

[Ray Vickson]

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.

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

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.

 

[D H]

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.

 

[Ray Vickson]

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.

 

[D H]

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.

 

[stvoffutt]

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?

 

[D H]

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.

 

[stvoffutt]

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

 

[D H]

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:

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

 

[AlephZero]

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.