该代码为通用代码，代码使用方式可见代码中示例代码。需安装numpy和sympy库。
牛顿法简单讲解：
其中J为雅克比行列式，不了解牛顿法的朋友可以先去了解一下计算过程。
from math import *
from sympy import *
import sympy
import numpy as np
class Newton:
def __init__(self, equations):
self.equations = equations
self.n = len(equations)
self.x = sympy.symbols(" ".join("x{}".format(i) for i in range(self.n)) + " x{}".format(self.n), real=True)
self.equationsSymbol = [equations[i](self.x) for i in range(self.n)]
# 初始化 Jacobian 矩阵
self.J = np.zeros(self.n * self.n, dtype=sympy.core.add.Add).reshape(self.n, self.n)
for i in range(self.n):
for j in range(self.n):
self.J[i][j] = sympy.diff(self.equationsSymbol[i], self.x[j])
def cal_J(self, x):
dict = {self.x[i]: x[i] for i in range(self.n)}
J = np.zeros(self.n * self.n).reshape(self.n, self.n)
for i in range(self.n):
for j in range(self.n):
J[i][j] = self.J[i][j].subs(dict)
return J
def cal_f(self, x):
f = np.zeros(self.n)
for i in range(self.n):
f[i] = self.equations[i](x)
f.reshape(self.n, 1)
return f
def interationsByStep(self, x0, step):
x0 = np.array(x0)
for i in range(step):
x0 = x0 - np.linalg.pinv(self.cal_J(x0)) @ self.cal_f(x0)
print("Step {}:".format(i + 1), ", ".join(["x{} = {}".format(j + 1, x0[j]) for j in range(self.n)]))
return x0
def interationsByEpsilon(self, x0, epsilon):
error = float("inf")
while error >= epsilon:
cal = np.linalg.pinv(self.cal_J(x0)) @ self.cal_f(x0)
error = max(abs(cal))
x0 = x0 - cal
print(x0)
return x0
if __name__ == "__main__":
# 多元非线性使用方法
equations = [lambda x: cos(0.4 * x[1] + x[0] ** 2) + x[0] ** 2 + x[1] ** 2 - 1.6,
lambda x: 1.5 * x[0] ** 2 - x[1] ** 2 / 0.36 - 1,
lambda x: 3 * x[0] + 4 * x[1] + 5 * x[2]]
newton = Newton(equations)
newton.interationsByStep([1, 1, 1], 5)
newton.interationsByEpsilon([1, 1, 1], 0.001)
# 一元非线性使用方法
equations = [lambda x: cos(x[0]) + sin(x[0])]
newton = Newton(equations)
newton.interationsByStep([1], 5)
newton.interationsByEpsilon([1], 0.001)
个人博客地址：Python 牛顿法求解多元非线性方程组 – Pancake's Personal Website
