Мне никогда не нравилось использовать решатели «вслепую», то есть без какой-то приличной схемы выбора начального значения. По моему опыту, ценности, которые вы обнаружите, делая что-то вслепую, также будут вне контекста. Это означает, что вы часто будете пропускать решения, думать, что что-то является решением, в то время как на самом деле решатель взорвался и т. д.
В этом конкретном случае важно понимать, что fzero
использует численные производные, чтобы находить все более точные приближения. Но производные для f(x) = x · tan(x) - 1
становится все труднее точно вычислить при увеличении x
:
![введите здесь описание изображения](https://i.stack.imgur.com/01Cp8.png)
Как видите, чем больше становится x
, тем лучше f(x)
приближается к вертикальной линии; fzero
просто взорвется! Поэтому крайне важно получить оценку как можно близкую к решению, прежде чем даже вводить fzero
.
Итак, вот способ получить хорошие начальные значения.
Рассмотрим функцию
f(x) = x · tan(x) - 1
Зная, что tan(x)
имеет разложение Тейлора:
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·π
Как видите, аппроксимация в основном равна решению. Соответствующий сюжет:
![введите здесь описание изображения](https://i.stack.imgur.com/N5VOh.png)
person
Rody Oldenhuis
schedule
21.11.2016