急求此非线性方程组的解法？谢谢！

zzkingg

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

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

zzkingg

回复 2 # qibbxxt 的帖子

1. function y=myfun(X)
   
2. h=5;
   
3. d=4;
   
4. R=0.7;
   
5. U=R*2*pi;
   
6. c=55000;
   
7. fi=35;
   
8. K0=0.95-sin(fi*pi/180);
   
9. Es=18000000;
   
10. Ep=20000000;
    
11. v=19000;
    
12. t=v*h*K0*tan(fi*pi/180)+c;
    
13. Ae=d^2;
    
14. Ap=pi*R^2;
    
15. As=Ae-Ap;
    
16. q=2000000;
    
17. Kvp=150000000;
    
18. Kvs=5000000;
    
19. ap=v/Ep;
    
20. bp=U/(Ep*Ap);
    
21. as=v/Es;
    
22. bs=U/(Es*As);
    
23. m=Ap/Ae;
    
24. y=[(1-m)*Kvs*X(7)+m*Kvp*X(6)-v*h-q;
    
25.     -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);
    
26.     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);
    
27.     -(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);
    
28.     (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);
    
29.     -(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);
    
30.     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))];

复制代码

1. X0=[3;-0.2;-0.2;0.1;0.1;0.1;0.05];
   
2. options=optimset('Display','iter');
   
3. [x,fval] = fsolve(@myfun,X0,options)

复制代码

qibbxxt

1. x =
   
2. 
   
3.     0.9917
   
4.    -0.2944
   
5.    -0.2132
   
6.     0.1889
   
7.     0.1255
   
8.     0.1227
   
9.     0.0718
   
10. 
    
11. 
    
12. fval =
    
13. 
    
14.    -0.0000
    
15.     0.7020
    
16.    -0.2071
    
17.    -1.0711
    
18.    -1.0566
    
19.     0.0155
    
20.     0.0171
    

复制代码

这是我运行的结果

zzkingg

回复 4 # qibbxxt 的帖子

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

Solver stopped prematurely.

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

zhouyang664

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

qibbxxt

回复 5 # zzkingg 的帖子

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

lvhaiwei007

回复 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)[exp(i.*(a11+a21))+exp(i.*(a31+a41))]+(a21.^2-a41.^2)(a31.^2-a11.^2)[exp(i.*(a21+a41))+exp(i.*(a31+a11))]+(a11.^2-a41.^2)(a21.^2-a31.^2)[exp(i.*(a11+a41))+exp(i.*(a21+a31))]=0

huaijuliu

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

对于楼上的问题，最起码有一个笨方法，在程序中添加x=c, c is a constant。然后解出y。然后再换一个x=c2，解出y。最后通过数据拟合得到公式。