zzkingg 发表于 2010-11-28 23:37

急求此非线性方程组的解法?谢谢!

h=5;                           
d=4;                        
R=0.7;                        
U=R*2*pi;                     
c=55000;                     
fi=35;                        
K0=0.95-sin(fi*pi/180);      
Es=18000000;                  
Ep=20000000;                  
v=19000;                     
t=v*h*K0*tan(fi*pi/180)+c;   
Ae=d^2;                        
Ap=pi*R^2;                     
As=Ae-Ap;                  
q=2000000;                  
Kvp=150000000;               
Kvs=5000000;                  
ap=v/Ep;
bp=U/(Ep*Ap);                  
as=v/Es;                     
bs=U/(Es*As);                  
m=Ap/Ae;                     

(1-m)*Kvs*X(7)+m*Kvp*X(6)-v*h-q=0
-0.5*bp*t*h*h+bp*t*h*(h-X(1))-ap*h*X(1)+X(2)*X(1)+(Kvp/Ep)*X(1)*X(6)=0
0.5*bs*t*h*h-bs*t*h*(h-X(1))-as*h*X(1)+X(3)*X(1)+(Kvs/Es)*X(1)*X(7)=0
-(bp*t*h^3)/6+0.5*bp*t*h*h*(h-X(1))-0.5*ap*h*h*X(1)+X(2)*h*X(1)+X(4)*X(1)-X(1)*X(6)=0
(bs*t*h^3)/6-0.5*bs*t*h*h*(h-X(1))-0.5*as*h*h*X(1)+X(3)*h*X(1)+X(5)*X(1)-X(1)*X(7)=0
-(bp*t*(h-X(1))^3)/6+0.5*(bp*t*(h-X(1))-ap*X(1))*(h-X(1))^2+X(2)*(h-X(1))*X(1)+X(1)*X(4)-(bs*t*(h-X(1))^3)/6
+0.5*(bs*t*(h-X(1))+as*X(1))*(h-X(1))^2-X(3)*(h-X(1))*X(1)-X(1)*X(5)=0
Ep*(-0.5*bp*t*(h-X(1))^2+(bp*t*(h-X(1))-ap*X(1))*(h-X(1))+X(2)*X(1))
-Es*(0.5*bs*t*(h-X(1))^2-(bs*t*(h-X(1))+as*X(1))*(h-X(1))+X(3)*X(1))]=0


在Matlab用fsolve无法求解,在1stOpt里迭代求解的结果误差太大,急求各位高手指点!


附上1stOpt代码:
Constant h=5,d=4,R=0.7,U=R*2*pi,c=55000,fi=35,K0=0.95-sin(fi*pi/180),Es=18000000,Ep=20000000,v=19000,t=v*h*K0*tan(fi*pi/180)+c,Ae=d^2,
Ap=pi*R^2,Ass=Ae-Ap,q=2000000,Kvp=150000000,Kvs=5000000,ap=v/Ep,bp=U/(Ep*Ap),as1=v/Es,bs=U/(Es*Ass),m=Ap/Ae;
Parameter x1,x2,x3,x4,x5,x6,x7;
Function   (1-m)*Kvs*x7+m*Kvp*x6-v*h-q=0;
-0.5*bp*t*h*h+bp*t*h*(h-x1)-ap*h*x1+x2*x1+(Kvp/Ep)*x1*x6=0;
0.5*bs*t*h*h-bs*t*h*(h-x1)-as1*h*x1+x3*x1+(Kvs/Es)*x1*x7=0;
-(bp*t*(h^3))/6+0.5*bp*t*h*h*(h-x1)-0.5*ap*h*h*x1+x2*h*x1+x4*x1-x1*x6=0;
(bs*t*(h^3))/6-0.5*bs*t*h*h*(h-x1)-0.5*as1*h*h*x1+x3*h*x1+x5*x1-x1*x7=0;
-(bp*t*((h-x1)^3))/6+0.5*(bp*t*(h-x1)-ap*x1)*((h-x1)^2)+x2*(h-x1)*x1+x1*x4-(bs*t*((h-x1)^3))/6+0.5*(bs*t*(h-x1)+as1*x1)*((h-x1)^2)-x3*(h-x1)*x1-x1*x5=0;
Ep*(-0.5*bp*t*((h-x1)^2)+(bp*t*(h-x1)-ap*x1)*(h-x1)+x2*x1)-Es*(0.5*bs*t*((h-x1)^2)-(bs*t*(h-x1)+as1*x1)*(h-x1)+x3*x1)=0;

qibbxxt 发表于 2010-11-29 09:34

请你把你的matlab代码贴出来,大家帮你分析下

zzkingg 发表于 2010-11-29 16:05

回复 2 # qibbxxt 的帖子

