不知道為啥連接配接傳不上去
clear all
clc
tic
%定義各參數
syms t
w0=3;
epsR=0.001;
m=[1 0；0 1];
epsilon=0.24;
r=0.4;
delta=0.56;
k=[1+epsilon*r -epsilon*r;-epsilon*r 1+epsilon*delta];
Cs=[cos(t) cos(3*t) sin(t) sin(3*t)];
Cs1=diff(Cs，t，1);
S=[cos(t) cos(3*t) sin(t) sin(3*t) 0 0 0 0;0 0 0 0 cos(t) cos(3*t) sin(t) sin(3*t)];
A1=[1 1 1 1]';
A2=[1 1 1 1]';
A0=[A1；A2];
T1=[eye(4，4) zeros(4，4)];
T2=[zeros(4，4) eye(4，4)];
S2=diff(S，t，2);
fm=inline(S'*m*S2);
M=quadv(fm，0，2*pi);
fk=inline(S'*k*S);
K=quadv(fk，0，2*pi);
S1=diff(S，t，1);
c=[1 0；0 1];
fc=inline(S'*c*S1);
C=quadv(fc，0，2*pi);
c3=diag(S*A0).^2;
fc3=inline(S'*c3*S1);
C3=quadv(fc3，0，2*pi);
k2=diag(S*A0).*diag(S1*A0);
fk2=inline(S'*k2*S);
K2=quadv(fk2，0，2*pi);
%代入推導出的公式
Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;
R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;
Rmc=-(2*w0*M+epsilon*(C3-C))*A0;
%AA首元素已知a1=0.0，求ww
a1=0.0;
%變換矩陣，使ww變量代替a1
Kmc11=-Rmc(:，1);
Kmcr=[Kmc11 Kmc(:，2:size(Kmc，2))];
%求未知變量
AA=inv(Kmcr)*R;
%drtA1(1)
ww=AA(1);
%drtW(1)
%賦予新變量新值
A01=A0+[a1; AA(2:length(A0)，1)];
%A(1)+drtA(1)
% Aw0=AA+A00;
%A1(0)+drtA1(1)=A1(1)
w01=w0+ww;
%W+drtW(1)
n=1;
tol=1;
while tol>epsR
A0=A01;
w0=w01;
c3=diag(S*A0).^2;
fc3=inline(S'*c3*S1);
C3=quadv(fc3，0，2*pi);
k2=diag(S*A0).*diag(S1*A0);
fk2=inline(S'*k2*S);
K2=quadv(fk2，0，2*pi);
%帶入推導出的公式
Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;
R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;
Rmc=-(2*w0*M+epsilon*(C3-C))*A0;
%%%%%
tol=norm(R);
if(n>1000)
disp('疊代步數太多，可能不收斂')
return;
end
Kmc11=-Rmc(:，1);
Kmcr=[Kmc11 Kmc(:，2:size(Kmc，2))];
AA=inv(Kmcr)*R;
ww=AA(1);
%A00=[w0；A0(2:6，1)];
A01=A0+[a1；AA(2:length(A0)，1)];
w01=w0+ww;
n=n+1;
end
X0=S*A0;
dX0=S1*A0;
%繪範德波圖
tt=0:.1:10;
xo1=subs(X0(1)，tt);
xo2=subs(X0(2)，tt);
dxo1=subs(dX0(1)，tt);
dxo2=subs(dX0(2)，tt);
figure(1)
plot(xo1，dxo1，'b'，'linewidth'，2)
hold on
plot(xo2，dxo2，'b'，'linewidth'，2)
axis([-3 3 -3 3])
title('範德波極限環')
xlabel('x0')
ylabel('dx0')
toc
這是一個典型例子
聲振之家網站copy過來的，裡面有很多大牛，給學習振動、非線性的同學安利一下這個網站‘聲振之家’