如何在 MATLAB 中获得该方程的所有解？答案

【问题标题】：How can I get all solutions to this equation in MATLAB?如何在 MATLAB 中获得该方程的所有解？
【发布时间】：2016-11-21 16:25:08
【问题描述】：

我想求解以下方程：tan(x) = 1/x

我做了什么：

syms x
eq = tan(x) == 1/x;
sol = solve(eq,x)

但这给了我解决方案的一个数值近似值。之后，我阅读了以下内容：

[sol, params, conds] = solve(eq, x, 'ReturnConditions', true)

但这告诉我它找不到明确的解决方案。

如何在给定范围内找到这个方程的数值解？

【问题讨论】：

“获取所有值”是什么意思？这个方程有无限的解（只是plot it and see）。
对不起，我的意思不止一个值。例如在一个范围内。
我的意思是，一定范围内的数值解

标签： matlab numeric equation numerical-methods equation-solving

【解决方案1】：

我从不喜欢“盲目地”使用求解器，也就是说，没有某种体面的初始值选择方案。根据我的经验，你在盲目做事时会发现的价值观也没有上下文。意思是，你经常会错过解决方案，认为某些东西是解决方案，而实际上解决方案爆炸了，等等。

对于这种特殊情况，重要的是要意识到fzero 使用数值导数来找到越来越好的近似值。但是，随着x 的增加，f(x) = x · tan(x) - 1 的导数越来越难以准确计算：

如您所见，x 越大，f(x) 越接近垂直线； fzero 会爆炸！因此，在输入fzero 之前，必须尽可能接近得到一个估计值。

所以，这是一种获得良好初始值的方法。

考虑函数

f(x) = x · tan(x) - 1

知道tan(x) 有Taylor expansion：

tan(x) ≈ x + (1/3)·x³ + (2/15)·x⁵ + (7/315)·x⁷ + ...

我们可以用它来近似函数f(x)。在第二项之后截断，我们可以写成：

f(x) ≈ x · (x + (1/3)·x³) - 1

现在，要实现的关键是 tan(x) 与句点 π 重复。因此，考虑函数族是最有用的：

fₙ(x) ≈ x · ( (x - n·π) + (1/3)·(x - n·π)³) - 1

评估几个倍数并收集术语得出以下概括：

f₀(x) = x⁴/3 - 0π·x³ + ( 0π² + 1)x² - (0π +   (0π³)/3)·x  - 1
f₁(x) = x⁴/3 - 1π·x³ + ( 1π² + 1)x² - (1π +   (1π³)/3)·x  - 1
f₂(x) = x⁴/3 - 2π·x³ + ( 4π² + 1)x² - (2π +   (8π³)/3)·x  - 1
f₃(x) = x⁴/3 - 3π·x³ + ( 9π² + 1)x² - (3π +  (27π³)/3)·x  - 1
f₄(x) = x⁴/3 - 4π·x³ + (16π² + 1)x² - (4π +  (64π³)/3)·x  - 1
                              ⋮
fₙ(x) = x⁴/3 - nπ·x³ + (n²π² + 1)x² - (nπ +  (n³π³)/3)·x  - 1

在一个简单的 MATLAB 测试中实现所有这些：

% Replace this with the whole number of pi's you want to 
% use as offset 
n = 5;

% The coefficients of the approximating polynomial for this offset
C = @(npi) [1/3
            -npi
            npi^2 + 1
            -npi - npi^3/3
            -1];

% Find the real, positive polynomial roots
R = roots(C(n*pi));
R = R(imag(R)==0);
R = R(R > 0);

% And use these as initial values for fzero()
x_npi = fzero(@(x) x.*tan(x) - 1, R)

在一个循环中，这可以产生下表：

% Estimate (polynomial)    Solution (fzero)
 0.889543617524132          0.860333589019380    0·π
 3.425836967935954          3.425618459481728    1·π
 6.437309348195653          6.437298179171947    2·π
 9.529336042900365          9.529334405361963    3·π
12.645287627956868         12.645287223856643
15.771285009691695         15.771284874815882    
18.902410011613000         18.902409956860023
22.036496753426441         22.036496727938566    ⋮
25.172446339768143         25.172446326646664    
28.309642861751708         28.309642854452012
31.447714641852869         31.447714637546234
34.586424217960058         34.586424215288922   11·π

