MATLAB求解非线性方程组fsolve源程序代码

function equation()
global sigma mu T lambda
sigma=5;   %定义sigma的值
mu=0.4;    %定义mu的值
T=1.7;     %定义T的值
N=1;
for lambda=0:0.05:1   %lambda从0到1，隔0.05计算一个点
    x0=[0.5000 0.5000 1.1817 1.1817 1.0000 1.0000 0.9354 0.9354]*1;    %定义迭代初值
    x=fsolve(@myfun,x0);    %采用fsolve函数解对应lambda下的方程组，结果保存在x里
    value(N)=x(7)-x(8);     %求出对应lambda下的omega1-omega2的值，保存在value里
    N=N+1;
end
lambda=0:0.05:1;
plot(lambda,value)   %绘图
title(['T=',num2str(T)])   %给出图的标题
msgbox('MATLAB编程答疑，请加QQ: 1530497909','MATLAB答疑','help')
web http://url.cn/TKcdXk -browser
end

function F=myfun(x)
global sigma mu T lambda
%x(1)~x(8)分别对应8个未知数
Y1=x(1);    
Y2=x(2);
G1=x(3);
G2=x(4);
w1=x(5);
w2=x(6);
omega1=x(7);
omega2=x(8);
%定义8个方程
eq1=Y1-mu*lambda*w1-(1-mu)/2;
eq2=Y2-mu*(1-lambda)*w2-(1-mu)/2;
eq3=G1-(lambda*w1^(1-sigma)+(1-lambda)*(w2*T)^(1-sigma))^(1/(1-sigma));
eq4=G2-(lambda*(w1*T)^(1-sigma)+(1-lambda)*w2^(1-sigma))^(1/(1-sigma));
eq5=w1-(Y1*G1^(sigma-1)+Y2*(G2^(sigma-1))*T^(1-sigma))^(1/sigma);
eq6=w2-(Y1*(G1^(sigma-1))*T^(1-sigma)+Y2*G2^(sigma-1))^(1/sigma);
eq7=omega1-w1*G1^(-mu);
eq8=omega2-w2*G2^(-mu);
%返回方程组
F=[eq1;eq2;eq3;eq4;eq5;eq6;eq7;eq8];
end