Метод ньютона реализация python

Нелинейные системы и уравнения

Ищется решение нелинейного уравнения $$ \begin \tag f(x) = 0, \end $$ где \( f(x) \) — заданная функция. Корни уравнения (1) могут быть комплексными и кратными. Выделяют как самостоятельную проблему разделения корней, когда проводится выделение области в комплексной плоскости, содержащей один корень. После этого на основе тех или иных итерационных процессов при выбранном начальном приближении находится решение нелинейного уравнения (1).

В более общем случае мы имеем не одно уравнение (1), а систему нелинейных уравнений $$ \begin \tag f_i(x_1, x_2, \ldots, x_n) = 0, \quad i = 1, 2, \ldots n. \end $$ Обозначим через \( \mathbf = (x_1, x_2, \ldots, x_n) \) вектор неизвестных и определим вектор-функцию \( \mathbf(\mathbf) = (f_1(\mathbf), f_2(\mathbf), \ldots, f_n(\mathbf)) \). Тогда система (2) записывается в виде $$ \begin \tag \mathbf(\mathbf) = 0. \end $$ Частным случаем (3) является уравнение (1) (\( n = 1 \)). Второй пример (3) — система линейных алгебраических уравнений, когда \( \mathbf (\mathbf) = A \mathbf — \mathbf \).

Метод Ньютона

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

При итерационном решении уравнений (1), (3) задается некоторое начальное приближение, достаточно близкое к искомому решению \( x^* \). В одношаговых итерационных методах новое приближение \( x_ \) определяется по предыдущему приближению \( x_k \). Говорят, что итерационный метод сходится с линейной скоростью, если \( x_ — x^* = O(x_k — x^*) \) и итерационный метод имеет квадратичную сходимость, если \( x_ — x^* = O(x_k — x^*)^2 \).

В итерационном методе Ньютона (методе касательных) для нового приближения имеем $$ \begin \tag x_ = x_k + \frac, \quad k = 0, 1, \ldots, \end $$

Вычисления по (4) проводятся до тех пор, пока \( f(x_k) \) не станет близким к нулю. Более точно, до тех пор, пока \( |f_(x_k)| > \varepsilon \), где \( \varepsilon \) — малая величина.

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

def naive_Newton(f, dfdx, x, eps): while abs(f(x)) > eps: x = x - float(f(x))/dfdx(x) return x

Чтобы найти корень уравнения \( x^2 = 9 \) необходимо реализовать функции

def f(x): return x**2 - 9 def dfdx(x): return 2*x print naive_Newton(f, dfdx, 1000, 0.001)

Данная функция хорошо работает для приведенного примера. Однако, в общем случае могут возникать некоторые ошибки, которые нужно отлавливать. Например: пусть нужно решить уравнение \( \tanh(x) = 0 \), точное решение которого \( x = 0 \). Если \( |x_0| \leq 1.08 \), то метод сходится за шесть итераций.

-1.0589531343563485 0.9894042072982367 -0.784566773085775 0.3639981611100014 -0.03301469613719421 2.3995252668003453e-05 Число вызовов функций: 25 Решение: 3.0000000001273204

Теперь зададим \( x_0 \) близким к \( 1.09 \). Возникнет переполнение

-1.0933161820201083 1.104903543244409 -1.1461555078811896 1.3030326182332865 -2.064923002377556 13.473142800575976 -126055892892.66042

Возникнет деление на ноль, так как для \( x_7 = -126055892892.66042 \) значение \( \tanh(x_7) \) при машинном округлении равно \( 1.0 \) и поэтому \( f^\prime(x_7) = 1 — \tanh(x_7)^2 \) становится равной нулю в знаменателе.

Проблема заключается в том, что при таком начальном приближении метод Ньютона расходится.

Еще один недостаток функции naive_Newton заключается в том, что функция f(x) вызывается в два раза больше, чем необходимо.

1. обрабатывать деление на ноль
2. задавать максимальное число итераций в случае расходимости метода
3. убрать лишний вызов функции f(x)

def Newton(f, dfdx, x, eps): f_value = f(x) iteration_counter = 0 while abs(f_value) > eps and iteration_counter  100: try: x = x - f_value/dfdx(x) except ZeroDivisionError as err: print("Ошибка: <>".format(err)) sys.exit(1) f_value = f(x) iteration_counter += 1 if abs(f_value) > eps: iteration_counter = -1 return x, iteration_counter

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

