声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1613|回复: 4

[综合讨论] 一个迭代求解的收敛问题

 关闭 [复制链接]
发表于 2008-5-4 15:43 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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了,不知道合理不合理!
回复
分享到:

使用道具 举报

发表于 2008-5-4 16:29 | 显示全部楼层
由于方程形如AX=b,所以我就直接写成X=inv(A)*b了,不知道合理不合理!

似乎写成X=A\b执行效率更高。
 楼主| 发表于 2008-5-4 16:30 | 显示全部楼层
这个我可以试一下,呵呵,不过怎么解决这个迭代不收敛的问题呢?感谢了!
 楼主| 发表于 2008-5-4 21:25 | 显示全部楼层
自己顶起来,继续求助!
 楼主| 发表于 2008-5-5 20:26 | 显示全部楼层
问题已解决,请版主关闭本帖!

有想做增量谐波平衡法的朋友请和我联系,一起探讨学习吧!呵呵!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-5-16 11:37 , Processed in 0.050722 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表