Python 牛顿法求解多元非线性方程组

Dilkople

已于 2022-11-02 13:04:02 修改

阅读量2.9k

点赞数 3

文章标签： python 算法

于 2022-04-22 22:38:41 首次发布

本文链接：https://blog.csdn.net/dilkople/article/details/124356524

版权

该代码为通用代码，代码使用方式可见代码中示例代码。需安装numpy和sympy库。

牛顿法简单讲解：

$\left. \begin{array} { c } { x _ { n + 1 } = x _ { n } - \frac { \Delta _ { x } ^ { ( n ) } } { J ( x _ { n } , y _ { n } ) } } \\ { y _ { n + 1 } = y _ { n } - \frac { \Delta _ { y } ^ { ( n ) } } { J ( x _ { n } , y _ { n } ) } } \end{array} \right.$

其中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("x0:{}".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 byStep(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 byEpsilon(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 的方程组
    # byStep(x0, step): x0 为初始值 step 为迭代次数
    # byEpsilon(x0, epsilon): x0 为初始值 epsilon 为精度(当|前一步计算值-后一步计算值|<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.byStep([1, 1, 1], 5)
    newton.byEpsilon([1, 1, 1], 0.001)

    # 一元非线性使用方法
    equations = [lambda x: cos(x[0]) + sin(x[0])]

    newton = Newton(equations)
    newton.byStep([1], 5)
    newton.byEpsilon([1], 0.001)

个人博客地址：Python 牛顿法求解多元非线性方程组 – Pancake's Personal Website