How can I get values for n=500, 1000 and 2000 ? only getting values for n=4000!

Question

MG 2020-10-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/616298-how-can-i-get-values-for-n-500-1000-and-2000-only-getting-values-for-n-4000

评论： MG 2020-10-16

采纳的回答： Alan Stevens

在 MATLAB Online 中打开

m Euler_1 Euler_2 RK_4

____ ________ __________ __________

500 0 0 0

1000 0 0 0

2000 0 0 0

4000 0.093706 2.0728e-05 8.9537e-13

clc
clear all
close all
%% ----------1-STAGE-EULER----------
n = [500,1000,2000,4000];
T = table();
for p = 1:4
    
    steps = n(p);
    T.m(p) = steps; 
end
init_cond = [1 1i];
G = 1;
M = 1;
t_f = 2*pi/norm(1i);
dt = t_f/steps;
t = 0:dt:t_f;
z = init_cond;
f = @(t,z) [z(2),-G*M*z(1)/norm(z(1))^3];
for i = 1:steps
    
    k1 = f(t(i),z(i,:));
    z(i+1,:) = z(i,:) + dt.*k1; 
end
T.Euler_1(p) = abs(init_cond(1) - z(end,1));
clearvars -except T steps n p Euler_1
%% ----------2-STAGE-EULER----------
init_cond = [1 1i];
G = 1;
M = 1;
t_f = 2*pi/norm(1i);
dt = t_f/steps;
t = 0:dt:t_f;
z = init_cond;
f = @(t,z) [z(2),-G*M*z(1)/norm(z(1))^3];
for i = 1:steps
    
    k1 = f(t(i),z(i,:));
    k2 = f(t(i)+dt,z(i,:)+dt*k1);
    z(i+1,:) = z(i,:) + dt/2.*(k1 + k2); 
end
T.Euler_2(p) = abs(init_cond(1) - z(end,1));
clearvars -except T steps n p Euler_2
%% ----------4-STAGE-CLASSICAL-RK----------
init_cond = [1 1i];
G = 1;
M = 1;
t_f = 2*pi/norm(1i);
dt = t_f/steps;
t = 0:dt:t_f;
z = init_cond;
f = @(t,z) [z(2),-G*M*z(1)/norm(z(1))^3];
for i = 1:steps
    k1 = f(t(i),z(i,:));
    k2 = f(t(i) + 0.5*dt,z(i,:) + 0.5*dt.*k1);
    k3 = f((t(i) + 0.5*dt),(z(i,:) + 0.5*dt.*k2));
    k4 = f((t(i) + dt),(z(i,:) + k3.*dt));
    z(i+1,:) = z(i,:) + (1/6)*(k1 + 2*k2 + 2*k3 + k4).*dt;
    
end
T.RK_4(p) = abs(init_cond(1) - z(end,1));
clearvars -except T steps n p
figure
loglog(T.m,T.Euler_1)
hold on
loglog(T.m,T.Euler_2)
loglog(T.m,T.RK_4)
xlabel('n'), ylabel('Error'), grid on
title('Error Mag. vs. Number of Points (n)')
legend('1 Step Euler','2 Step Euler','Runge-Kutta 4')
T

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Alan Stevens 2020-10-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/616298-how-can-i-get-values-for-n-500-1000-and-2000-only-getting-values-for-n-4000#answer_515968

在 MATLAB Online 中打开

Your assignments

T.Euler_1(p) = abs(init_cond(1) - z(end,1));

etc. are all outside of the for p = ... loop, hence they just retain the last value of p, namely 4, and only store values in T.Euler_1(4) etc. overwriting anything previously stored there, hence ending with only the n = 4000 values.