代码使用方式可见代码中示例代码。需安装numpy和sympy库。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 |
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 newtonByStep(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 newtonByEpsilon(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 为 equation = 0 的方程组 # x0 为初始值 step 为迭代次数 epsilon 为精度 # 多元非线性使用方法 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.newtonByStep([1, 1, 1], 5) newton.newtonByEpsilon([1, 1, 1], 0.001) # 一元非线性使用方法 equations = [lambda x: cos(x[0]) + sin(x[0])] newton = Newton(equations) newton.newtonByStep([1], 5) newton.newtonByEpsilon([1], 0.001) # 考试题 equations = [lambda x: 2 * x[0] ** 3 - x[1] ** 2 - 1, lambda x: x[0] * x[1] ** 3 - x[1] - 4] newton = Newton(equations) newton.newtonByStep([1.2, 1.7], 5) |
浙公网安备 33010902002953号