在Matlab中解非线性方程组遇到的问题！

zgq1004118

刚刚接触MATLAB不久，属于 菜菜鸟 级的。在求解一组非线性方程时无法收敛！是否有其他解法？请高人指教！！
function non_linearfun
x0=[0 0 0];
answer_x_y_z=fsolve(@fun,x0)
function y=fun(x)
%----------------------Variable-------------------
%x=sigma_r./A;
%y=sigma_a./A
%z=Theta/A  
%

Kn      =   434912.7857;
A       =   0.9922;
B       =   0.05;
D       =   19.844;
Ri      =   40.4027718;
dm      =   80;
n       =   1.5;
alpha0  =   40/180*pi;
%
%---------------------------------------------------
sinao   =   sin(alpha0);
cosao   =   cos(alpha0);

phi0    =   0;
phi1    =   2/11*pi;
phi2    =   2*2/11*pi;
cosphi0 =   cos(phi0);
cosphi1 =   cos(phi1);
cosphi2 =   cos(phi2);

%------------------------------------------------
B0=sinao+x(1)+Ri*x(3)*cosphi0;
B1=sinao+x(1)+Ri*x(3)*cosphi1;
B2=sinao+x(1)+Ri*x(3)*cosphi2;
%------------------------------------
%------------------------------------
C0=cosao+x(2)*cosphi0;
C1=cosao+x(2)*cosphi1;
C2=cosao+x(2)*cosphi2;
%----------------------------------
%----------------------------------
D0=(cosao+x(2)*cosphi0)*cosphi0;
D1=(cosao+x(2)*cosphi1)*cosphi1;
D2=(cosao+x(2)*cosphi2)*cosphi2;

%-----------------------------------
Fa      =   10000;
Fr      =   0;
My      =   0;                                    
%---------------------------------------

y(1)=   Fa-Kn*(A^n)*(((B0^2+C0^2)^.5-1)^1.5*B0/(B0^2+C0^2)^0.5+2*((B1^2+C1^2)^0.5-1)^n*B1/(B1^2+C1^2)^0.5+2*((B2^2+C2^2)^0.5-1)^n*B2/(B2^2+C2^2)^.5);
y(2)=   Fr-Kn*(A^n)*(((B0^2+C0^2)^.5-1)^1.5*D0/(B0^2+C0^2)^0.5+2*((B1^2+C1^2)^0.5-1)^n*D1/(B1^2+C1^2)^0.5+2*((B2^2+C2^2)^0.5-1)^n*D2/(B2^2+C2^2)^.5);
y(3)=   My-0.5*dm*Kn*(A^n)*(((B0^2+C0^2)^.5-1)^1.5*B0*cosphi0/(B0^2+C0^2)^0.5+2*((B1^2+C1^2)^0.5-1)^n*B1*cosphi1/(B1^2+C1^2)^0.5+2*((B2^2+C2^2)^0.5-1)^n*B2*cosphi2/(B2^2+C2^2)^.5);

zgq1004118

我再解释一下:@) ，三个未知数x(1),x(2),x(3).
B0-B2,C0-C2,D0-D2是含有未知数的变量。
其它都是常量，其数值已经给出。这是一个很复杂的方程组，已被我简化写出。
高人啊，快来帮帮兄弟啊:@L

xjzuo

请把你的＂无法收敛＂贴出来看看．
另：我运行了一下，似乎可以得到结果．

zgq1004118

> Warning: Cannot determine from calling sequence whether to use new
  (2.0 or later) FSOLVE function or grandfathered FSOLVE function.
  Assuming new syntax; if call was grandfathered FSOLVE syntax,
  this may give unexpected results.
  To force new syntax and avoid warning, add options structure argument:
      x = fsolve(@sin,3,optimset('fsolve'));
  To force grandfathered syntax and avoid warning, add options array argument:
      x = fsolve(@sin,3,foptions);

> In E:\Matlab\toolbox\optim\fsolve.m (parse_call) at line 389
  In E:\Matlab\toolbox\optim\fsolve.m at line 101
  In E:\Matlab\work\non_linearfunz6.m at line 3
Optimization terminated successfully:
Norm of relative change in X is less than max(options.TolX^2,eps) and
sum-of-squares of function values is less than sqrt(options.TolFun)

answer_x_y_z =

  Column 1

          0.261630241614871 +    0.0432136503412044i

  Column 2

         -0.147501542528093 -    0.0461512387637423i

  Column 3

       -0.00424579357031622 -   0.00155455859937376i
这是我运行的结果，是完全数解，得不到实数解，请问高人的解怎样？

xjzuo

由于你的问题背景没有讲清楚,所以不知道你要求什么样的解.
我试了试以下指令,就你给的这段简化程序而言,似乎并非你说的不收敛:
%%%-------------------------------------------------------%%%
x0=[0 0 0];
x = fsolve(@myfunt1,x0,optimset('fsolve'))
%%%-------------------------------------------------------%%%
x0=[1 1 1];
x = fsolve(@myfunt1,x0,optimset('fsolve'))
%%%-------------------------------------------------------%%%

