Compare two different size arrays

My code is listed below. I receive an error of Arrays have incompatible sizes for this operation.

Error (line 92)

err_mp = trapz(t, abs(x1_analytical - x1_midpoint));

Id like to compare the values however it seems they both output in different sizes. Im assuming I will have the same problem with the next few lines. How would I fix this? The rest of the code functions fine to my understanding.

% Define the ODEs
f = @(t, x) [-1/2*x(1); 1/2*x(1)-1/4*x(2); 1/4*x(2)-1/6*x(3)];
% Define the initial conditions
x0 = [10; 10; 10];
% Define the time span
tspan = [0 10];

% Analytical solution
t = linspace(0, 10);
x1_analytical = x0(1)*exp(-t/2);
x2_analytical = -2*x0(1)*exp(-t/2) + (x0(2)+2*x0(1)).*exp(-t/4);
x3_analytical = 3/2*x0(1)*exp(-t/2) - (x0(2)+2*x0(1)).*exp(-t/4) + (x0(3)-3/2*x0(1)+3*(x0(2)+2*x0(1)).*exp(-t/6));

% Mid-point method
h = 0.1;
N = (tspan(2) - tspan(1)) / h;
t_mp = linspace(tspan(1), tspan(2), N+1);
x_mp = zeros(length(x0), N+1);
x_mp(:,1) = x0;
for i = 1:N
    k1 = f(t_mp(i), x_mp(:,i));
    k2 = f(t_mp(i)+h/2, x_mp(:,i)+h/2*k1);
    x_mp(:,i+1) = x_mp(:,i) + h*k2;
end
x1_midpoint = x_mp(1,:);
x2_midpoint = x_mp(2,:);
x3_midpoint = x_mp(3,:);


% Heun's predictor-corrector method
t_hpc = linspace(tspan(1), tspan(2), N+1);
x_hpc = zeros(length(x0), N+1);
x_hpc(:,1) = x0;
for i = 1:N
    k1 = f(t_hpc(i), x_hpc(:,i));
    k2 = f(t_hpc(i)+h, x_hpc(:,i)+h*k1);
    x_hpc(:,i+1) = x_hpc(:,i) + h/2*(k1+k2);
end
x1_heuns = x_hpc(1,:);
x2_heuns = x_hpc(2,:);
x3_heuns = x_hpc(3,:);

% RK45 method
[t_rk45, x_rk45] = ode45(f, tspan, x0);
x1_rk45 = x_rk45(:,1);
x2_rk45 = x_rk45(:,2);
x3_rk45 = x_rk45(:,3);

% RK45 adaptive
tol = 1e-5;
options = odeset('RelTol',tol,'AbsTol',tol);
[t_rk45_adap,x_rk45_adap] = ode45(@(t,x) [(-1/2)*x(1); (1/2)*x(1)-(1/4)*x(2); (1/4)*x(2)-(1/6)*x(3)], tspan, x0, options);
x1_rk45_adap = x_rk45_adap(:,1);
x2_rk45_adap = x_rk45_adap(:,2);
x3_rk45_adap = x_rk45_adap(:,3);

% Plot results
figure;
plot(t, x1_analytical, t_mp, x1_midpoint, t_hpc, x1_heuns, t_rk45, x1_rk45, t_rk45_adap, x1_rk45_adap);
title('x1');
legend('Analytical', 'Mid-point', 'Heun''s', 'RK45', 'RK45 adaptive');
xlabel('t');
ylabel('x1');

figure;
plot(t, x2_analytical, t_mp, x2_midpoint, t_hpc, x2_heuns, t_rk45, x2_rk45, t_rk45_adap, x2_rk45_adap);
title('x2');
legend('Analytical', 'Mid-point', 'Heun''s', 'RK45', 'RK45 adaptive');
xlabel('t');
ylabel('x2');

figure;
plot(t, x3_analytical, t_mp, x3_midpoint, t_hpc, x3_heuns, t_rk45, x3_rk45, t_rk45_adap, x3_rk45_adap);
title('x3');
legend('Analytical', 'Mid-point', 'Heun''s', 'RK45', 'RK45 adaptive');
xlabel('t');
ylabel('x3');

% Plot the results
figure;
plot(t, x1_analytical, 'k-', 'LineWidth', 2);
hold on;
plot(t_mp, x1_midpoint, 'r--', 'LineWidth', 1.5);
plot(t_hpc, x1_heuns, 'g:', 'LineWidth', 1.5);
plot(t_rk45, x1_rk45, 'b-.', 'LineWidth', 1.5);
plot(t_rk45_adap, x1_rk45_adap(:,1), 'm-.', 'LineWidth', 1.5);
xlabel('Time');
ylabel('Concentration');
legend('Analytical', 'Mid-point', 'Heun''s', 'RK45', 'RK45 adaptive', 'Location', 'northwest');

% Calculate the cumulative error for each method
err_mp = trapz(t, abs(x1_analytical - x1_midpoint));

Arrays have incompatible sizes for this operation.

err_heuns = trapz(t, abs(x1_analytical - x1_heuns));
err_rk45 = trapz(t_rk45, abs(x1_analytical - x1_rk45));
err_rk45_adap = trapz(t_rk45_adap, abs(x1_analytical - x1_rk45_adap));

% Report the results
fprintf('Cumulative error for Mid-point method: %f\n', err_mp);
fprintf('Cumulative error for Heun''s method: %f\n', err_heuns);
fprintf('Cumulative error for RK45 method: %f\n', err_rk45);
fprintf('Cumulative error for RK45 adaptive method: %f\n', err_rk45_adap);

Answer 1

[Walter Roberson]

Remember when you use a tspan that consists of exactly two elements, that the ode*() functions choose which times to report for. By default they report 4 times whenever there is an "accepted" point -- the accepted point itself and 3 interpolated points.

You should pass in the times that you used for the analytic solution, so that it will report only for those times.

[Paul]

The error shows up for two cases that don't use any ode functions. I suspect the problem is that numel(t) = 100 for the analytical case, but the midpoint case uses 101 time points. I didn't check the code too carefully.