I have symmetric Euler method. Error ratio with Runge - Kuta rule should be approximately 4. But I get 2. The order is p=2.

Symmetric_Euler_Method()

ans = 4×6 table
    N       h       Approximate_Solution    Exact_Solution      Error       Error_Ratio
    __    ______    ____________________    ______________    __________    ___________

    10       0.1          0.54778              0.54784         0.0060561       2.176   
    20      0.05          0.54782              0.54784         0.0027832      2.1003   
    40     0.025          0.54783              0.54784         0.0013251      2.0536   
    80    0.0125          0.54784              0.54784        0.00064526      2.0536   

function results_table = Symmetric_Euler_Method()
y0 = newton(@(y) y - 1, 1, 1e-6);
x0 = 0;
xN = 1;
f = @(x, y) cos(2*x + y) + (3/2)*(x - y);
[x_exact, y_exact] = ode45(f, [0, 1], y0);
N0_values = [10, 20, 40, 80];
p = 2;
results_table = table('Size', [length(N0_values), 6], 'VariableTypes', {'double', 'double', 'double', 'double', 'double', 'double'}, 'VariableNames', {'N', 'h', 'Approximate_Solution', 'Exact_Solution', 'Error', 'Error_Ratio'});
figure;
hold on;
for i = 1:length(N0_values)
    N0 = N0_values(i);
    N = N0;
    h = (xN - x0) / N;
    % Symmetric Euler Method
    t = 0:h:1;
    y = zeros(1, N+1);
    y(1) = y0;
    for j = 1:N
        F = @(y_new) y_new - y(j) - (h/2) * (f(t(j), y(j)) + f(t(j+1), y_new));
        y_new = newton(F, y(j), 1e-6);
        y(j+1) = y_new;
    end
    Yh = y(N+1);
    Y2h = y(N);
    error = abs(Y2h-Yh) / (2^p - 1);
    exact_solution = interp1(x_exact, y_exact, xN);
    if i > 1
        previous_error = results_table.Error(i-1);
        error_ratio = previous_error/error;
    else
        error_ratio = NaN;
    end
    results_table(i, :) = {N, h, Yh, exact_solution, error, error_ratio};
    
    if i == 1
        plot(x_exact, y_exact, 'r', 'DisplayName', 'Exact Solution');
    end
    plot(t, y, 'DisplayName', ['Approximate Solution, h = ' num2str(h)]);
end

% Calculating error ratio with one previous error
for k = 1:length(N0_values) - 1
    error_ratio = results_table.Error(k) / results_table.Error(k+1);
    results_table.Error_Ratio(k) = error_ratio;
end
xlabel('x');
ylabel('y');
title('Solutions with Different N');
legend;
end
function result = newton(F, x0, tol)
max_iter = 100;
x = x0;
iter = 0;

while iter < max_iter
    f_val = F(x);
    f_prime = (F(x + tol) - f_val) / tol;
    x = x - f_val / f_prime;
    
    if abs(f_val) < tol
        result = x;
        return;
    end
    
    iter = iter + 1;
end

error('Newton''s method did not converge');
end

[Torsten]

I don't understand why you define the error of the integration as

Yh = y(N+1);
Y2h = y(N);
error = abs(Y2h-Yh) / (2^p - 1);

I think you must use the difference of the results with stepsize h and h/2 in x=1, not the difference of the last time step with stepsize h.

[E]

How can I fix it? To count error ratio for every h step?

[Torsten]

If you do it like in your other two codes you posted in the forum, it will work.

Answer 1

[Shivansh]

Hi E,

It seems like you are using Symmetric Euler method and the expected error ratio with Runge - Kuta rule should be approximately 4. The output of expected error from the code snippet provided by you is close to 2 in all the cases. As mentioned in the comments, issue with the error ratio being approximately 2 instead of 4 seems to stem from the way the error is calculated.

The error should be calculated using the results obtained with different step sizes, not just the difference between the last two steps of the same run. The error estimation should compare the solutions at the same point ( x = x_N ) for two different step sizes, typically ( h ) and ( h/2 ), and use Richardson's extrapolation to estimate the error.

Here's how you can correct the error calculation:

1. Calculate the solution ( Y_h ) using the step size ( h ).

2. Calculate the solution ( Y_{h/2} ) using the step size ( h/2 ).

