声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1347|回复: 6

[编程技巧] 求助,请各位同学、朋友帮帮忙看看这个方程组怎么解

[复制链接]
发表于 2006-1-7 15:37 | 显示全部楼层 |阅读模式

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

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

x
这个方程组解了好久了,一直都没有解决,实在是想不通,
这是第一次用matlab编程求解方程组,请各位帮忙看看,我的问题出在哪里,有什么更好的方法解这种方程么,主要是想求出未知数ae,ai两个角随外力Fa的变化关系
f1=Ki*(Si^1.5)*sin(ai)-Ke*(Se^1.5)*sin(ae)-(pi*p*((Db/1000)^4)*(w^2)/(30))*(sin(ae)*cos(ae)/r)*(((1-r*cos(ai))/(1+cos(ai-ae)))^2);
f2=Ki*(Si^1.5)*cos(ai)-Ke*(Se^1.5)*cos(ae)+(p*pi*((Db/1000)^4)*(w^2)/(30*r))*((sin(ae))^2)*(((1-r*cos(ai))/(1+cos(ai-ae)))^2)+(pi*p*((Db/1000)^3)*(Dm/1000)*(w^2)/12)*(((1-r*cos(ai))/(1+cos(ai-ae)))^2);
f3=((fe-0.5)*Db+Se)*sin(ae)+((fi-0.5)*Db+Si)*sin(ai)-(fe+fi-1)*Db*sin(a)-Sa;
f4=((fe-0.5)*Db+Se)*cos(ae)+((fi-0.5)*Db+Si)*cos(ai)-(fe+fi-1)*Db*cos(a);
f5=Fa-20*Ki*(Si^1.5)*sin(ai);
以下是我编的牛顿迭代,请求各位帮帮忙,愁死我了
为什么Fa变化时,五个未知数的值只有Sa一个变化,其他的变化很小,前四个方程的值在零附近,而第五个方程的值随Fa的变化一直增大,
Fa按道理说应该能加到很大的,可是我算的Fa到1000多点时,就不行了
而且有时迭代时,五个未知数中只有一个变化,其他的只在初值附近,变化不大,怎么改进

function x=xiugai( )
各参数初值
fe=0.52;
fi=0.515;
Fa=250;
Ki=381447.0905;
Ke=329428.2257;
a=15*pi/180;
Db=11;
Dm=82;
w=1046;
p=7800;
r=11/82;

%五个未知数的初解
ae=9*pi/180; %角度
ai=24.5*pi/180; %角度
Se=0.0037; %变形
Si=0.00188; %变形
Sa=0.0045; %变形

for j=1:1:90

for i=1:1:100
%J是Jacobi矩阵
y1=(1000*p*pi*((Db/1000)^5)*(w^2)/(30*r*Db))*((((cos(ae))^2)-((sin(ae))^2))*(((1-r*cos(ai))/(1+cos(ai-ae)))^2)-sin(2*ae)*sin(ai-ae)*((1-r*cos(ai))^2)/((1+cos(ai-ae))^3));
y2=(1000*p*pi*((Db/1000)^5)*(w^2)/(30*r*Db))*sin(2*ae)*(r*sin(ai)*(1-r*cos(ai))/((1+cos(ai-ae))^2)+((1-r*cos(ai))^2)*sin(ai-ae)/((1+cos(ai-ae))^3));
y3=(1000*p*pi*((Db/1000)^5)*(w^2)/(30*r*Db))*(sin(2*ae)*((((1-r*cos(ai))^2)/(1+cos(ai-ae)))^2)-2*((sin(ae))^2)*((1-r*cos(ai))^2)*sin(ai-ae)/((1+cos(ai-ae))^3))-(pi*p*((Db/1000)^3)*(Dm/1000)*(w^2)/6)*((1-r*cos(ai))^2)*sin(ai-ae)/((1+cos(ai-ae))^3);
y4=(1000*p*pi*((Db/1000)^5)*(w^2)/(30*r*Db))*((sin(ae))^2)*(2*r*sin(ai)*(1-r*cos(ai))/((1+cos(ai-ae))^2)+2*sin(ai-ae)*((1-r*cos(ai))^2)/((1+cos(ai-ae))^3))+(pi*p*((Db/1000)^3)*(Dm/1000)*(w^2)/6)* (r*sin(ai)*(1-r*cos(ai))/((1+cos(ai-ae))^2)+ sin(ai-ae)*((1-r*cos(ai))^2)/((1+cos(ai-ae))^3));

