guidance on code structure improvement
I would like to know if this code can be improved/optimized. is it well-written? Any help will be appreciated. Thanks
function main_solver_optimized()
params.alpha = 0.04; params.beta = 1e4; params.gamma = 3e7;
y0 = [1; 0; 0];
tspan = [0, 1e6];
h0 = 1e-4;
tol = 1e-2;
fprintf(‘Running Part 1: Radau IIA (2-stage vs 3-stage)…n’);
[T1, Y1, s1, r1] = solve_robertson(y0, tspan, h0, tol, params, ‘Radau’);
fprintf(‘Running Part 2: SDIRK 3(4)…n’);
[T2, Y2, s2, r2] = solve_robertson(y0, tspan, h0, tol, params, ‘SDIRK’);
fprintf(‘n================ RESULTS COMPARISON ================n’);
fprintf(‘Method | Successful Steps | Rejected Stepsn’);
fprintf(‘—————————————————-n’);
fprintf(‘Radau IIA | %16d | %14dn’, s1, r1);
fprintf(‘SDIRK 3(4) | %16d | %14dn’, s2, r2);
fprintf(‘====================================================n’);
figure(‘Name’, ‘Robertson Problem Solutions’);
subplot(2,1,1);
semilogx(T1, Y1(:,1), T1, Y1(:,2)*1e4, T1, Y1(:,3), ‘LineWidth’, 1.5);
title(‘Part 1: Radau IIA Method’); xlabel(‘Time (s)’); ylabel(‘Concentration’);
legend(‘y_1’, ‘y_2 times 10^4’, ‘y_3’, ‘Location’, ‘Best’); grid on;
subplot(2,1,2);
semilogx(T2, Y2(:,1), T2, Y2(:,2)*1e4, T2, Y2(:,3), ‘LineWidth’, 1.5);
title(‘Part 2: SDIRK 3(4) Method’); xlabel(‘Time (s)’); ylabel(‘Concentration’);
legend(‘y_1’, ‘y_2 times 10^4’, ‘y_3’, ‘Location’, ‘Best’); grid on;
end
function [T, Y, success, rejected] = solve_robertson(y0, tspan, h, tol, p, method)
Nmax = 1e6;
T = zeros(Nmax,1); Y = zeros(Nmax,3);
t = tspan(1); y = y0;
T(1) = t; Y(1,:) = y’;
stepCount = 1;
success = 0; rejected = 0;
if strcmp(method,’Radau’)
A2 = [5/12, -1/12; 3/4, 1/4]; c2 = [1/3, 1]’;
sq6 = sqrt(6);
A3 = [(88-7*sq6)/360, (296-169*sq6)/1800, (-2+3*sq6)/225; …
(296+169*sq6)/1800, (88+7*sq6)/360, (-2-3*sq6)/225; …
(16-sq6)/36, (16+sq6)/36, 1/9];
c3 = [(4-sq6)/10, (4+sq6)/10, 1]’;
solver_func = @(y,h) solve_radau_adaptive(y,h,A2,c2,A3,c3,p);
order_p = 3;
else
gam = 1/4;
A = [gam, 0, 0, 0, 0;
1/2-gam, gam, 0, 0, 0;
2*gam, 1-4*gam, gam, 0, 0;
1/10, 1/10, 1/10+gam, gam, 0;
1/4, 1/4, 0, 1/4, gam];
c = [1/4, 3/4, 11/20, 1/2, 1]’;
b = [1/4, 1/4, 0, 1/4, 1/4];
bhat = [-3/16, 5/8, 3/16, 1/16, 5/16];
solver_func = @(y,h) solve_sdirk_sequential(y,h,A,b,bhat,c,p);
order_p = 3;
end
while t < tspan(2)
if t + h > tspan(2), h = tspan(2) – t; end
[y_n, y_hat, converged] = solver_func(y,h);
if converged
err = max(abs(y_n – y_hat));
if err <= tol
t = t + h; y = y_n;
stepCount = stepCount + 1;
T(stepCount) = t;
Y(stepCount,:) = y’;
success = success + 1;
h = h * min(5, max(0.1, 0.9*(tol/err)^(1/(order_p+1))));
else
h = h/4; rejected = rejected + 1;
end
else
h = h/4; rejected = rejected + 1;
end
end
T = T(1:stepCount); Y = Y(1:stepCount,:);
end
function [ynext, yhat, conv] = solve_radau_adaptive(y, h, A2, c2, A3, c3, p)
ynext = y + 0;
yhat = y + 0;
conv = true;
end
function [ynext, yhat, conv] = solve_sdirk_sequential(y, h, A, b, bhat, c, p)
s = length(c);
K = zeros(3,s);
conv = true;
for i = 1:s
ki = zeros(3,1);
rhs = y + h*(K(:,1:i-1)*A(i,1:i-1)’);
for iter = 1:5
val = rhs + h*A(i,i)*ki;
f_val = [-p.alpha*val(1) + p.beta*val(2)*val(3); …
p.alpha*val(1) – p.beta*val(2)*val(3) – p.gamma*val(2)^2; …
p.gamma*val(2)^2];
jac = [-p.alpha, p.beta*val(3), p.beta*val(2); …
p.alpha, -p.beta*val(3)-2*p.gamma*val(2), -p.beta*val(2); …
0, 2*p.gamma*val(2), 0];
F = ki – f_val;
DF = eye(3) – h*A(i,i)*jac;
step = DFF;
ki = ki – step;
if max(abs(step)) < 1e-8, break; end
if iter == 5, conv = false; end
end
K(:,i) = ki;
end
ynext = y + h*(K*b’);
yhat = y + h*(K*bhat’);
endI would like to know if this code can be improved/optimized. is it well-written? Any help will be appreciated. Thanks
function main_solver_optimized()
params.alpha = 0.04; params.beta = 1e4; params.gamma = 3e7;
y0 = [1; 0; 0];
tspan = [0, 1e6];
h0 = 1e-4;
tol = 1e-2;
fprintf(‘Running Part 1: Radau IIA (2-stage vs 3-stage)…n’);
[T1, Y1, s1, r1] = solve_robertson(y0, tspan, h0, tol, params, ‘Radau’);
fprintf(‘Running Part 2: SDIRK 3(4)…n’);
[T2, Y2, s2, r2] = solve_robertson(y0, tspan, h0, tol, params, ‘SDIRK’);
fprintf(‘n================ RESULTS COMPARISON ================n’);
fprintf(‘Method | Successful Steps | Rejected Stepsn’);
fprintf(‘—————————————————-n’);
fprintf(‘Radau IIA | %16d | %14dn’, s1, r1);
fprintf(‘SDIRK 3(4) | %16d | %14dn’, s2, r2);
fprintf(‘====================================================n’);
figure(‘Name’, ‘Robertson Problem Solutions’);
subplot(2,1,1);
semilogx(T1, Y1(:,1), T1, Y1(:,2)*1e4, T1, Y1(:,3), ‘LineWidth’, 1.5);
title(‘Part 1: Radau IIA Method’); xlabel(‘Time (s)’); ylabel(‘Concentration’);
legend(‘y_1’, ‘y_2 times 10^4’, ‘y_3’, ‘Location’, ‘Best’); grid on;
subplot(2,1,2);
semilogx(T2, Y2(:,1), T2, Y2(:,2)*1e4, T2, Y2(:,3), ‘LineWidth’, 1.5);
title(‘Part 2: SDIRK 3(4) Method’); xlabel(‘Time (s)’); ylabel(‘Concentration’);
legend(‘y_1’, ‘y_2 times 10^4’, ‘y_3’, ‘Location’, ‘Best’); grid on;
end
function [T, Y, success, rejected] = solve_robertson(y0, tspan, h, tol, p, method)
Nmax = 1e6;
T = zeros(Nmax,1); Y = zeros(Nmax,3);
t = tspan(1); y = y0;
T(1) = t; Y(1,:) = y’;
stepCount = 1;
success = 0; rejected = 0;
if strcmp(method,’Radau’)
A2 = [5/12, -1/12; 3/4, 1/4]; c2 = [1/3, 1]’;
sq6 = sqrt(6);
A3 = [(88-7*sq6)/360, (296-169*sq6)/1800, (-2+3*sq6)/225; …
(296+169*sq6)/1800, (88+7*sq6)/360, (-2-3*sq6)/225; …
(16-sq6)/36, (16+sq6)/36, 1/9];
c3 = [(4-sq6)/10, (4+sq6)/10, 1]’;
solver_func = @(y,h) solve_radau_adaptive(y,h,A2,c2,A3,c3,p);
order_p = 3;
else
gam = 1/4;
A = [gam, 0, 0, 0, 0;
1/2-gam, gam, 0, 0, 0;
2*gam, 1-4*gam, gam, 0, 0;
1/10, 1/10, 1/10+gam, gam, 0;
1/4, 1/4, 0, 1/4, gam];
c = [1/4, 3/4, 11/20, 1/2, 1]’;
b = [1/4, 1/4, 0, 1/4, 1/4];
bhat = [-3/16, 5/8, 3/16, 1/16, 5/16];
solver_func = @(y,h) solve_sdirk_sequential(y,h,A,b,bhat,c,p);
order_p = 3;
end
while t < tspan(2)
if t + h > tspan(2), h = tspan(2) – t; end
[y_n, y_hat, converged] = solver_func(y,h);
if converged
err = max(abs(y_n – y_hat));
if err <= tol
t = t + h; y = y_n;
stepCount = stepCount + 1;
T(stepCount) = t;
Y(stepCount,:) = y’;
success = success + 1;
h = h * min(5, max(0.1, 0.9*(tol/err)^(1/(order_p+1))));
else
h = h/4; rejected = rejected + 1;
end
else
h = h/4; rejected = rejected + 1;
end
end
T = T(1:stepCount); Y = Y(1:stepCount,:);
end
function [ynext, yhat, conv] = solve_radau_adaptive(y, h, A2, c2, A3, c3, p)
ynext = y + 0;
yhat = y + 0;
conv = true;
end
function [ynext, yhat, conv] = solve_sdirk_sequential(y, h, A, b, bhat, c, p)
s = length(c);
K = zeros(3,s);
conv = true;
for i = 1:s
ki = zeros(3,1);
rhs = y + h*(K(:,1:i-1)*A(i,1:i-1)’);
for iter = 1:5
val = rhs + h*A(i,i)*ki;
f_val = [-p.alpha*val(1) + p.beta*val(2)*val(3); …
p.alpha*val(1) – p.beta*val(2)*val(3) – p.gamma*val(2)^2; …
p.gamma*val(2)^2];
jac = [-p.alpha, p.beta*val(3), p.beta*val(2); …
p.alpha, -p.beta*val(3)-2*p.gamma*val(2), -p.beta*val(2); …
0, 2*p.gamma*val(2), 0];
F = ki – f_val;
DF = eye(3) – h*A(i,i)*jac;
step = DFF;
ki = ki – step;
if max(abs(step)) < 1e-8, break; end
if iter == 5, conv = false; end
end
K(:,i) = ki;
end
ynext = y + h*(K*b’);
yhat = y + h*(K*bhat’);
end I would like to know if this code can be improved/optimized. is it well-written? Any help will be appreciated. Thanks
function main_solver_optimized()
params.alpha = 0.04; params.beta = 1e4; params.gamma = 3e7;
y0 = [1; 0; 0];
tspan = [0, 1e6];
h0 = 1e-4;
tol = 1e-2;
fprintf(‘Running Part 1: Radau IIA (2-stage vs 3-stage)…n’);
[T1, Y1, s1, r1] = solve_robertson(y0, tspan, h0, tol, params, ‘Radau’);
fprintf(‘Running Part 2: SDIRK 3(4)…n’);
[T2, Y2, s2, r2] = solve_robertson(y0, tspan, h0, tol, params, ‘SDIRK’);
fprintf(‘n================ RESULTS COMPARISON ================n’);
fprintf(‘Method | Successful Steps | Rejected Stepsn’);
fprintf(‘—————————————————-n’);
fprintf(‘Radau IIA | %16d | %14dn’, s1, r1);
fprintf(‘SDIRK 3(4) | %16d | %14dn’, s2, r2);
fprintf(‘====================================================n’);
figure(‘Name’, ‘Robertson Problem Solutions’);
subplot(2,1,1);
semilogx(T1, Y1(:,1), T1, Y1(:,2)*1e4, T1, Y1(:,3), ‘LineWidth’, 1.5);
title(‘Part 1: Radau IIA Method’); xlabel(‘Time (s)’); ylabel(‘Concentration’);
legend(‘y_1’, ‘y_2 times 10^4’, ‘y_3’, ‘Location’, ‘Best’); grid on;
subplot(2,1,2);
semilogx(T2, Y2(:,1), T2, Y2(:,2)*1e4, T2, Y2(:,3), ‘LineWidth’, 1.5);
title(‘Part 2: SDIRK 3(4) Method’); xlabel(‘Time (s)’); ylabel(‘Concentration’);
legend(‘y_1’, ‘y_2 times 10^4’, ‘y_3’, ‘Location’, ‘Best’); grid on;
end
function [T, Y, success, rejected] = solve_robertson(y0, tspan, h, tol, p, method)
Nmax = 1e6;
T = zeros(Nmax,1); Y = zeros(Nmax,3);
t = tspan(1); y = y0;
T(1) = t; Y(1,:) = y’;
stepCount = 1;
success = 0; rejected = 0;
if strcmp(method,’Radau’)
A2 = [5/12, -1/12; 3/4, 1/4]; c2 = [1/3, 1]’;
sq6 = sqrt(6);
A3 = [(88-7*sq6)/360, (296-169*sq6)/1800, (-2+3*sq6)/225; …
(296+169*sq6)/1800, (88+7*sq6)/360, (-2-3*sq6)/225; …
(16-sq6)/36, (16+sq6)/36, 1/9];
c3 = [(4-sq6)/10, (4+sq6)/10, 1]’;
solver_func = @(y,h) solve_radau_adaptive(y,h,A2,c2,A3,c3,p);
order_p = 3;
else
gam = 1/4;
A = [gam, 0, 0, 0, 0;
1/2-gam, gam, 0, 0, 0;
2*gam, 1-4*gam, gam, 0, 0;
1/10, 1/10, 1/10+gam, gam, 0;
1/4, 1/4, 0, 1/4, gam];
c = [1/4, 3/4, 11/20, 1/2, 1]’;
b = [1/4, 1/4, 0, 1/4, 1/4];
bhat = [-3/16, 5/8, 3/16, 1/16, 5/16];
solver_func = @(y,h) solve_sdirk_sequential(y,h,A,b,bhat,c,p);
order_p = 3;
end
while t < tspan(2)
if t + h > tspan(2), h = tspan(2) – t; end
[y_n, y_hat, converged] = solver_func(y,h);
if converged
err = max(abs(y_n – y_hat));
if err <= tol
t = t + h; y = y_n;
stepCount = stepCount + 1;
T(stepCount) = t;
Y(stepCount,:) = y’;
success = success + 1;
h = h * min(5, max(0.1, 0.9*(tol/err)^(1/(order_p+1))));
else
h = h/4; rejected = rejected + 1;
end
else
h = h/4; rejected = rejected + 1;
end
end
T = T(1:stepCount); Y = Y(1:stepCount,:);
end
function [ynext, yhat, conv] = solve_radau_adaptive(y, h, A2, c2, A3, c3, p)
ynext = y + 0;
yhat = y + 0;
conv = true;
end
function [ynext, yhat, conv] = solve_sdirk_sequential(y, h, A, b, bhat, c, p)
s = length(c);
K = zeros(3,s);
conv = true;
for i = 1:s
ki = zeros(3,1);
rhs = y + h*(K(:,1:i-1)*A(i,1:i-1)’);
for iter = 1:5
val = rhs + h*A(i,i)*ki;
f_val = [-p.alpha*val(1) + p.beta*val(2)*val(3); …
p.alpha*val(1) – p.beta*val(2)*val(3) – p.gamma*val(2)^2; …
p.gamma*val(2)^2];
jac = [-p.alpha, p.beta*val(3), p.beta*val(2); …
p.alpha, -p.beta*val(3)-2*p.gamma*val(2), -p.beta*val(2); …
0, 2*p.gamma*val(2), 0];
F = ki – f_val;
DF = eye(3) – h*A(i,i)*jac;
step = DFF;
ki = ki – step;
if max(abs(step)) < 1e-8, break; end
if iter == 5, conv = false; end
end
K(:,i) = ki;
end
ynext = y + h*(K*b’);
yhat = y + h*(K*bhat’);
end code improvement MATLAB Answers — New Questions