3. Estimate the error using the difference between ( Y_h ) and ( Y_{h/2} ), scaled by ( 2^p - 1 ).

Here is the updated code snippet:

Symmetric_Euler_Method()

ans = 4×6 table
    N       h       Approximate_Solution    Exact_Solution      Error       Error_Ratio
    __    ______    ____________________    ______________    __________    ___________

    10       0.1          0.54778              0.54784               NaN         NaN   
    20      0.05          0.54782              0.54784        3.5931e-06      4.0238   
    40     0.025          0.54783              0.54784        8.9297e-07      4.0059   
    80    0.0125          0.54784              0.54784        2.2291e-07      4.0059   

function results_table = Symmetric_Euler_Method()
y0 = newton(@(y) y - 1, 1, 1e-6);
x0 = 0;
xN = 1;
f = @(x, y) cos(2*x + y) + (3/2)*(x - y);
[x_exact, y_exact] = ode45(f, [0, 1], y0);
N0_values = [10, 20, 40, 80];
p = 2;
results_table = table('Size', [length(N0_values), 6], 'VariableTypes', {'double', 'double', 'double', 'double', 'double', 'double'}, 'VariableNames', {'N', 'h', 'Approximate_Solution', 'Exact_Solution', 'Error', 'Error_Ratio'});
figure;
hold on;
for i = 1:length(N0_values)
    N = N0_values(i);
    h = (xN - x0) / N;
    t = x0:h:(x0 + N*h);
    
    % Symmetric Euler Method with step size h
    y = symmetric_euler(f, x0, y0, h, N);
    Yh = y(end);  % Solution with step size h
    
    % Symmetric Euler Method with step size h/2 (if not the first iteration)
    if i > 1
        N_half = N * 2;
        h_half = (xN - x0) / N_half;
        y_half = symmetric_euler(f, x0, y0, h_half, N_half);
        Yh_half = y_half(end);  % Solution with step size h/2
        
        % Error estimation using Richardson's extrapolation
        error = abs(Yh - Yh_half) / (2^p - 1);
        previous_error = results_table.Error(i-1);
        error_ratio = previous_error / error;
    else
        error = NaN;
        error_ratio = NaN;
    end
    
    exact_solution = interp1(x_exact, y_exact(:,1), xN);
    results_table(i, :) = {N, h, Yh, exact_solution, error, error_ratio};
    
    
    if i == 1
        plot(x_exact, y_exact, 'r', 'DisplayName', 'Exact Solution');
    end
    plot(t, y, 'DisplayName', ['Approximate Solution, h = ' num2str(h)]);
end
for k = 1:length(N0_values) - 1
    error_ratio = results_table.Error(k) / results_table.Error(k+1);
    results_table.Error_Ratio(k) = error_ratio;
end
xlabel('x');
ylabel('y');
title('Solutions with Different N');
legend;
end
function result = newton(F, x0, tol)
max_iter = 100;
x = x0;
iter = 0;
while iter < max_iter
    f_val = F(x);
    f_prime = (F(x + tol) - f_val) / tol;
    x = x - f_val / f_prime;
    
    if abs(f_val) < tol
        result = x;
        return;
    end
    
    iter = iter + 1;
end
error('Newton''s method did not converge');
end

function y = symmetric_euler(f, x0, y0, h, N)
y = zeros(1, N+1);
y(1) = y0;
for j = 1:N
    t_j = x0 + (j-1)*h; % Calculate current time step
    t_jp1 = t_j + h;    % Calculate next time step
    F = @(y_new) y_new - y(j) - (h/2) * (f(t_j, y(j)) + f(t_jp1, y_new));
    y_new = newton(F, y(j), 1e-6);
    y(j+1) = y_new;
end
end

I defined a new function "symmetric_euler" that performs the Symmetric Euler method and returns the array of ( y ) values. We use this function to calculate ( Y_h ) and ( Y_{h/2} ). The "interp1" function is used to obtain the exact solution at ( x = x_N ) from the ode45 results.

The error ratio column shows values that are close to 4, which is consistent with the expected behavior of the method. The NaN in the first row for the Error and Error_Ratio columns is expected because there is no previous step size to compare with for the first calculation.

You can refer to the following documentation to learn more about the "interp1" function: https://www.mathworks.com/help/matlab/ref/interp1.html.

I hope it helps!