如您所见，近似值基本上等于解。对应情节：

【讨论】：

为什么不简单地乘以违规项 cos(x) 来找到无极方程 x*sin(x)-cos(x)==0？
@LutzL 嗯....然后呢？仍然很难找到一个好的初步估计。这是一个与角度相关的相量，意思是，它与√(x²+1) · sin(x + atan2(-1,x)) 相同。如您所见，周期不是“干净”的 2π，而是取决于 x。您可以将n·π 添加到您找到的第一个解决方案中以获得后续估计，对于这种情况，我认为这确实足够好。但是，经验表明，对于周期与您在初始估计中假设的完全不同的其他种情况，您最终将两次收敛到相同的解决方案，并错过下一个解决方案一个……
@LutzL 好吧，你可以随意添加一个让我看起来很可笑的答案，当这种情况发生时我真的很喜欢 :)
做到了，把我的观点变成答案。 -- 您是否正在考虑修改问题tan(a*x)=1/x 或其他附加条款？

【解决方案2】：

要在某个范围内找到函数的数值解，可以像这样使用fzero：

fun = @(x)x*tan(x)-1; % Multiplied by x so fzero has no issue evaluating it at x=0.
range = [0 pi/2];
sol = fzero(fun,range);

以上将只返回一个解决方案 (0.8603)。如果您需要其他解决方案，则必须多次致电fzero。例如，这可以在循环中完成：

fun = @(x)tan(x)-1/x;

RANGE_START = 0;
RANGE_END = 3*pi;
RANGE_STEP = pi/2;

intervals = repelem(RANGE_START:RANGE_STEP:RANGE_END,2);
intervals = reshape(intervals(2:end-1),2,[]).';

sol = NaN(size(intervals,1),1);

for ind1 = 1:numel(sol)
  sol(ind1) = fzero(fun, mean(intervals(ind1,:)));
end

sol = sol(~isnan(sol)); % In case you specified more intervals than solutions.

这给出了：

[0.86033358901938;
 1.57079632679490; % Wrong
 3.42561845948173;
 4.71238898038469; % Wrong
 6.43729817917195; 
 7.85398163397449] % Wrong

注意：

函数是对称的，它的根也是对称的。这意味着您可以仅求解正区间（例如）并“免费”获得负根。
sol 中的所有其他条目都是错误的，因为这是我们有渐近不连续性的地方（tan 从+Inf 转换到-Inf），这被 MATLAB 错误地识别为解决方案。所以你可以忽略它们（即sol = sol(1:2:end);。

【讨论】：

谢谢，但它只显示一个解决方案。
将等式更改为 x*sin(x)-cos(x) == 0 以进一步去奇异化。

【解决方案3】：

将等式乘以x 和cos(x) 以避免任何可能具有值0 的分母，

f(x)=x*sin(x)-cos(x)==0

考虑归一化函数

h(x)=(x*sin(x)-cos(x)) / (abs(x)+1)

对于大x，这将越来越接近sin(x)（或-sin(x)，对于大负x）。实际上，对于x>pi，绘制此图在视觉上已经是真实的，直到幅度因子。

对于[0,pi/2] 中的第一个根，使用二阶x^2-(1-0.5x^2)==0 处x=0 处的泰勒近似得到x[0]=sqrt(2.0/3) 作为根近似，对于更高的根，取正弦根x[n]=n*pi，n=1,2,3,... 作为牛顿迭代中的初始近似值xnext = x - f(x)/f'(x) 得到

 n        initial              1. Newton         limit of Newton

 0    0.816496580927726    0.863034004302817    0.860333589019380
 1    3.141592653589793    3.336084918413964    3.425618459480901
 2    6.283185307179586    6.403911810682199    6.437298179171945
 3    9.424777960769379    9.512307014150883    9.529334405361963
 4   12.566370614359172   12.635021895208379   12.645287223856643
 5   15.707963267948966   15.764435036320542   15.771284874815882
 6   18.849555921538759   18.897518573777646   18.902409956860023
 7   21.991148575128552   22.032830614521892   22.036496727938566
 8   25.132741228718345   25.169597069842926   25.172446326646664
 9   28.274333882308138   28.307365162331923   28.309642854452012
10   31.415926535897931   31.445852385744583   31.447714637546234
11   34.557519189487721   34.584873343220551   34.586424215288922

【讨论】：