Как я могу получить все решения этого уравнения в MATLAB?

Я хотел бы решить следующее уравнение: tan(x) = 1/x

Что я сделал:

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

Но это дает мне только одно численное приближение решения. После этого я прочитал о следующем:

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

Но это говорит мне, что он не может найти явное решение.

Как я могу найти численные решения этого уравнения в некотором заданном диапазоне?


person user137425    schedule 21.11.2016    source источник
comment
Что вы подразумеваете под получением всех значений? У этого уравнения есть бесконечное количество решений (просто постройте и посмотрите).   -  person Dev-iL    schedule 21.11.2016
comment
Извините, я имел в виду более одного значения. Например в ассортименте.   -  person user137425    schedule 21.11.2016
comment
Я имел в виду численные решения в некотором диапазоне   -  person user137425    schedule 21.11.2016


Ответы (3)


Чтобы найти численное решение функции в некотором диапазоне, вы можете использовать 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

Обратите внимание, что:

  1. Функция симметрична, как и ее корни. Это означает, что вы можете решить только положительный интервал (например) и получить отрицательные корни «бесплатно».
  2. Каждая другая запись в sol неверна, потому что именно здесь у нас есть асимптотические разрывы (tan переходы от +Inf к -Inf), которые MATLAB ошибочно распознает как решение. Так что вы можете просто игнорировать их (например, sol = sol(1:2:end);.
person Dev-iL    schedule 21.11.2016
comment
Спасибо, но он показывает только одно решение. - person user137425; 21.11.2016
comment
Измените уравнение на x*sin(x)-cos(x) == 0, чтобы еще больше избавиться от единственного числа. - person Lutz Lehmann; 22.11.2016

Мне никогда не нравилось использовать решатели «вслепую», то есть без какой-то приличной схемы выбора начального значения. По моему опыту, ценности, которые вы обнаружите, делая что-то вслепую, также будут вне контекста. Это означает, что вы часто будете пропускать решения, думать, что что-то является решением, в то время как на самом деле решатель взорвался и т. д.

В этом конкретном случае важно понимать, что fzero использует численные производные, чтобы находить все более точные приближения. Но производные для f(x) = x · tan(x) - 1 становится все труднее точно вычислить при увеличении x:

введите здесь описание изображения

Как видите, чем больше становится 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·π 

Как видите, аппроксимация в основном равна решению. Соответствующий сюжет:

введите здесь описание изображения

person Rody Oldenhuis    schedule 21.11.2016
comment
Почему бы просто не умножить на оскорбительный член cos(x), чтобы найти бесполюсное уравнение x*sin(x)-cos(x)==0? - person Lutz Lehmann; 22.11.2016
comment
@LutzL хм .... и что потом? До сих пор трудно найти хорошие первоначальные оценки для этого. Это вектор, зависящий от угла, то есть он идентичен √(x²+1) · sin(x + atan2(-1,x)). Как видите, период не является чистым 2π, а зависит от x. Вы можете добавить n·π к первому решению, которое вы найдете, чтобы получить последующие оценки, и для этого случая я думаю, что это действительно достаточно хорошо. Но опыт показывает, что для других случаев, когда период точно не совпадает с тем, что вы предполагаете в первоначальных оценках, вы в конечном итоге дважды сойдетесь к одному и тому же решению и пропустите следующий один... - person Rody Oldenhuis; 22.11.2016
comment
@LutzL хорошо, вы можете добавить ответ, который сделает мой смехотворным, мне искренне нравится, когда это происходит :) - person Rody Oldenhuis; 22.11.2016
comment
Сделал именно это, выразил свое мнение в ответе. -- Вы имеете в виду модифицированные задачи tan(a*x)=1/x или еще какие-то дополнительные термины? - person Lutz Lehmann; 22.11.2016

Умножьте уравнение на 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=0 второй степени x^2-(1-0.5x^2)==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
person Lutz Lehmann    schedule 22.11.2016