我从来不喜欢“盲目地”使用求解器,也就是说,没有某种合适的初始值选择方案 . 根据我的经验,盲目做事时你会发现的 Value 观也是没有背景的 . 意思是,你经常会错过解决方案,认为某些东西是解决方案,而实际上解算器爆炸了,等等 .
对于这种特殊情况,重要的是要认识到 fzero 使用数值导数来找到越来越好的近似值 . 但是, f(x) = x · tan(x) - 1 的衍生物变得越来越难以准确计算以增加 x :
如您所见, x 越大, f(x) 越接近垂直线; fzero 会简直爆炸!因此,在进入 fzero 之前,必须尽可能接近解决方案 .
所以,这是一种获得 good 初始值的方法 .
考虑这个功能
f(x) = x · tan(x) - 1
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·π
如您所见,近似值基本上等于解决方案 . 对应情节: