Eulers method to Huegns Method

Question

0 votos

Hello everyone, I am trying to get my former Eulers code to now represent Huen's method. I am having trouble doing this as it seems its giving me the same answers as the Eulers, and it is needing to be more accurate. I have everything commented out as the 'end' in the for loop stops the display of code in matlab central.

if true
  % code
% % This is the work of Dr. Gretarson
% % Damped Simple Harmonic Oscillator solved by Euler's Method
% %
% % This routine implements a numerical (Euler's Method) solution to the
% % damped simple harmonic oscillator.
% % Physical Parameters of the Oscillator
% w0=10; beta=1;
% w=sqrt(w0^2-(beta/2)^2);
% phi=atan(-1/20);
% x0=1/cos(phi);
% % Parameters for the Integration
% h = 0.01; % Time step in seconds
% duration=15; % Duration of solution (seconds)
% Nstep=round(duration/h); % Resulting number of steps needed.
% t=0:h:Nstep*h; % Resulting time vector (seconds)
% % Reserve some memory for the solution vector (these lines are optional but
% % do speed up the code.)
% x=zeros(size(t)); % x is position (meters)
% y=zeros(size(t)); % y is velocity (meters)
% sloperx = zeros(size(t));
% % Set the initial conditions to be the first elements of x and y.
% x(1)=x0; % Initial position (meters)
% y(1)=0; % Initial velocity (meters/second)
% % Iterate to find solution with Euler's method
% for n=2:Nstep
% 
%  f=y(n-1); % f(x,y,t) is dx/dt
%  x(n)=x(n-1)+f*h; % Use the last position to find the new position
%  g=-beta*y(n-1)-w0^2*x(n-1); % g(x,y,t) is dy/dt.
%  y(n)=y(n-1)+g*h; % Use the last velocity to find the new velocity
%  
%  slopery(n) = (y(n)+h*(x(n)));
%  sloperx(n) = (x(n)+h); 
%  
%  
%  ideal(n) = ((1/2)*(sloperx(n) + slopery(n)));
%  
% 
% end
% % The exact solution (solved on paper) is
% xexact=x(1)*exp(-beta/2*t).*cos(w*t+phi);
% yexact=x(1)*(-beta/2*exp(-beta/2*t).*cos(w*t+phi)-w*exp(-beta/2*t).*sin(w*t+phi));
% GDE=max(xexact-x);
% % Compare the exact solution and the numerical solution
% figure(1);
% plot(t,x,'bo',t,xexact,'r','linewidth',1,'markersize',2); hold on;
% title(['GDE = ',num2str(GDE)]);
% grid on; hold off;
% % =======================================================================
% figure(2)
% plot(t,sloperx,'bo',t,xexact,'r')
end