clear all; % ''pagkosmies'' statheres global gamma om gamma = 0.1; om = 1; % arxikes synthikes x0 =1; y0 =0; % orismos ton xronikon stigmon T = linspace(0,30,201); dt = T(2)-T(1) % to dianysma sto opoio tha apothikevo ti lysi x(t),y(t) X = zeros(2,length(T)); % i proti colona tou X einai oi arxikes synthikes X(1,1) = x0; X(2,1) = y0; for it=2:length(T) x = X(1,it-1); y = X(2,it-1); t = T(it-1); % RK-4 [k1x , k1y ] = func2D(x,y,t); [k2x , k2y ] = func2D(x+k1x*dt/2 , y+k1y*dt/2 , t+dt/2); [k3x , k3y ] = func2D(x+k2x*dt/2 , y+k2y*dt/2 , t+dt/2); [k4x , k4y ] = func2D(x+k3x*dt , y+k3y*dt , t+dt); xnew = x + (k1x+2*k2x+2*k3x+k4x)*dt/6; ynew = y + (k1y+2*k2y+2*k3y+k4y)*dt/6; % apothikevo tis nees times x,y stin epomeni kolona tou X X(1,it)=xnew; X(2,it)=ynew; end % i analytiki lysi gia ton talantoti me aposvesi OM = sqrt(om^2-gamma^2); xtheor = x0*exp(-gamma*T).*( cos(OM*T) + gamma/OM*sin(OM*T) ) + y0*exp(-gamma*T)/OM.*sin(OM*T); figure(1) plot(T,X(1,:),'ob','LineWidth',2); hold on; plot(T,xtheor,'--r') hold off; xlabel('t','fontsize',20); ylabel('x','fontsize',20); h=legend('numerical','analytical'); figure(3) plot(X(1,:),X(2,:),'-k') axis square; xlabel('x','fontsize',20); ylabel('xdot','fontsize',20); drawnow; pause(1)