clear all; global om0 gamma gamma=0.1; om0=1; % initial y(0) and v(0) y0=1; v0=2; % to check % [ydot,vdot ]=twoD(1,2) %how many periods? nper=5; T=linspace(0,nper*2*pi,nper*40); d= T(2)-T(1); Y=zeros(1,length(T)); V=zeros(1,length(T)); Y(1)=y0; V(1)=v0; for k=2:length(T); y=Y(k-1); v=V(k-1); [ky1,kv1] = twoD( y , v ); [ky2,kv2] = twoD( y+d*ky1/2 , v+d*kv1/2 ); [ky3,kv3] = twoD( y+d*ky2/2 , v+d*kv2/2 ); [ky4,kv4] = twoD( y+d*ky3 , v+d*kv3 ); Y(k) = y + d/6*(ky1+2*ky2+2*ky3+ky4); V(k) = v + d/6*(kv1+2*kv2+2*kv3+kv4); end % i theoritiki lysi omg = sqrt(om0^2-gamma^2); A = -1i*(v0+y0*(gamma+1i*omg))/(2*omg); B = 1i*(v0+y0*(gamma-1i*omg))/(2*omg); Yth = A*(exp(-gamma*T + 1i*omg*T))+B*(exp(-gamma*T - 1i*omg*T)); % i troxia y(t) figure(14);clf; plot(T,Y,'b',T,Yth,'--r','LineWidth',2) hold on; plot(T,0*T,'--k'); hold off; xlabel('t','FontSize',18); ylabel('y(t)','FontSize',18); title(['damped oscillator with \omega_0=' num2str(om0) ' and \gamma=' num2str(gamma)],'FontSize',18); xlim([T(1) T(end)]); h=legend('numerical','theoretical'); maxd = sqrt(y0^2 + om0^2*v0^2); figure(15);clf; plot(Y,V,'k') hold on; plot(Y(1),V(1),'or','MarkerSize',8) hold off; axis square; axis([-1.2 1.2 -1.2 1.2]*maxd) xlabel('y','FontSize',18); ylabel('v','FontSize',18); title('phase space','FontSize',18);