zgq1004118

多谢版主关心。我现在把方程换个写法。三个未知数:x(1),x(2),x(3);三个方程是y(1),y(2),y(3);
欲求出一组实数解！！这样写出来便于理解方程组的结构，就是方程本身太复杂，不便于在matlab中编程，因此就没有写出原始表达式。
不知版主上次选取x0=[0 0 0]做初值，在运行时的结果是否跟我的一样也是完全数解？

function non_linearfun
x0=[0 0 0];
answer_x_y_z=fsolve(@fun,x0)
function y=fun(x)
%----------------------function-------------------
y(1)=10000-429834.2416*((((.642787609+x(1)+40.304027718*x(3))^2+(x(2)+.766044443)^2)^.5-1)^1.5*(.642787609+x(1)+40.304027718*x(3))/((.642787609+x(1)+40.304027718*x(3))^2+(x(2)+.766044443)^2)^.5+2*(((.642787609+x(1)+33.9059057*x(3))^2+(.841253532*x(2)+.766044443)^2)^.5-1)^1.5*(.642787609+x(1)+33.9059057*x(3))/((.642787609+x(1)+33.9059057*x(3))^2+(.841253532*x(2)+.766044443)^2)^.5+2*(((.642787609+x(1)+16.7428982*x(3))^2+(.415415013*x(2)+.766044443)^2)^.5-1)^1.5*(.642787609+x(1)+16.7428982*x(3))/((.642787609+x(1)+16.7428982*x(3))^2+(.415415013*x(2)+.766044443)^2)^.5);
y(2)=0-429834.2416*((((.642787609+x(1)+40.304027718*x(3))^2+(x(2)+.766044443)^2)^.5-1)^1.5*(.766044443+x(2))/((.642787609+x(1)+40.304027718*x(3))^2+(x(2)+.766044443)^2)^.5+2*(((.642787609+x(1)+33.9059057*x(3))^2+(.841253532*x(2)+.766044443)^2)^.5-1)^1.5*(.766044443+.841253532*x(2))*.841253532/((.642787609+x(1)+33.9059057*x(3))^2+(.841253532*x(2)+.766044443)^2)^.5+2*(((.642787609+x(1)+16.7428982*x(3))^2+(.415415013*x(2)+.766044443)^2)^.5-1)^1.5*(.766044443+.415415013*x(2))*.415415013/((.642787609+x(1)+16.7428982*x(3))^2+(.415415013*x(2)+.766044443)^2)^.5);
y(3)=0-40*429834.2416*((((.642787609+x(1)+40.304027718*x(3))^2+(x(2)+.766044443)^2)^.5-1)^1.5*(.642787609+x(1)+40.304027718*x(3))/((.642787609+x(1)+40.304027718*x(3))^2+(x(2)+.766044443)^2)^.5+2*(((.642787609+x(1)+33.9059057*x(3))^2+(.841253532*x(2)+.766044443)^2)^.5-1)^1.5*(.642787609+x(1)+33.9059057*x(3))*.841253532/((.642787609+x(1)+33.9059057*x(3))^2+(.841253532*x(2)+.766044443)^2)^.5+2*(((.642787609+x(1)+16.7428982*x(3))^2+(.415415013*x(2)+.766044443)^2)^.5-1)^1.5*(.642787609+x(1)+16.7428982*x(3))*.415415013/((.642787609+x(1)+16.7428982*x(3))^2+(.415415013*x(2)+.766044443)^2)^.5);

xjzuo

我认为还是一楼的表达式比较清楚．
就你给的表达式，似乎没有实数解．
另：你说的＂完全数＂，应该不是指数论中的完全数吧？

[ 本帖最后由 xjzuo 于 2007-3-1 08:59 编辑 ]

dingd

第一个公式中：

"10000-429834.2416*(...)"变为“1000000-429834.2416*()"或”100000-42983.42416*()'就会有实数解了，否则，似乎无实数解。确认代码没输错吗？

zgq1004118

8楼的高手我前不久在本论坛读过你的关于非线性方程组解法的帖子，从中收获很多，在这谢了。
下面是把"10000-429834.2416*(...)"变为“1000000-429834.2416*()"后的结果，
仍然不是实数解。不知高手的结如何，能否贴出来？另外是否有其它解法？:@)
answer_x_y_z =

  Column 1

           3.64616251216352 -     0.644111407183738i

  Column 2

         -0.742140780488419 +     0.147938243941897i

  Column 3

         -0.100008046487466 +    0.0360334493768073i

dingd

"10000-429834.2416*(...)"变为“1000000-429834.2416*(...)"，结果：

x1: 5.44662954384634
x3: -0.210122988803927
x2: -1.06537137387729

dingd

再精确点：

x1: 5.44663122246003
x3: -0.210123048907879
x2: -1.06536950937869

zgq1004118

多谢这位高人，但10楼和11楼的解也不在取值范围内，另外我用MATLAB是得不到这个解的。郁闷啊:@( :@o

grta

用solve函数试试那