Решить уравнение методом ньютона python

Метод Ньютона, решение нелинейных уравнений на Python

Метод Ньютона, также известный как метод касательных, является численным методом для нахождения корня нелинейного уравнения. Он основан на идее использования локального линейного приближения функции в окрестности предполагаемого корня для последующего уточнения его координат.

Алгоритм метода Ньютона для решения нелинейных уравнений

Процесс решения нелинейного уравнения методом Ньютона состоит из следующих шагов:

1. Задать начальное приближение корня x0.
2. Выразить функцию f(x) в виде локального линейного приближения в окрестности x0, используя формулу Тейлора: f(x) ≈ f(x0) + f'(x0) * (x – x0)
3. Решить полученное линейное уравнение относительно x, положив его равным нулю: 0 = f(x0) + f'(x0) * (x – x0)
4. Вычислить новое приближение корня, используя полученное решение: x1 = x0 – f(x0) / f'(x0)
5. Повторять шаги 2-4, пока не будет достигнута необходимая точность или заданное количество итераций.

Преимущества и недостатки метода Ньютона

Заметим, что метод Ньютона может сходиться к корню очень быстро, если начальное приближение близко к корню и функция имеет достаточно гладкую поверхность. Однако, если начальное приближение далеко от корня или функция имеет слишком сложную поверхность, метод может сходиться очень медленно или даже расходиться. Также важно учитывать, что метод Ньютона может находить только один корень в заданной окрестности начального приближения.

Пошаговая реализация метода Ньютона на Python

• Определение функции newton_method , которая принимает на вход функцию f , её производную df , начальное приближение x0 , допустимую погрешность tol и максимальное число итераций max_iter . Функция решает нелинейное уравнение f(x) = 0 с помощью метода Ньютона.

def newton_method(f, df, x0, tol=1e-6, max_iter=100):

• В теле функции newton_method создается переменная x , которая инициализируется значением начального приближения x0 .

• В цикле for выполняются итерации метода Ньютона. Цикл будет продолжаться до тех пор, пока не будет достигнута максимальная допустимая погрешность tol или не будет достигнуто максимальное число итераций max_iter .

• Затем проверяется, достигнута ли требуемая точность. Если значение fx меньше, чем tol , значит, мы нашли корень уравнения, и функция возвращает значение x и число итераций i .

• Если требуемая точность не достигнута, то вычисляется производная df в точке x и сохраняется в переменной dfx .

• Затем проверяется, равна ли производная dfx нулю. Если равна, значит, мы не можем продолжать итерации, и функция возвращает None и число итераций i .

• Если производная dfx не равна нулю, то вычисляется новое значение x с помощью формулы метода Ньютона: x = x — fx / dfx .

• Цикл завершается, и если мы не нашли корень уравнения, то функция возвращает None и число итераций i .

Метод Ньютона, готовый код

import math def newton_method(f, df, x0, tol=1e-6, max_iter=100): """ Решает нелинейное уравнение f(x) = 0 с помощью метода Ньютона. f - функция, для которой ищется корень df - производная функции f x0 - начальное приближение tol - допустимая погрешность max_iter - максимальное число итераций Возвращает корень уравнения и число итераций. """ x = x0 for i in range(max_iter): fx = f(x) if abs(fx) < tol: return x, i dfx = df(x) if dfx == 0: break x = x - fx / dfx return None, i # Пример использования функции newton_method # Решение уравнения x^2 - 2 = 0 f = lambda x: x**2 - 2 df = lambda x: 2*x x0 = 1.0 root, num_iter = newton_method(f, df, x0) if root is not None: print(f"Корень уравнения: ") print(f"Число итераций: ") else: print("Решение не найдено.")

Если функция newton_method нашла корень уравнения, то программа выводит его и число итераций. Если корень не найден, то программа выводит сообщение о том, что решение не найдено.

Результат выполнения программы в консоли:

Источник

Решение нелинейных уравнений#

Подмодуль scipy.optimize также содержит в себе методы для поиска корней нелинейных уравнений и их систем.

from matplotlib import pyplot as plt def configure_matplotlib(): plt.rc('text', usetex=True) plt.rcParams["axes.titlesize"] = 28 plt.rcParams["axes.labelsize"] = 24 plt.rcParams["legend.fontsize"] = 24 plt.rcParams["xtick.labelsize"] = plt.rcParams["ytick.labelsize"] = 18 plt.rcParams["text.latex.preamble"] = r""" \usepackage[utf8]  \usepackage[english,russian]  \usepackage  """ configure_matplotlib() 

Поиск корня скалярной функции одного аргумента#

Функция scipy.optimize.root_scalar позволяет искать корни функции \(f\colon \mathbb \to \mathbb\) :

Функция root_scalar предоставляет доступ к разным методам поиска корней, таким как newton , bisect , secant и многим другим. Какие-то из этих методов ищут корень внутри отрезка bracket , а другие ищут корень, начиная с какого-то начального приближения x0 .

import numpy as np from scipy import optimize from matplotlib import pyplot as plt def f(x): return x**3 - 1 solution = optimize.root_scalar(f, bracket=[-10, 10], method="bisect") print(solution) x = np.linspace(-3, 3, 100) fig, ax = plt.subplots(figsize=(8, 8), layout="tight") ax.plot(x, f(x), label="$f(x)$") ax.scatter(solution.root, f(solution.root) , color="red", label="найденный корень") ax.axvline(solution.root, color="black") ax.axhline(0, color="black") ax.set_xlabel("$x$") ax.set_ylabel("$y$") ax.legend() 

converged: True flag: 'converged' function_calls: 46 iterations: 44 root: 1.0000000000002274

Метод бисекции сошелся за 44 итерации. Проверим метод Ньютона.

