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
clcclear allclose all%% ----------1-STAGE-EULER----------
n = [500,1000,2000,4000];T = table();for p = 1:4 steps = n(p); T.m(p) = steps; endinit_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; endT.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); endT.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; endT.RK_4(p) = abs(init_cond(1) - z(end,1));clearvars -except T steps n pfigureloglog(T.m,T.Euler_1)hold onloglog(T.m,T.Euler_2)loglog(T.m,T.RK_4)xlabel('n'), ylabel('Error'), grid ontitle('Error Mag. vs. Number of Points (n)')legend('1 Step Euler','2 Step Euler','Runge-Kutta 4')T
Best Answer