Download
<< Back
OutputCode
%% Question 2 (Modified ode12)
%% Test cases
% Compare ode12 (Euler & Huen) to ode12_modified (Euler & Runge)
function q2_solution
[t,x] = ode12(@linsystem, 0, 0.1, [2 1]);
[t1,x1] = ode12_modified(@linsystem, 0, 0.1, [2 1]);
plot(t, x(:,1), '-ms', 'DisplayName', 'x_1 (ode12)');
hold on;
plot(t, x(:,2), '-ks', 'DisplayName', 'x_2 (ode12)');
plot(t1, x1(:,1), '-bx', 'DisplayName', 'x_1 (ode12\_modified)');
plot(t1, x1(:,2), '-rx', 'DisplayName', 'x_2 (ode12\_modified)');
legend('Location', 'east');
xlabel('t');
%% Modified ode12 - Euler ($p=1$) and Runge ($p+1=2$)
function [t,x] = ode12_modified(field, t0, T, x0)
t = [t0];
x = zeros(1, length(x0));
x(1,:) = x0;
hmax = (T - t0)/16; % change from ode12_2
h = hmax/8; % change from ode12_2
tol = 0.001;
power = 1/2; % = 1/(p+1)
while t(end) <= T
k1 = feval(field, t(end), x(end,:).');
k2 = feval(field, t(end)+h/2, x(end,:).'+h*k1/2);
% Euler: x(i,:).' + h*k1;
% Runge: x(i,:).' + h*k2;
% => Runge-Euler = h*(k2-k1);
delta = norm(h*(k2-k1), 'inf');
if delta <= tol
% delta less than tolerance, so next timepoint and solution
% computed by Runge
t = [t;t(end)+h];
x = [x; x(end,:)+h*k2.'];
end
if delta ~= 0.0
h = min(hmax, 0.9*h*(tol/delta)^power);
end
end
end
end