При реализации метода Ньютона нужно знать аналитическое выражение для производной \( f^\prime(x) \). Python содержит пакет SymPy, который можно использовать для создания функции dfdx . Для нашей задачи это можно реализовать следующим образом:

from sympy import symbol, tanh, diff, lambdify x = symbol.Symbol('x') # определяем математический символ x f_expr = tanh(x) # символьное выражение для f(x) dfdx_expr = diff(f_expr, x) # вычисляем символьно f'(x) # Преобразуем f_expr и dfdx_expr в обычные функции Python f = lambdify([x], # аргумент f f_expr) dfdx = lambdify([x], dfdx_expr) tol = 1e-1 sol, no_iterations = Newton(f, dfdx, x=1.09, eps=tol) if no_iterations > 0: print("Уравнение: tanh(x) = 0. Число итераций: <>".format(no_iterations)) print("Решение: <>, eps = <>".format(sol, tol)) else: print("Решение не найдено!") def F(x): y = np.zeros_like(x) y[0] = (3 + 2*x[0])*x[0] - 2*x[1] - 3 y[1:-1] = (3 + 2*x[1:-1])*x[1:-1] - x[:-2] - 2*x[2:] -2 y[-1] = (3 + 2*x[-1])*x[-1] - x[-2] - 4 return y def J(x): n = len(x) jac = np.zeros((n, n)) jac[0, 0] = 3 + 4*x[0] jac[0, 1] = -2 for i in range(n-1): jac[i, i-1] = -1 jac[i, i] = 3 + 4*x[i] jac[i, i+1] = -2 jac[-1, -2] = -1 jac[-1, -1] = 3 + 4*x[-1] return jac n = 10 guess = 3*np.ones(n) sol, its = Newton_system(F, J, guess) if its > 0: print("x = <>".format(sol)) else: print("Решение не найдено!")

Решение нелинейных систем

Идея метода Ньютона для приближенного решения системы (2) заключается в следующем: имея некоторое приближение \( \pmb^ \), мы находим следующее приближение \( \pmb^ \), аппроксимируя \( \pmb(\pmb^) \) линейным оператором и решая систему линейных алгебраических уравнений. Аппроксимируем нелинейную задачу \( \pmb(\pmb^) = 0 \) линейной $$ \begin \tag \pmb(\pmb^) + \pmb(\pmb^)(\pmb^ — \pmb^) = 0, \end $$ где \( \pmb(\pmb^) \) — матрица Якоби (якобиан): $$ \pmb(\pmb^) = \begin \frac<\partial f_1(\pmb^)> <\partial x_1>& \frac<\partial f_1(\pmb^)> <\partial x_2>& \ldots & \frac<\partial f_1(\pmb^)> <\partial x_n>\\ \frac<\partial f_2(\pmb^)> <\partial x_1>& \frac<\partial f_2(\pmb^)> <\partial x_2>& \ldots & \frac<\partial f_2(\pmb^)> <\partial x_n>\\ \vdots & \vdots & \ldots & \vdots \\ \frac<\partial f_n(\pmb^)> <\partial x_1>& \frac<\partial f_n(\pmb^)> <\partial x_2>& \ldots & \frac<\partial f_n(\pmb^)> <\partial x_n>\\ \end $$ Уравнение (5) является линейной системой с матрицей коэффициентов \( \pmb \) и вектором правой части \( -\pmb(\pmb^) \). Систему можно переписать в виде $$ \pmb(\pmb^)\pmb = — \pmb(\pmb^), $$ где \( \pmb = \pmb^ — \pmb^ \).

Таким образом, \( k \)-я итерация метода Ньютона состоит из двух стадий:

1. Решается система линейных уравнений (СЛАУ) \( \pmb(\pmb^)\pmb = -\pmb(\pmb^) \) относительно \( \pmb \).

2. Находится значение вектора на следующей итерации \( \pmb^ = \pmb^ + \pmb \).

Для решения СЛАУ можно использовать приближенные методы. Можно также использовать метод Гаусса. Пакет numpy содержит модуль linalg , основанный на известной библиотеке LAPACK, в которой реализованы методы линейной алгебры. Инструкция x = numpy.linalg.solve(A, b) решает систему \( Ax = b \) методом Гаусса, реализованным в библиотеке LAPACK.

