code improvement/optimization
显示 更早的评论
I would like to know if this code can be improved or it well written? any comments will be greatly appreciated. Thanks
u0 = [1; 0; 0];
tspan = [0, 1e-3]; % small final time for explicit solver
tol = 1e-2;
h0 = 1e-6;
% Execute solver
[t, u, err, h_steps, stats] = dormand_prince_solver(@robertson_ode, tspan, u0, tol, h0);
% Part (a): Plot Solutions
figure;
subplot(2,1,1);
plot(t, u(:,1), 'b', t, u(:,3), 'g', 'LineWidth', 1.5);
legend('y1 (Reactant)', 'y3 (Product)');
xlabel('Time t'); ylabel('Concentration');
%title('Robertson Problem: Concentrations y1 and y3');
grid on;
subplot(2,1,2);
semilogx(t, u(:,2), 'r', 'LineWidth', 1.5);
xlabel('Time t (Log Scale)'); ylabel('Concentration y2');
%title('Robertson Problem: Intermediate Species y2');
grid on;
% Part (b): Error and Stepsize Plots
figure;
subplot(2,1,1);
semilogy(t(2:end), err);
xlabel('Time t'); ylabel('||y_n^5 - y_n^4||_{\infty} / h');
%title('Part (b): Estimated Error vs Time');
grid on;
subplot(2,1,2);
semilogy(t(2:end), h_steps);
xlabel('Time t'); ylabel('Stepsize h');
%title('Part (b): Stepsize h vs Time');
grid on;
figure;
subplot(2,1,1);
semilogy(t(2:end), err, 'LineWidth', 1.5);
xlabel('Time t');
ylabel('Error');
title('Estimated Local Error vs Time');
grid on;
subplot(2,1,2);
semilogy(t(2:end), h_steps, 'LineWidth', 1.5);
xlabel('Time t');
ylabel('Stepsize h');
title('Stepsize vs Time');
grid on;
function dydt = robertson_ode(~, y)
% Coefficients from Screenshot 2026-02-28 at 3.38.10 PM.png
alpha = 0.04;
beta = 1e4;
gamma = 3e7;
dydt = [ -alpha*y(1) + beta*y(2)*y(3);
alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)^2;
gamma*y(2)^2 ];
end
function [t_all, u_all, error_history, h_history, stats] = ...
dormand_prince_solver(f, tspan, u0, tol, h0)
% Initialization
t = tspan(1);
u = u0(:);
h = h0;
t_all = t;
u_all = u.';
error_history = [];
h_history = [];
n_accepted = 0;
n_rejected = 0;
% ---- Dormand–Prince 4(5) Coefficients ----
a = [1/5, 0, 0, 0, 0, 0;
3/40, 9/40, 0, 0, 0, 0;
44/45, -56/15, 32/9, 0, 0, 0;
19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0;
9017/3168, -355/33, 46732/5247, 49/176, -5103/18656, 0;
35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];
b4 = [5179/57600, 0, 7571/16695, 393/640, ...
-92097/339200, 187/2100, 1/40];
b5 = [35/384, 0, 500/1113, 125/192, ...
-2187/6784, 11/84, 0];
c = [0, 1/5, 3/10, 4/5, 8/9, 1, 1];
safety = 0.9;
h_min = 1e-12;
h_max = 1.0;
while t < tspan(2)
if h < h_min
error('Step size underflow: problem likely stiff.');
end
if t + h > tspan(2)
h = tspan(2) - t;
end
% ---- Stage calculations ----
k = zeros(length(u), 7);
k(:,1) = f(t, u);
for i = 2:7
k(:,i) = f(t + c(i)*h, ...
u + h * (k(:,1:i-1) * a(i-1,1:i-1)'));
end
% ---- 4th and 5th order solutions ----
y4 = u + h * (k * b4');
y5 = u + h * (k * b5');
error_val = norm(y5 - y4, inf);
if error_val <= tol
t = t + h;
u = y5;
t_all = [t_all; t];
u_all = [u_all; u.'];
error_history = [error_history; error_val];
h_history = [h_history; h];
n_accepted = n_accepted + 1;
else
n_rejected = n_rejected + 1;
end
if error_val == 0
h_new = 2*h;
else
h_new = h * safety * (tol / error_val)^(1/5);
end
h = min(max(h_new, h_min), h_max);
end
stats.accepted = n_accepted;
stats.rejected = n_rejected;
end
1 个评论
Torsten
about 4 hours 前
Did you test the code for a non-stiff system of differential equations against ode45 for a reasonable timespan ? Are the number of steps and the results similar for equally chosen tolerances ?
回答（1 个）
Hi @Cesar,
Thanks for sharing your code — it is clear you have a solid grasp of the method and the structure is easy to follow. That said, I went through it carefully and found a few things worth fixing before you take this further.
The most important issue is that your b4 and b5 labels are swapped. In the published Dormand-Prince tableau, b5 should hold the coefficients [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0] and b4 the others — but in your code it is the other way around. The numerical values are all correct, so your error estimate still works, but the line u = y5 ends up advancing the integration with the lower-accuracy 4th-order solution instead of the 5th-order one. Dormand and Prince specifically designed the method so the 5th-order solution propagates forward — that is the whole point of DOPRI over the earlier Fehlberg method — so this is worth fixing first.
The second issue is your error criterion. You are using norm(y5 - y4, inf) <= tol which is a flat absolute threshold. MATLAB's ode45 uses a mixed relative and absolute formula: sc = AbsTol + RelTol * max(abs(y4), abs(y5)), then normalises the error against sc. In the Robertson problem, y2 drops to around 1e-5, so with your tol = 1e-2 the solver is essentially blind to inaccuracies in that component. This is exactly what Torsten is getting at when he asks you to compare step counts against ode45 — the two solvers would diverge immediately because they are not measuring the same thing when deciding whether to accept a step.
On that note, Robertson is actually the wrong problem to validate this solver on. It is a canonically stiff system — MathWorks uses it as a benchmark for ode15s, not ode45. On your short tspan the step size is being squeezed by stability constraints, not accuracy, so you are not really testing your error control at all. A much better validation test, which is what Torsten is suggesting, is to run your solver against ode45 on something simple and non-stiff like y' = -y or a harmonic oscillator, with equal tolerances, and check that your step counts are in the same ballpark. If they are not, Issues 1 and 2 above are the reason.
A few more things worth tidying up:
Your step size growth is using fixed absolute bounds (h_min and h_max) rather than the relative growth cap that ode45 uses. MATLAB's formula ensures the step never grows by more than a factor of roughly 5 or shrinks below 10% in a single update — your current approach does not have that guardrail.
The FSAL property is not implemented. Dormand-Prince was designed so that the last stage of one step (k7) is reused as the first stage of the next, saving one function evaluation per accepted step — about 14% less work overall. Right now you are recomputing k(:,1) = f(t, u) fresh every iteration.
The Robertson system satisfies y1 + y2 + y3 = 1 at all times. It is worth adding a quick check after each accepted step to make sure sum(u) stays close to 1 — it is a free sanity check and with your current loose error control it could be flagging something.
Finally, your error subplot label says the quantity is divided by h but the code never actually divides by h — either fix the label or fix the plot. Also, Figures 2 and 3 are plotting the same data with only minor cosmetic differences, so those can be merged into one.
To summarise the priority order: fix the coefficient swap first, then the error normalisation, then run Torsten's benchmark on a non-stiff problem. Everything else is secondary to those three.
The foundation is genuinely solid — the Butcher tableau values are correct and the adaptive loop logic is coherent. These are targeted fixes, not a rewrite.
Good luck with it.
类别
在 帮助中心 和 File Exchange 中查找有关 Interpolation 的更多信息
about 15 hours 前
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
