function [ dy ] = harmonic( t,y,m,k,l,c )
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=-(c*y(2)+k*y(1))/m;
end
function [ dy ] = anharmonicA( t,y,m,k,l,c )
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=-(c*y(2)+k*y(1)+l*y(1)^2)/m;
end
function [ dy ] = anharmonicB( t,y,m,k,l,c )
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=-(c*y(2)+k*y(1)+l*y(1)^3)/m;
end
m=1;
k=1;
l=0.2;
c=0;
x0=1;
x01=2.4;
x02=3;
p0=0;
dy=@(t,y) anharmonicA(t,y,m,k,l,c);
dyB=@(t,y) anharmonicB(t,y,m,k,l,c);
dyh=@(t,y) harmonic(t,y,m,k,l,c);
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);
PSRange=[-10 10];
subplot(3,2,1);
[T,Y] = ode45(dy,[0 40],[x0, p0],options);
[Th,Yh] = ode45(dyh,[0 40],[x0, p0],options);
[TB,YB] = ode45(dyB,[0 40],[x0, p0],options);
plot(Th,Yh(:,1), 'b-', T,Y(:,1),'r-', TB,YB(:,1),'g-');
ylim([-1.5 1.5]*x0);
xlabel('time t');
ylabel('position x(t)');
title(['solution for x_0=' num2str(x0) ', p_0=0']);
subplot(3,2,2);
for x0=-5:0.5:10
[Th,Yh] = ode45(dyh,[0 100],[x0, p0],options);
plot(Yh(:,1),Yh(:,2), 'b-');
if x0==0.5
hold on
end
end
hold off
xlim(PSRange);
ylim(PSRange);
title(['phase space plot, harmonic oscillator']);
subplot(3,2,3);
[T,Y] = ode45(dy,[0 40],[x01, p0],options);
[Th,Yh] = ode45(dyh,[0 40],[x01, p0],options);
[TB,YB] = ode45(dyB,[0 40],[x01, p0],options);
plot(Th,Yh(:,1), 'b-', T,Y(:,1),'r-', TB,YB(:,1),'g-');
xlabel('time t');
ylabel('position x(t)');
title(['solution for x_0=' num2str(x01) ', p_0=0']);
ylim([-2 2]*x01);
subplot(3,2,4);
for x0=-5:0.5:10
[T,Y] = ode45(dy,[0 100],[x0, p0],options);
plot(Y(:,1),Y(:,2), 'r-');
if x0==0.5
hold on
end
end
hold off
xlim(PSRange);
ylim(PSRange);
title(['phase space plot, asymmetric force law']);
subplot(3,2,5);
[T,Y] = ode45(dy,[0 40],[x02, p0],options);
[Th,Yh] = ode45(dyh,[0 40],[x02, p0],options);
[TB,YB] = ode45(dyB,[0 40],[x02, p0],options);
plot(Th,Yh(:,1), 'b-', T,Y(:,1),'r-', TB,YB(:,1),'g-');
xlabel('time t');
ylabel('position x(t)');
title(['solution for x_0=' num2str(x02) ', p_0=0']);
ylim([-3 3]*abs(x02));
subplot(3,2,6);
for x0=-5:0.5:10
[TB,YB] = ode45(dyB,[0 100],[x0, p0],options);
plot(YB(:,1),YB(:,2), 'g-');
if x0==0.5
hold on
end
end
hold off
xlim(PSRange);
ylim(PSRange);
xlabel('position x');
ylabel('impulse p');
title(['phase space plot, symmetric force law']);