def fprime(x): return 3 * x**2 solution = optimize.root_scalar(f, fprime=fprime, x0=-2, method="newton") print(solution) fig, ax = plt.subplots(figsize=(8, 8), layout="tight") ax.plot(x, f(x), label="$f(x)$") ax.scatter(solution.root, f(solution.root) , color="red", label="найденный корень") ax.axvline(solution.root, color="black") ax.axhline(0, color="black") ax.set_xlabel("$x$") ax.set_ylabel("$y$") ax.legend() 

converged: True flag: 'converged' function_calls: 20 iterations: 10 root: 1.0

Метод Ньютона сошелся за 11 итераций.

Решение системы нелинейных уравнений.#

Функция optimize.root предназначена для поиска корней уравнений вида

где \(x\in\mathbb^n\) и \(F\colon \mathbb^n \to \mathbb^n\) многомерны, т.е. optimize.root решает системы вида

Рассмотрим поиск корня на примере функции

Матрица Якоби этого уравнения имеет вид

\[\begin \begin 1 + \frac(x_2 - x_1)^2 & \frac(x_2 - x_1)^2 \\ \frac(x_2 - x_1)^2 & 1 + \frac(x_2 - x_1)^2 \end, \end\]

а единственный действительный корень которой

В самом простом варианте достаточно передать методу optimize.root функцию левой части уравнения \(F\) и начальное приближение к корню.

import numpy as np from scipy import optimize from matplotlib import pyplot as plt def f(x): x = np.array(x) return np.array([ x[0] + 0.5 * (x[0] - x[1])**3 - 1.0, 0.5 * (x[1] - x[0])**3 + x[1] - 1.0 ]) solution = optimize.root(f, x0 = [0, 0]) print(solution) 

fjac: array([[-1.00000000e+00, -1.09073861e-12], [ 1.09073861e-12, -1.00000000e+00]]) fun: array([0., 0.]) message: 'The solution converged.' nfev: 5 qtf: array([-2.18114415e-12, -2.18114415e-12]) r: array([-1.00000000e+00, -2.18119967e-12, -1.00000000e+00]) status: 1 success: True x: array([1., 1.])

Если известно аналитическое выражение для производно, то лучше задействовать его и передать по параметру jac .

def jac(x): return np.array([ [1 + 1.5 * (x[0] - x[1])**2, -1.5 * (x[0] - x[1])**2], [-1.5 * (x[1] - x[0])**2, 1 + 1.5 * (x[1] - x[0])**2] ]) solution = optimize.root(f, jac=jac, x0 = [0, 0], method="hybr") print(solution) 

fjac: array([[-1., 0.], [ 0., -1.]]) fun: array([0., 0.]) message: 'The solution converged.' nfev: 2 njev: 1 qtf: array([1., 1.]) r: array([-1., -0., -1.]) status: 1 success: True x: array([1., 1.])

Источник

How to use the Newton's method in python ?

In numerical analysis, Newton's method (also known as the Newton–Raphson method), named after Isaac Newton and Joseph Raphson, is a method for finding successively better approximations to the roots (or zeroes) of a real-valued function. wikipedia. Example of implementation using python:

Solution 1

from scipy import misc def NewtonsMethod(f, x, tolerance=0.000001): while True: x1 = x - f(x) / misc.derivative(f, x) t = abs(x1 - x) if t < tolerance: break x = x1 return x def f(x): return (1.0/4.0)*x**3+(3.0/4.0)*x**2-(3.0/2.0)*x-2 x = 4 x0 = NewtonsMethod(f, x) print('x: ', x) print('x0: ', x0) print("f(x0) = ", ((1.0/4.0)*x0**3+(3.0/4.0)*x0**2-(3.0/2.0)*x0-2 )) 

x: 4 x0: 2.0000002745869883 f(x0) = 1.2356416165815176e-06 

Solution 2 (scipy)

from scipy.optimize import newton def f(x): return (1.0/4.0)*x**3+(3.0/4.0)*x**2-(3.0/2.0)*x-2 x = 4 x0 = newton(f, x, fprime=None, args=(), tol=1.48e-08, maxiter=50, fprime2=None) print('x: ', x) print('x0: ', x0) print("f(x0) = ", ((1.0/4.0)*x0**3+(3.0/4.0)*x0**2-(3.0/2.0)*x0-2 )) 

x: 4 x0: 2.000000000000008 f(x0) = 3.597122599785507e-14 

Plot the above figure using matplotlib

#!/usr/bin/env python from pylab import * t = arange(-6.0, 4.0, 0.01) s= t*t*t/4.0+3.0*t*t/4.0-3*t/2.0-2.0 ax = subplot(111) ax.plot(t, s) ax.scatter([-4,-1,2],[0,0,0]) ax.grid(True) ax.spines['left'].set_position('zero') ax.spines['right'].set_color('none') ax.spines['bottom'].set_position('zero') ax.spines['top'].set_color('none') ax.set_xlim(-6,6) ax.set_ylim(-20,20) text(-3.0, 12, r"$f(x)=(1/4)*X^3+(3/4)*X^2-(3/2)*X-2$", horizontalalignment='center', fontsize=8) plt.title("How to implement the Newton's method using python \n for finding the zeroes of a real-valued function",fontsize=10) plt.xticks(fontsize=8) plt.yticks(fontsize=8) plt.savefig('NewtonMethodExemple.png') show() 

Note: with numpy it is also possible to find the root of a polynomial with root, example the following polynomial has 3 roots:

>>> import numpy as np >>> coeff = [1.0/4.0,3.0/4.0,-3.0/2.0,-2] >>> np.roots(coeff) array([-4., 2., -1.]) 

References

Источник