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

1 view (last 30 days)
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

Accepted Answer

Alan Stevens
Alan Stevens on 16 Oct 2020
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.

More Answers (0)

Categories

Find more on Entering Commands in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!