问题背景

求解方程组的时候，对某些未知数的求解结果的取值范围有要求。例如在某些物理问题求解中，要求待求解量大于0。

代码

一共两个文件:

my_fun.m
main.m

my_fun.m

function F=my_fun(x)
R = 10;
yinzi = pi/180;
alpha_ON = 50*yinzi;
alpha_OM = 10*yinzi; 
alpha_MN = 40*yinzi;
F(1)=sin(alpha_ON)*sin(x(2))-sin(alpha_OM)*sin(x(3));
F(2)=2*sin((x(1)-1)*pi/9)*sin(alpha_ON)*cos((x(1)-1)*pi/9+x(2)) - sin(alpha_MN)*sin(x(3)+alpha_ON);
F(3)=x(2)+ 2*(x(1)-1)*pi/9- x(3) - alpha_MN;
difference = R/sin(alpha_OM)*sin(x(2)) - R;
tolerance = 1e-5;
if abs(difference) < toleranceF(4) = 0;
elseF(4) = difference;
end
end

核心是这里:

difference = R/sin(alpha_OM)*sin(x(2)) - R;
tolerance = 1e-5;
if abs(difference) < toleranceF(4) = 0;
elseF(4) = difference;

设置difference和tolerance，定义F(4),在求解过程中，会让F(4)等于0，也就是令
$$\frac{R}{sin(alph{a}_{OM})}\times sin(x(2))=R$$
这个条件逐渐满足。

main.m

%x0=[0,pi/18,5*pi/18];
x0 = zeros(1,3)
options.Display = 'iter';
R = 10;
result_x=fsolve(@my_fun,x0,options)
yinzi = pi/180;
alpha_ON = 50*yinzi;
alpha_OM = 10*yinzi; 
alpha_MN = 40*yinzi;
pho_x = R/sin(alpha_OM)*sin(result_x(2))

结果展示:

不加入F(4)

result_x =1.6379   -0.0725   -0.3253pho_x =-4.1722

加入F(4)

result_x =2.9988    0.1745    0.8718pho_x =10.0000

这个解正确。