|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
近期用一个方法(增量谐波平衡法)求解非线性微分方程,用duffing方程做例子,很容易就收敛了,但是用到自己的系统的时候,就得不到收敛解,请高手帮忙解决一下这个问题,或者给一些迭代求解线性方程组的建议,不胜感谢啊!
%------------------ parameter definition
m=240;
k0=-2316.4;
k1=12394;
k2=-73696;
k3=3170400;
c0=1385.4;
c1=524.28;
B=0.05;
w1=8;
%------------------方程中各参数求解
w0=sqrt(k1/m);
d1=c0/sqrt(k1*m);
d2=c1*B/m;
s1=1;
s2=k2*B/k1;
s3=k3*B^2/k1;
f=(w1/w0)^2;
ww=w1/w0;
g=k0/(B*k1);
%------------------
syms t
eps=1e-3;
delta_w=0;
w=1;
%%*--------------------- 取2次谐波 A=[a1 a2 b1 b2]' Cs=[1 cos(t)sin(t) sin(2*t)]
A=[0 0 0.05 0]';
Cs=[1 cos(t) sin(t) sin(2*t)];
diff1_Cs=diff(Cs);
diff2_Cs=diff(Cs,2);
i=0;
%----------------------------------
while 1
%------------------------------------------------- 代数方程组参数求解
x0=Cs*A;
diff1_x0=diff1_Cs*A;
diff2_x0=diff2_Cs*A;
Kmc_int=w^2*Cs'*m*diff2_Cs+w*Cs'*(d1+2*w*d2*diff1_x0)*diff1_Cs+Cs'*(s1+2*s2*x0+s3*x0^2)*Cs;
R_int=Cs'*(f*sin(t)-g-(w^2*diff2_x0+d1*w*diff1_x0+d2*w^2*diff1_x0^2+s1*x0+s2*x0^2+s3*x0^3));
Rmc_int=2*w*diff2_x0+2*w*d2*diff1_x0^2+d1*diff1_x0;
Kmc=subs(int(Kmc_int,0,2*pi));
R=subs(int(R_int,0,2*pi));
Rmc=subs(int(Rmc_int,0,2*pi));
%上述为代数方程组中各参数的求解语句,代数方程组为:Kmc*delta_A=R-Rmc*delta_w
%--------------------------------------------------
delta_A=inv(Kmc)*(R-Rmc*delta_w); %此为求解代数方程组的语句
A=A+delta_A;
%-------------------------------------------------
% 跳出循环条件计算
d_R=abs(norm(delta_A)/max(A))
i=i+1
if d_R<=eps
break;
end
end
A
说明:在求解代数方程的时候,由于方程形如AX=b,所以我就直接写成X=inv(A)*b了,不知道合理不合理! |
|