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

1 次查看（过去 30 天）
E
E 2024-2-6
回答： Shivansh 2024-2-22
Ran in:
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
  3 个评论
显示 1更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Shivansh
Shivansh 2024-2-22
Ran in:
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!

类别

在 Help Center 和 File Exchange 中查找有关 Graphics Object Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by