43..利用fsolve函数解对应lambda下的方程组（matlab程序）

1.简述

fsolve的基本用法 : x = fsolve(fun,x0)

其中fun应为函数句柄，x0为搜索的种子，即预估的fun=0的解的大致位置。

函数句柄的定义方式主要有两种：

1.定义函数文件，使用@操作符

定义function文件root2d.m, 如下：

function F = root2d(x)
    F(1) = exp(-exp(-(x(1)+x(2)))) - x(2)*(1+x(1)^2);
    F(2) = x(1)*cos(x(2)) + x(2)*sin(x(1)) - 0.5;
end

命令行运行：

fun = @root2d;      % 获得函数句柄 
x0 = [0,0];         % 定义种子 
x = fsolve(fun,x0)  % 求解

或者

x0 = [0,0];         % 定义种子 
x = fsolve(@root2d,x0)  % 求解

即可获得目标解。

这时会有提示：

fsolve completed because the vector of function values is near zero

意思是：fsolve完成，因为函数值的向量接近零。这代表求解成功，其他情况的提示一般为求解失败，各种报错陆续更新。

fsolve的进阶用法

上面的定义方式可以很简单的获得非线性方程的根，但是需要定义function文件，使用起来不是很方便，尤其是当需要求解 �(�)=�0 时，反复修改function是非常不划算的。这时候，采用匿名函数进行求解就非常方便了。

以求解 �∗�∗�=[1,2;3,4] 为例，运行

fun = @(x)x*x*x - [1,2;3,4]; 
x0 = ones(2);
fsolve(fun,x0)

即可获得想要的结果。

以求解上一节中的问题为例

fun = @(x)[exp(-exp(-(x(1)+x(2)))) - x(2)*(1+x(1)^2);...
x(1)*cos(x(2)) + x(2)*sin(x(1)) - 0.5;] - [0.1 0.1];
x0 = [0,0];  
fsolve(fun,x0)

2.代码

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)]) %给出图的标题

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