作業

clc;close all;clear all;
 
global m c k 
m = 100;
c = 5;
k = 10;
 
[t,y] = ode45(@mck,[0,10*pi],[0,50]); %% ODE主程式
 
dt=0.0001;
n=10*pi/dt;
 
y1(1) = 0;
y2(1) = 50;
time(1) = 0; 
 
for i=1:n
    time(i+1) = dt*i;
    y1(i+1) = dt*y2(i)+y1(i);
    y2(i+1) = dt*((2*sin(time(i))+sin(time(i)-pi/3))/m-(k/m*y1(i))-((c/m)*y2(i)))+y2(i);
end

subplot(2,1,1)
plot(t,y(:,1),'r')
hold on
plot(time,y1,'b')
title('ODE45&差分法(位移)')
xlabel('t(s)')
ylabel('S(m)')
 
subplot(2,1,2)
plot(t,y(:,2),'r')
hold on
plot(time,y2,'b')
title('ODE45&差分法(速度)')
xlabel('t(s)')
ylabel('V(m/s)')
 
function ydot=mck(t,y)
    global m c k 
    u=2*sin(t)+sin(t-pi/3);
    ydot=zeros(2,1);
    ydot(1)=y(2);
    ydot(2) = u/100-c/m*y(2)-k/m*y(1);
end