function y=myfun(X)
h=5;
d=4;
R=0.7;
U=R*2*pi;
c=55000;
fi=35;
K0=0.95-sin(fi*pi/180);
Es=18000000;
Ep=20000000;
v=19000;
t=v*h*K0*tan(fi*pi/180)+c;
Ae=d^2;
Ap=pi*R^2;
As=Ae-Ap;
q=2000000;
Kvp=150000000;
Kvs=5000000;
ap=v/Ep;
bp=U/(Ep*Ap);
as=v/Es;
bs=U/(Es*As);
m=Ap/Ae;
y=[(1-m)*Kvs*X(7)+m*Kvp*X(6)-v*h-q;
    -0.5*bp*t*h*h+bp*t*h*(h-X(1))-ap*h*X(1)+X(2)*X(1)+(Kvp/Ep)*X(1)*X(6);
    0.5*bs*t*h*h-bs*t*h*(h-X(1))-as*h*X(1)+X(3)*X(1)+(Kvs/Es)*X(1)*X(7);
    -(bp*t*h^3)/6+0.5*bp*t*h*h*(h-X(1))-0.5*ap*h*h*X(1)+X(2)*h*X(1)+X(4)*X(1)-X(1)*X(6);
    (bs*t*h^3)/6-0.5*bs*t*h*h*(h-X(1))-0.5*as*h*h*X(1)+X(3)*h*X(1)+X(5)*X(1)-X(1)*X(7);
    -(bp*t*(h-X(1))^3)/6+0.5*(bp*t*(h-X(1))-ap*X(1))*(h-X(1))^2+X(2)*(h-X(1))*X(1)+X(1)*X(4)-(bs*t*(h-X(1))^3)/6+0.5*(bs*t*(h-X(1))+as*X(1))*(h-X(1))^2-X(3)*(h-X(1))*X(1)-X(1)*X(5);
    Ep*(-0.5*bp*t*(h-X(1))^2+(bp*t*(h-X(1))-ap*X(1))*(h-X(1))+X(2)*X(1))-Es*(0.5*bs*t*(h-X(1))^2-(bs*t*(h-X(1))+as*X(1))*(h-X(1))+X(3)*X(1))];X0=;
options=optimset('Display','iter');
= fsolve(@myfun,X0,options)

qibbxxt 发表于 2010-11-29 16:22

x =

    0.9917
   -0.2944
   -0.2132
    0.1889
    0.1255
    0.1227
    0.0718


fval =

   -0.0000
    0.7020
   -0.2071
   -1.0711
   -1.0566
    0.0155
    0.0171
这是我运行的结果

zzkingg 发表于 2010-11-29 17:24

回复 4 # qibbxxt 的帖子

对啊,我也是这个,但是Matlab说是

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunEvals = 700 (the default value).

zhouyang664 发表于 2010-11-29 17:57

把MaxFunEvals设置大一点就可以了!

qibbxxt 发表于 2010-11-29 21:02

回复 5 # zzkingg 的帖子

这也许和设置有关系,你可以仔细看看帮助文档,调整一下参数,再试一试

lvhaiwei007 发表于 2010-12-2 21:21

回复 3 # zzkingg 的帖子

不知道下面问题的Matlab程序会不会?
方程组中有五个方程,六个未知量,要求画x与y的关系。
(0.8-0.01*x.^2).*a11^4-0.02.*x.*y.*a11^3-(x.^2+0.01.*y^2-5).*a11^2-2.*x.*y.*a11-y.^2=0
(0.8-0.01*x.^2).*a21^4-0.02.*x.*y.*a21^3-(x.^2+0.01.*y^2-5).*a21^2-2.*x.*y.*a21-y.^2=0
(0.8-0.01*x.^2).*a31^4-0.02.*x.*y.*a31^3-(x.^2+0.01.*y^2-5).*a31^2-2.*x.*y.*a31-y.^2=0
(0.8-0.01*x.^2).*a41^4-0.02.*x.*y.*a41^3-(x.^2+0.01.*y^2-5).*a41^2-2.*x.*y.*a41-y.^2=0
(a11.^2-a21.^2)(a31.^2-a41.^2)+(a21.^2-a41.^2)(a31.^2-a11.^2)+(a11.^2-a41.^2)(a21.^2-a31.^2)=0

huaijuliu 发表于 2010-12-9 08:03

如果matlab提供的fsolve不够合适,可以参考一些网上可以搜索到的迭代法程序进行方程组的求解,如低松弛雅克比法或者高斯赛戴尔法求解。

对于楼上的问题,最起码有一个笨方法,在程序中添加x=c, c is a constant。然后解出y。然后再换一个x=c2,解出y。最后通过数据拟合得到公式。
页: [1]
查看完整版本: 急求此非线性方程组的解法?谢谢!