Метод Ньютона – Рафсона невозможен, когда функция почти плоская - PullRequest
2 голосов
/ 12 октября 2019

Я пытаюсь вычислить значение х ^ х. Я рассчитал, что это производная, и яиспользуя метод Ньютона-Рафсона, чтобы найти корень f (x) = x ^ xa для заданного a.

Это функция, которую я использую для вычисления следующего приближения:

double xkp1(double xk, double a){
    double lnxp1 = log(xk) +1;
    return xk -( (fx(xk, a))  /  (exp(xk*log(xk)) * (log(xk) + 1)  ));
}

, где функция fx определяется следующим образом:

double fx(double x, double a){
    return exp(x*log(x))-a;
}

Теперь этот подход работает правильнотолько если начальное значение xk уже близко к корню. Если он отличается на + -0,5, xk просто взрывается до действительно высокого значения, а f (x) становится бесконечностью. Теперь я думаю, что здесь может быть не так - производная от x ^ x очень мала по сравнению с фактическим значением x ^ x, поэтому целое (fx(xk, a)) / (exp(xk*log(xk)) * (log(xk) + 1) ) просто становится + бесконечностью - касательная пересекает корень. Это означает, что мне нужна более высокая двойная точность, но возможно ли это как-нибудь? Если нет, можно ли что-то изменить в методе, чтобы сделать его применимым в этой ситуации?

1 Ответ

3 голосов
/ 12 октября 2019

Я предполагаю, что вас интересуют только положительные x значения, а не отрицательные целые числа (которые также имеют действительные значения).

У меня есть два решения, которые могут помочь вам в вашей проблеме. Во-первых, численная схема может быть более стабильной, если вы сначала используете преобразование журнала. Тогда ваше уравнение должно найти a такой, что x log(x) = log(a). Производная от x log(x) равна log(x) + 1, которая может быть более стабильной, чем ваша экспоненциальная функция.

В противном случае вы можете использовать свои знания в том смысле, что f(x) = x^x монотонно увеличивается на x>1/e для выполнения деления пополампоиск. Начните с x_l = 1/e и x_r = 2/e. Если f(x_l) > a, то решения нет. Затем, пока f(x_r) < a, установите x_l = x_r и x_r = c * x_r с c > 0 (или любой другой схемой увеличения x_r). Если у вас есть пара, x_l, x_r, которая следует за f(x_l) < a < f(x_r), вы можете начать поиск пополам, чтобы найти небольшой интервал, содержащий решение для f(x) = a. Как только этот интервал достаточно мал, вы можете начать с итераций метода Ньютона.

Вы можете даже объединить два вышеупомянутых метода и выполнить биссектрису для поиска x log(x) = log(a).

...