主程序

InitiationVal=[1;0;2;0];
SolveSpan=[0,30];
[T,X]=ode45('weifenequotion',SolveSpan,InitiationVal);
figure(1);
keyboard;
subplot(5,1,1),plot(T,X(:,1),'r'),xlabel('t'),ylabel('x(1)'), grid on
subplot(5,1,2),plot(T,X(:,2),'k'),xlabel('t'),ylabel('x(2)'),grid on
subplot(5,1,3),plot(T,X(:,3),'y'),xlabel('t'),ylabel('x(2)'),grid on
subplot(5,1,4),plot(T,X(:,4),'g'),xlabel('t'),ylabel('x(2)'),grid on
subplot(5,1,3),plot(X(:,1),X(:,2),'g'),xlabel('x(1)'),ylabel('x(2)'),grid on

 

子程序：

function dy=weifenequotion(t,y)
% 1
% dy=[y(2);0.4];
%
% 2
% dy=[y(2)+cos(t);sin(2*t)];

%3
M=1;
Fr=5;
C=300e-3;
kn=4.5829e10;
Gr=2e-5;
j=1:15;%滚动体数目
Omega=8000;%转速
t0=2;
theta=2*pi*(j-1)/15+Omega*t0;
T=(y(1)*cos(theta)+y(3)*sin(theta)-Gr);

Q1=sum((T.^1.5*kn).*(cos(theta)));
Q2=sum((T.^1.5*kn).*(sin(theta)));
dy(1,1)=y(2);
dy(2,1)=1/M*(Fr-C*y(2)-Q1);
dy(3,1)=y(4);
dy(4,1)=1/M*(-C*y(4)-Q2);

% 4
% global mu;
% dy=[y(2);mu*(1-y(1)^2)*y(2)-y(1)];