马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
这是一个比较复杂的四元非线性方程组
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
求高手帮助,谢谢! |