Heun's method stiff ODE

6 次查看（过去 30 天）
Aneesa Shahbaz
Aneesa Shahbaz 2021-6-11
回答： Paul 2021-6-15
% 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.
  3 个评论
显示 1更早的评论
Jan
Jan 2021-6-11
As far as I see, this comment is missleading:
% solving the equation numerically using explicit Euler's method
Paul
Paul 2021-6-11
Is there a reference for the equation
x(i+1) = ...
Also, what makes this first order differential equation with exponential solution stiff?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Paul
Paul 2021-6-15
Ran in:
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');

更多回答（1 个）

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;

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by