Improved Euler Method with embedded error estimate and variable time step
1 view (last 30 days)
Show older comments
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)
0 Comments
Accepted Answer
Richard Brown
on 5 Apr 2012
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
0 Comments
More Answers (0)
See Also
Categories
Find more on Language Fundamentals in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!