问题背景

求解方程组的时候，对某些未知数的求解结果的取值范围有要求。例如在某些物理问题求解中，要求待求解量大于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) < tolerance
    F(4) = 0;
else
    F(4) = difference;
end
end

核心是这里:

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

设置difference和tolerance，定义F(4),在求解过程中，会让F(4)等于0，也就是令
R s i n ( a l p h a O M ) × s i n ( x ( 2 ) ) = R \frac{R}{sin(alpha_{OM})} \times sin(x(2)) = R sin(alphaOM)R×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.3253


pho_x =

   -4.1722

加入F(4)

result_x =

    2.9988    0.1745    0.8718


pho_x =

   10.0000

这个解正确。

yscript

matlab求解方程组-求解过程中限制解的取值范围