Когда система нелинейных уравнений возникает при решении задач для нелинейных уравнений в частных производных, матрица Якоби часто бывает разреженной. В этом случае целесообразно использовать специальные методы для разреженных матриц или итерационные методы.

Можно также воспользоваться методами, реализованными для систем линейных уравнений.

Источник

Русские Блоги

[Численный анализ] Ньютон Метод — Реализация Python

Эта статья реализована для численного анализа кода курса — Newtona

# -*- coding: utf-8 -*- """ Ньютон @author: hhuaf """ import numpy as np import matplotlib.pyplot as plt # input ''' X0: начальное значение Тета: порог ''' X0 = float (вход ('входная начальная точка: (например, 5, 10, 15, 20 . ) \ n')) theta=1e-5 # Вы можете отобразить китайский plt.rcParams["font.sans-serif"] = ["SimHei"] plt.rcParams['axes.unicode_minus'] = False # Установить стиль plt.style.use('ggplot') # Определить функции init_fun = lambda x: x**2-4*x # Производное deri_fun = lambda x: 2*x-4 # Установите изображение fig_1 = plt.figure(figsize = (8, 6)) plt.hlines(0,-1,x0,'black','--') plt.xlabel('X') plt.ylabel('Y') Plt.title ('$ f (x) = x ^ 2-4x $ Image') # Функция изображения x=[] if x0>0: x = np.arange(-1,x0,0.05) plt.hlines(0,-1,x0,'black','--') else: x = np.arange(x0,10,0.05) plt.hlines(0,x0,10,'black','--') y = init_fun(x) # def Newton(func = init_fun, df = deri_fun, x0 = x0,theta = theta): number=0 xi = x0 while True and number 

Добро пожаловать всех к критике!

Интеллектуальная рекомендация

Просмотр номера порта

windows+R Затем введите CMD Затем введите netstat -a Или введите Netstat -ano Нажмите Диспетчер задач Нажмите Подробную информацию, чтобы просмотреть номер PID Согласно номеру порта, PID соответствует.

Веб-путешествие - знакомство с HTML5 первым (1)

Чтобы начать фронтенд-путешествие, мы должны сначала понять некоторые базовые знания о резерве. В этом разделе будет продемонстрирована простейшая композиция страницы HTML5 и представлены исторические.

Установка RocketMQ

1. Подготовка окружающей среды 64 -битный OS Linux JDK выше 1,8 Установите Maven Команда установки заключается в следующем: wget https://mirrors.tuna.tsinghua.edu.cn/apache/maven/maven-3/3.6.3/binarie.

java.net.BindException: Address already in use: bind

Исключение: Причина: порт занят. В моем проекте это происходит потому, что одна и та же программа сокетов запускается дважды. Это исключение отличается от того, что tomcat запускал неско.

Установить svn под linux

Выполните следующую команду, чтобы установить SVN. Выполните следующую команду, чтобы просмотреть версию SVN. Выполните следующие действия, чтобы создать репозиторий: Выполните следующую команду, чтоб.

Вам также может понравиться

iOS12.1 использует UINavigationController + UITabBarController (UITabBar matte), после установки hidesBottomBarWhenPasted, .

Если вы используете систему OS12.1 UINavigationController + UITabBarController (UITabBar matte), вы столкнетесь с проблемой беспорядка макета панели вкладок в popViewControllerAnimated: Код, вызывающи.

переменные среды Java

jvm/jdk/jre 》 Jvm (виртуальная машина, компилировать код в класс) 》 Jdk (укажите пакет классов, который должен вызывать Java) 》 Jre (среда, необходимая для выполнения) Переменные среды: 》 Путь: путь, .

Springboot интегрирует шаблон MongoDB

Предисловие Основным проектом в прошлом году была разработка и применение springboot, и я собираюсь немного использовать данные springboot. Дело не в том, чтобы записывать использование api (возможно.

CF 1165D почти все делители воды

http://codeforces.com/contest/1165/problem/D Главная мысль: Введите количество раз, для каждого, введите n, затем следуйте N номер N, запрашивайте их минимальные общепринятые несколько, а затем опреде.

Rand (), Srand () Случайный номер сборки

Система предоставляет две функции в функции библиотеки для генерации случайных чисел: SRAND () и RAND (); Файл заголовка Функция определения: int rand (void), Функция функции: генерир.

Источник