MATLAB: Improved Euler Method with embedded error estimate and variable time step

euler

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)
% subplot(4,1,1)
% plot(t,X)
% subplot(4,1,2)
% plot(t,Y)
% subplot(4,1,3)
% plot(t,Z)
% subplot(4,1,4)
% plot(t,E)
% 





% figure

% 
% subplot(3,1,1)

% plot(t,Error_x)
% subplot(3,1,2)

% plot(t,Error_y)
% subplot(3,1,3)

% plot(t,Error_z)
% 
figure
% 
% subplot(3,1,1)
%plot(t,nx)
% subplot(3,1,2)
plot(Error)
% subplot(3,1,3)
% plot(t,nz)
% 
% figure
% 
% plot(X,Y)

Best Answer

you could include a while loop inside your inner loop that's loosely like this:
e_norm = inf;
while e_norm > E_max
   do some stuff
   e_norm = something better, hopefully
   if "I've done this too many times and it's not working"
      error('bad')
   end
end
Related Question