python求解非线性方程组

非线性方程组概念

非线性方程组是指包含两个或多个非线性方程的方程组。非线性方程是那些不是线性的方程，即它们不能表示为 $a x + b y = c$ 的形式，其中 $a$ 、 $b$ 和 $c$ 是常数， $x$ 和 $y$ 是变量。非线性方程可以包含变量的平方、立方、乘积、商、指数、对数等。

例如，以下是一个非线性方程组：

$\begin{cases} x^2 + y^2 = 4 \\ xy = 1 \end{cases}$

非线性方程组的数值解法

解非线性方程组通常比解线性方程组更复杂，因为它们可能有多个解，或者没有解。解非线性方程组的方法包括代入法、消元法、图形法、数值方法等。

牛顿法（Newton’s Method）：
- 牛顿法是一种迭代方法，它使用函数的一阶泰勒展开来逼近方程的根。
- 对于非线性方程组，牛顿法需要计算雅可比矩阵（Jacobian matrix）并求解线性方程组来更新解的近似值。
- 牛顿法的收敛速度通常很快，但需要一个良好的初始猜测，并且当雅可比矩阵接近奇异时可能不收敛。
拟牛顿法（Quasi-Newton Methods）：
- 拟牛顿法是牛顿法的一种变体，它不需要计算雅可比矩阵的逆。
- 这类方法通过近似雅可比矩阵的逆来更新解的近似值，常见的拟牛顿法有 Broyden 法和 DFP 法。
高斯-赛德尔迭代法（Gauss-Seidel Iteration）：
- 高斯-赛德尔迭代法是一种简单的迭代方法，它通过逐个更新方程组中每个方程的解来逼近方程组的根。
- 这种方法适用于对角占优的方程组，但收敛速度可能较慢。
雅可比迭代法（Jacobi Iteration）：
- 雅可比迭代法与高斯-赛德尔迭代法类似，但它在每次迭代中同时更新所有方程的解。
- 雅可比迭代法的收敛速度通常比高斯-赛德尔迭代法慢。
Levenberg-Marquardt 方法：
- Levenberg-Marquardt 方法是一种用于求解非线性最小二乘问题的方法，它结合了牛顿法和梯度下降法。
- 这种方法通过调整一个参数来控制步长，以确保收敛。
信赖域方法（Trust Region Methods）：
- 信赖域方法通过在每一步中定义一个“信赖域”来控制步长，以确保解的质量和收敛性。
- 这类方法通常用于求解优化问题，但也可以用于求解非线性方程组。
子空间方法（Subspace Methods）：
- 子空间方法在解空间的一个子空间内寻找解，通常用于大规模问题。
割线法（Secant Method）：
- 割线法是一种不需要计算导数的迭代方法，它通过使用函数值的差分来近似导数。
不动点迭代法（Fixed-Point Iteration）：
- 不动点迭代法通过反复应用一个函数来逼近不动点，即方程 ( g(x) = x ) 的解。
多变量割线法（Multivariate Secant Method）：
- 多变量割线法是割线法在多变量情况下的推广，它使用多个点来近似雅可比矩阵。

这些方法各有优缺点，选择哪种方法取决于具体问题的特性和计算资源的要求。在实际应用中，可能需要结合问题的具体特点和计算条件来选择最合适的数值解法。

牛顿法求解非线性方程组

牛顿法（Newton’s Method），也称为牛顿-拉夫逊方法（Newton-Raphson Method），是一种用于数值求解非线性方程组的迭代方法。下面是牛顿法求解非线性方程组的详细介绍：

基本思想

对于非线性方程组 $ F(x) = 0 $，其中 $ F: \mathbb{R}^n \rightarrow \mathbb{R}^n ) 和 ( x \in \mathbb{R}^n $，牛顿法通过在每次迭代中线性化问题来逼近方程组的根。

局部线性化

在当前迭代点 $ x_k $ 处，我们可以使用泰勒展开式来近似 $F(x) $：
$\approx F(x_k) + J(x_k)(x - x_k)$
其中，$ J(x_k) $ 是 $ F$ 在 $ x_k $ 处的雅可比矩阵。

迭代公式推导

为了找到 $ F(x) = 0 $ 的根，我们希望找到一个 $ x $ 使得：
$F(x_k) + J(x_k)(x - x_k) = 0$
解这个线性方程组，我们得到：
$J(x_k)(x - x_k) = -F(x_k)$
$x - x_k = -J(x_k)^{-1} F(x_k)$
$x = x_k - J(x_k)^{-1} F(x_k)$
这就是牛顿法的迭代公式。

算法流程

取初始点 $ x^{(0)}$，设定最大迭代次数 $ N $ 和精度要求 $ \varepsilon $，置 $ k=0 $；
计算当前点 $ x_k $ 处的函数值 $ F(x_k) $ 和雅可比矩阵 $ J(x_k) $；
如果 $ |F(x_k)| < \varepsilon $，则停止迭代，输出 $ x_k $ 作为方程组的解；
否则，计算增量 $ \Delta x_k = -J(x_k)^{-1} F(x_k) $；
更新解 $ x_{k+1} = x_k + \Delta x_k $；
$ k = k + 1 $，返回步骤 2。

python 实现

下面是两个使用牛顿法的原理来手动实现求解非线性方程组的示例：
$\begin{cases} x^2 + y - 1 = 0 \\ x - y^2 - 1 = 0 \end{cases}$

import numpy as np# 定义非线性方程组
def F(x):return np.array([x[0]**2 + x[1] - 1,x[0] - x[1]**2 - 1])# 定义雅可比矩阵
def J(x):return np.array([[2*x[0], 1],[1, -2*x[1]]])# 牛顿法迭代求解非线性方程组
def newton_raphson(F, J, x0, tol=1e-8, max_iter=100):x = x0.copy()  # 创建初始猜测值的副本for i in range(max_iter):F_val = F(x)J_val = J(x)# 检查是否满足容差条件if np.linalg.norm(F_val) < tol:return x, i+1# 计算增量delta = np.linalg.solve(J_val, -F_val)# 更新解x = x + delta  # 这里应该是加法，而不是赋值return x, max_iter# 初始猜测值
x0 = np.array([1.0, 1.0])# 调用牛顿法求解函数
solution, iterations = newton_raphson(F, J, x0)print(f'解：{solution}')
print(f'迭代次数：{iterations}')

解：[1.00000000e+00 1.57081035e-11]
迭代次数：7

$\begin{cases} x^2 + y^2 - x = 0 \\ x^2 - y^2 - y = 0 \end{cases}$

import numpy as np# 定义非线性方程组
def F(x):return np.array([x[0]**2 + x[1]**2 - x[0],x[0]**2 - x[1]**2 - x[1]])# 定义雅可比矩阵
def J(x):return np.array([[2*x[0]-1, 2*x[1]],[2*x[0], -2*x[1]-1]])# 牛顿法迭代求解非线性方程组
def newton_raphson(F, J, x0, tol=1e-6, max_iter=100):x = x0.copy()  # 创建初始猜测值的副本for i in range(max_iter):F_val = F(x)J_val = J(x)# 检查是否满足容差条件if np.linalg.norm(F_val) < tol:return x, i+1# 计算增量delta = np.linalg.solve(J_val, -F_val)# 更新解x = x + delta  # 这里应该是加法，而不是赋值return x, max_iter# 初始猜测值
x0 = np.array([0.76, 0.4])# 调用牛顿法求解函数
solution, iterations = newton_raphson(F, J, x0)print(f'解：{solution}')
print(f'迭代次数：{iterations}')

解：[0.77184473 0.4196436 ]
迭代次数：3