求解非线性方程组遇到的问题

zijiu

这是一个比较复杂的四元非线性方程组

m文件代码如下：

function [x1,x2,x3,x4]=mf8(X)

rou=7;
kai=8.2;
niu=3;
you=0.1;
yeta=0.1;
n=0.2;
sigma=10;
theta=1.9;
a=0.3;
E=1;
Cj=1;
Ck=1;
Kj=1;
Kk=1;



Ij=you*Kj;
Ik=you*Kk;
MCh=((X(1)/a)^a)*(((1-a)/X(3))^(a-1));
MCf=((X(2)/a)^a)*(((1-a)/X(4))^(a-1));
ph=sigma*MCh/(sigma-1);
phx=sigma*MCh/((sigma-1)*E);
pfx=sigma*MCf/(sigma-1);
pf=E*pfx;
PH=ph;
PF=pf;
PHx=phx;
PFx=pfx;
P=(yeta*PH^(1-theta)+(1-yeta)*PF^(1-theta))^(1/(1-theta));
Px=(yeta*PHx^(1-theta)+(1-yeta)*PFx^(1-theta))^(1/(1-theta));
STj=niu*P*Cj^rou/X(3);
STk=niu*Px*Ck^rou/X(4);
Yh=(Kj^a)*(STj^(1-a));
Yf=(Kk^a)*(STk^(1-a));
yh=yeta*(ph/PH)^(-sigma)*(PH/P)^(-theta)*(n*Cj+n*Ij)/n;
yf=(1-yeta)*(pf/PF)^(-sigma)*(PF/P)^(-theta)*(n*Cj+n*Ij)/(1-n);
yhx=yeta*(phx/PHx)^(-sigma)*(PHx/Px)^(-theta)*((1-n)*Ck+(1-n)*Ik)/n;
yfx=(1-yeta)*(pfx/PFx)^(-sigma)*(PFx/Px)^(-theta)*((1-n)*Ck+(1-n)*Ik)/(1-n);
x1=Yh-yh-yhx;
x2=Yf-yf-yfx;
x3=yh*ph+yhx*phx-P*Cj-P*Ij;
x4=yf*pf+yfx*pfx-Px*Ck-Px*Ik;
end







%上述函数是一组非线性方程，X为自变量，包含4个元素：X(1),X(2),X(3),X(4)。我想用寻找到这样的X，使得最后四个方程中x1,x2,x3,x4等于0

%求解的代码：
options=optimset('Display','iter','MaxIter',200000, 'NonlEqnAlgorithm', 'gn','MaxFunEvals',10000000);
[X,FVAL,EXITFLAG,OUTPUT]=fsolve(@mf8,[4,3,2,9],options);
%反复测试后用了上面这一组初值

返回的结果：


                                                     Directional
Iteration  Func-count    Residual     Step-size      derivative
     0           5         22.1216
     1          15         3.53439         4.47           -1.77
Iteration matrix ill-conditioned - Switching to LM method.
     2          23          1.7281        0.235           -11.5       0.100634
     3          31       0.0453446        0.806            -1.3        2.93274
     4          39    2.82673e-006        0.906        -0.00088        1.44781
     5          47    1.60298e-019         1.01      -1.34e-012       0.757026
Optimization terminated: directional derivative along
search direction less than TolFun and infinity-norm of
gradient less than 10*(TolFun+TolX).

X =

    0.1304    2.9124    8.4168    8.9318


FVAL =

  4.0037e-010


EXITFLAG =

     1


OUTPUT =

       iterations: 6
        funcCount: 47
         stepsize: 1
     cgiterations: []
    firstorderopt: []
        algorithm: 'medium-scale: Levenberg-Marquardt, line-search'
          message: [1x147 char]



%根据返回结果看来，应该找到解了，因为EXITFLAG =     1

%但是把求得的解带入函数：
[x1,x2,x3,x4]=mf8([0.1304,2.9124,8.4168,8.9318]);

%得到的结果是：


x1 =

-3.8486e-004


x2 =

    1.6040


x3 =

    0.6598


x4 =

   -0.1650

%x1,x2,x3,x4都不等于0。为什么会出现这样矛盾的情况呢？另外，如果把初值略作改变，比如用[3,2,2,9]，也可以得到另一组解：[0.1596    2.0330    5.4434    9.0171]，同样地，带入函数，并不能使函数值为0


求高手帮助，谢谢！

zijiu

为了看看我对fsolve的理解有没有问题，我用了下面这个简单的方程组做个试验：

m函数：


function [y1,y2]=xxx(X)
y1=X(1)-X(2);
y2=X(1)^2-X(2);


用下面的代码求解
[X,FVAL,EXITFLAG,OUTPUT]=fsolve(@xxx,[0.1,0.1],options)

结果返回了


>> [X,FVAL,EXITFLAG,OUTPUT]=fsolve(@xxx,[0.1,0.1],options)

                                                     Directional
Iteration  Func-count    Residual     Step-size      derivative
     0           3               0
Optimization terminated: directional derivative along
search direction less than TolFun and infinity-norm of
gradient less than 10*(TolFun+TolX).

X =

    0.1000    0.1000


FVAL =

     0


EXITFLAG =

     1


OUTPUT =

       iterations: 1
        funcCount: 3
         stepsize: 1
     cgiterations: []
    firstorderopt: []
        algorithm: 'medium-scale: Gauss-Newton, line-search'
          message: [1x147 char]


很显然，两个正确的解应该是[0 0],[1 1]。这么简单的方程，都出错了，fsolve应该不至于这么弱智吧。是我对fsolve的理解和使用有问题吗？

chenguijia

建议你重新装一下matlab,我也用matlab解方程，之所以不为零，是因为迭代法判定结果收敛到所要求的精度，就不再进行迭代了，所以结果不为零，第二个问题是迭代法最还给定的初值要接近真实值这样的结果才具有可靠性。
个人建议，仅供参考

dingd

试下1stOpt，不需要初值。

zijiu

谢谢楼上两位！

我来找个matlab新版本试试看，是不是能解决。

1stopt我之前也用过，解相对比较简单的方程组没有问题，比matlab厉害，像我这样的复杂一点的方程组好像就不行了，每一次都是不同的解，代入结果测试，也不为0。

而且我故意用一个无解的方程组测试，1stopt还能像模像样地算出“结果”来，每次也给出不同的解。

如果matlab新版本也解决不了这个问题，那我该怎么办呢。难道说可以认定这个方程组无解？

zijiu

自己再顶顶

ChaChing

回复 2 # zijiu 的帖子

当然fsolve不至于这么弱智, 一定是使用错误造成!:@)
首先, 需承认个人时间有限&懒惰, 1F的内容太长不太想细看!
就仅试试2F的东东, 我想问题原因应该一个样
同样的, 需请LZ再看下fsolve的help, 我也是仅能如此

1. function ff=xxx(X)
   
2. ff=[X(1)-X(2);
   
3. X(1)^2-X(2)];

复制代码

1. [X,FVAL,EXITFLAG,OUTPUT]=fsolve(@xxx,[0.1,0.1],options)
   
2. [X,FVAL,EXITFLAG,OUTPUT]=fsolve(@xxx,[0.8,0.9],options)

复制代码

ChaChing

回复 2 # zijiu 的帖子

LZ知道为何错误了吗?
请将xxx.m改为下式再试试

1. function y1=xxx(X)
   
2. y1=X(1)-X(2);
   
3. y2=X(1)^2-X(2);

复制代码

是不是与LZ原来情况相同?
原因是仅y1被求解!
请LZ细看下fsolve的help, 裡头有例子:@)

小六儿

看看问题