Greetings all,
I am trying to simulate a system of equations using Improved Euler's method. I'm not really sure if I am doing the embedded error estimate or the variable time step correctly. But I want my code to repeat the previous step through my for loop when e_norm > Emax but am unsure how to do this. Any help is appreciated.
e = 0;
x = .1;
y = .1;
z = .1;
T = .1;
N = 30;
a = 15;
Emax = .1;
Emin = .01;
for n = 1:N+1
xdot1 = -a*(-0.07*e+0.085*( abs(e+1) - abs(e-1) ));
ydot1 = -1*(-0.07*e+0.085*( abs(e+1) - abs(e-1) )) - z;
zdot1 = y;
xh = x + T*xdot1;
yh = y + T*ydot1;
zh = z + T*zdot1;
xdot2 = -a*(-0.07*e+0.085*( abs(e+1) - abs(e-1) ));
ydot2 = -1*(-0.07*e+0.085*( abs(e+1) - abs(e-1) )) - zh;
zdot2 = yh;
x = x + .5*T*(xdot2 + xdot1);
y = y + .5*T*(ydot2 + ydot1);
z = z + .5*T*(zdot2 + zdot1);
e = y - x;
X(n) = x;
Y(n) = y;
Z(n) = z;
E(n) = e;
Error_x(n) = x - xh;
Error_y(n) = y - yh;
Error_z(n) = z - zh;
nx = norm(Error_x,1);
ny = norm(Error_y,1);
nz = norm(Error_z,1);
e_norm = sqrt(nx.^2 + ny.^2 + nz.^2);
if e_norm > Emax;
T = T/2;
elseif e_norm > Emin && e_norm < 0.8*Emax;
T = 1.05*T;
elseif e_norm < Emin;
T = 1.25*T;
elseif e_norm > 0.8*Emax && e_norm < Emax;
T = T;
end
timestep(n) = T
Error(n) = e_norm
end
plot(timestep)
figure
plot(Error)
