Euler's Method In Matlab

Question

I am working on a problem involves my using the Euler Method to approximate the differential equation df/dt= af(t)−b[f(t)]^2, both when b=0 and when b is not zero; and I am to compare the analytic solution to the approximate solution when b=0.

f(1) = 1000;
t(1)= 0;
a = 10;
b = 0 ;
dt = 0.01;
Nsteps = 10/dt;

for i = 2:Nsteps
    t(i) = dt + t(i-1);
    %f(i) = f(i-1)*(1 + dt*(a - b*f(i-1)));
    f(i) = f(i-1)*(1 + a*dt); 
end

plot(t,f,'r-')

hold on

fa= a*exp(a*t)

plot(t,fa,'bo')

When b=0, the solution to the differential equation is f(t)=c*exp(at). When I apply the initial condition, that f(0) = 1000, then the differential equation becomes f(t)=1000*exp(at). Now, my professor said that a differential equation has an analytic solution, no matter what time step you use, the graph of analytic solution and the approximation (Euler's Method) will coincide. So, I expected the two graphs to overlap. I attached a picture of what I got.

Why did this occur? In order to get the graphs to overlap, I changed 1000 to 10, which is a=10, just for the heck of it. When I did this, the two overlapped. I don't understand. What am I doing incorrectly?

Answer 1

Why should the numerical solution give the same answer as the analytical one? Looking at pixels overlapping on the screen is not a very precise way to discern anything. You should examine the error between the two (absolute and/or relative). You might also want to examine what happens when you change the step size. And you might want to play with a linear system as well. You don't need to integrate out so far to see these effects – just setting t equal 0.1 or 1 suffices. Here is some better-formatted code to work with:

t0 = 0;
dt = 0.01;
tf = 0.1;
t = t0:dt:tf;    % No need to integrate t in for loop for fixed time step
lt = length(t);
f = zeros(1,lt); % Pre-allocate f
f0 = 1000;       % Initial condition
f(1) = f0;
a = 10;
for i = 1:lt-1
    f(i+1) = f(i) + a*f(i)*dt;
    %f(i+1) = f(i) + a*dt; % Alternative linear system to try
end

% Analytic solution
fa = f0*exp(a*t);
%fa = f0+a*t; % Alternative linear system to try

figure;
plot(t,f,'r-',t,fa,'bo')

% Plot absolute error
figure;
plot(t,abs(f-fa))

% Plot relative error
figure;
plot(t,abs(f-fa)./fa)

You're also not preallocating any of your arrays which makes you code very inefficient. My code does. Read about that here.

Much more beyond this is really off-topic for this site, which is focussed on programming rather than mathematics. If you really have questions about the numerical details that aren't answered by reading your text book (or the Wikipedia page for the Euler method) then you should ask them at Math.StackExchange.

Answer 2

Numerical methods are not precise and there is always an error between numerical and analytical solution. As Euler's method is first order method, global truncation error is proportional to step of integration step.