J=[Ki*(Si^1.5)*cos(ai)-y2 -Ke*(Se^1.5)*cos(ae)-y1 1.5*Ki*(Si^0.5)*sin(ai) -1.5*Ke*(Se^0.5)*sin(ae) 0
-20*Ki*(Si^1.5)*cos(ai) 0 -30*Ki*(Si^0.5)*sin(ai) 0 0
-sin(ai)*((fi-0.5)*Db+Si) -sin(ae)*((fe-0.5)*Db+Se) cos(ai) cos(ae) 0
((fi-0.5)*Db+Si)*cos(ai) ((fe-0.5)*Db+Se)*cos(ae) sin(ai) sin(ae) -1
-Ki*(Si^1.5)*sin(ai)+y4 Ke*(Se^1.5)*sin(ae)+y3 1.5*Ki*(Si^0.5)*cos(ai) 1.5*Ke*(Se^0.5)*cos(ae) 0]
%f1-f5是五个方程
f1=Ki*(Si^1.5)*sin(ai)-Ke*(Se^1.5)*sin(ae)-(pi*p*((Db/1000)^4)*(w^2)/(30))*(sin(ae)*cos(ae)/r)*(((1-r*cos(ai))/(1+cos(ai-ae)))^2);
f2=Ki*(Si^1.5)*cos(ai)-Ke*(Se^1.5)*cos(ae)+(p*pi*((Db/1000)^4)*(w^2)/(30*r))*((sin(ae))^2)*(((1-r*cos(ai))/(1+cos(ai-ae)))^2)+(pi*p*((Db/1000)^3)*(Dm/1000)*(w^2)/12)*(((1-r*cos(ai))/(1+cos(ai-ae)))^2);
f3=((fe-0.5)*Db+Se)*sin(ae)+((fi-0.5)*Db+Si)*sin(ai)-(fe+fi-1)*Db*sin(a)-Sa;
f4=((fe-0.5)*Db+Se)*cos(ae)+((fi-0.5)*Db+Si)*cos(ai)-(fe+fi-1)*Db*cos(a);
f5=Fa-20*Ki*(Si^1.5)*sin(ai);
Fx=[f1 f5 f4 f3 f2]';


%迭代计算
E=eye(5:5);
x=x-inv(J)*Fx;
ai=x(1,1);
ae=x(2,1);
Si=abs(x(3,1));
Se=abs(x(4,1));
Sa=abs(x(5,1));
q=inv(J)*Fx;
if abs(q(1,1))<0.01&abs(q(2,1))<0.01,break,end
end

Fa=Fa+10;

end
先谢谢了
回复
分享到:

使用道具 举报

 楼主| 发表于 2006-1-7 15:57 | 显示全部楼层
弄了张图片,方便大家看清楚方程
 楼主| 发表于 2006-1-7 16:18 | 显示全部楼层

回复:(wyllzy)求助,请各位同学、朋友帮帮忙看看这...

就是这个方程
1.jpg
发表于 2006-1-7 17:46 | 显示全部楼层

回复:(wyllzy)求助,请各位同学、朋友帮帮忙看看这...

用fslove试试看
 楼主| 发表于 2006-1-7 19:46 | 显示全部楼层
<P><BR>如果想求出参数随Fa的变化,用fsolve可以么</P>
发表于 2006-1-8 10:18 | 显示全部楼层

回复:(wyllzy)求助,请各位同学、朋友帮帮忙看看这...

<P>应该是可以的</P>
发表于 2006-1-11 10:04 | 显示全部楼层

回复:(wyllzy)如果想求出参数随Fa的变化,用fsolve...

应该可以的,给一个Fa的序列
针对每个Fa的值,调用fsolve求一次
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-20 23:35 , Processed in 0.062577 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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