Heun's method stiff ODE

% solving the equation analyticaly

clear all

clc

s = dsolve('Dx = -2.3*x','x(0) = 1', 't');

sol = simplify(s)

% plot the analytical solution in the time interval

t0 = 0;

h = 0.01;

tend = 5;

t = [t0:h:tend];

sol = exp(-(2.3*t));

plot (t, sol, 'k-')

hold on;

% solving the equation numerically using explicit Euler's method

x0 = 1;

h1 = 1;

N = (tend - t0)/h1;

t = [t0:h1:tend];

x(1) = x0;

for i = 1:N

x(i+1) = x(i)+0.5*(h1*-2.3*x(i)+h1*-2.3*(h1*-2.3*x(i)));

end

plot(t,x,'o--')

title('Figure 3.5: Stiff ODE with Euler method when h = 1')

legend({'exact', 'Euler'},'location','southeast')

xlabel('t')

ylabel('x')

ax = gca;

ax.XTick = 0:1:5;

ax.XAxisLocation = 'origin';

ax.YAxisLocation = 'origin';

box off;

hold off;

Is there something wrong with my code??

As the correct graph is not showing.

[John D'Errico]

What did result? Do you recognize that since the problem is stiff, and that an explicit Euler method is often NOT expected to produce a viable result on stiff problems. That is why you are learning about stiff problems.

[Jan]

As far as I see, this comment is missleading:

% solving the equation numerically using explicit Euler's method

[Paul]

Is there a reference for the equation

x(i+1) = ...

Also, what makes this first order differential equation with exponential solution stiff?

Answer 1

[Paul]

Correct implementation of Heun's method, as I understand it, is given below and compare to the analytic solution and other equations given in this thread:

xdot = @(t,x)(-2.3*x);

c1 = 0; c2 = 1;
a11 = 0; a12 = 0; a21 = 1; a22 = 0;
b1 = 1/2; b2 = 1/2;

h1 = 0.25;

t  = 0:h1:5;
xc = 0*t;
x  = 0*t;
y  = 0*t;

xc(1) = 1;
x(1)  = 1;
y(1)  = 1;

for ii = 1:(numel(t)-1)
    k1 = xdot(t(ii),xc(ii));
    k2 = xdot(t(ii)+c2*h1,xc(ii)+h1*a21*k1);
    xc(ii+1) = xc(ii) + h1*(b1*k1 + b2*k2);                      
    x(ii+1) = x(ii)+0.5*(h1*-2.3*x(ii)+h1*-2.3*(h1*-2.3*x(ii))); %Original question
    y(ii+1) = y(ii)+0.5*h1*(-2.3*y(ii)+y(ii)+-2.3*h1*y(ii));     %Answer1
end

plot(t,exp(-2.3*t),'-x',t,x,'-o',t,y,'-o',t,xc,'-o'),grid
legend('analytic','xc','x','y');

Answer 2

[Sulaymon Eshkabilov]

There are a few crucial errs in your code. Here is the completely corrected code. Note that the step size of h1 = 1 gives very crude numerical soutions and thus, h1 has to be chosen carefully.

clearvars; clc                        % It is better start with clearvars instead of clear all that takes some time
syms x(t)                             % This one is recommended syntax of Symbolic Math toolbox     
Dx = diff(x, t);                      % This one is recommended syntax
s = dsolve(Dx == -2.3*x,x(0) == 1);   % This one is recommended syntax
sol = simplify(s)
% plot the analytical solution in the time interval
t0 = 0;
h = 0.01;                             % This is a reasonably small step size  
tend = 5;
t = t0:h:tend;
sol = subs(sol, t);                   % Use the above computed symbolic solution 
plot(t, sol, 'k-')
hold on
% solving the equation numerically using explicit Euler's method
x0 = 1;
h1 = 0.25;                        % Step size needs to be selected carefully
t = t0:h1:tend;                   
y = [x0, zeros(1, numel(t)-1)];   % Memory allocation
for i = 1:numel(t)-1
    y(i+1) = y(i)+0.5*h1*(-2.3*y(i)+y(i)+-2.3*h1*y(i));   % See the corrected loop
end
plot(t,y,'o--')
title('Figure 3.5: Stiff ODE with Euler method when h = 1')
legend({'exact', 'Euler'},'location','southeast')
xlabel('t')
ylabel('x')
ax = gca;
ax.XTick = 0:1:5;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
box off;
hold off;