求助，请各位同学、朋友帮帮忙看看这个方程组怎么解

wyllzy

这个方程组解了好久了，一直都没有解决，实在是想不通，
这是第一次用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
先谢谢了

wyllzy

弄了张图片，方便大家看清楚方程

wyllzy

就是这个方程

suffer

用fslove试试看

wyllzy

<P><BR>如果想求出参数随Fa的变化，用fsolve可以么</P>

glise

<P>应该是可以的</P>

suffer

应该可以的，给一个Fa的序列
针对每个Fa的值，调用fsolve求一次