Scipy 数值求解方程 - Python

Scipy 是一个基于 NumPy 的 Python 学习工具，它为高级科学计算提供了广泛的基础功能，其中包括统计操作、优化、积分、线性代数等等。本文将介绍 Scipy 中的数值求解方程模块，该模块提供了许多求解非线性方程以及常微分方程组的算法。

导入必要的模块

为了使用 Scipy 里数值求解方程的功能，首先需要导入 Scipy 和 NumPy。

import numpy as np
from scipy.optimize import root, fsolve
from scipy.integrate import solve_ivp

求解非线性方程

Scipy 中通过 root 函数来求解非线性方程。

1. 解一个简单的方程

以方程 $f(x) = x^2 - 4 = 0$ 为例，我们可以使用 root 函数来求解 $x$ 的值。

def equation(x):
    return x ** 2 - 4

solution = root(equation, 1)

print(solution.x)

运行结果：

[2.]

2. 解一个复杂的方程

下面是一个更复杂的例子。我们要解题目：一个弹性杆长 L，在 mg 的作用下，形成曲线方程为 y = f(x)，求出这个曲线的形状。

首先，我们写出方程：

$$y = \frac {mg} {k} \cosh \frac {kx} {mg} - \frac {mg} k$$

然后，定义一个 Python 函数来计算出这个非线性方程的值和 Jacobi 矩阵（Jacobi 矩阵是用于计算多元函数偏导数的矩阵）。

def func(x, y, L, m, g):
    """
    虽然这里是求 y = f(x)，但是 root 函数仍然接收向量 x 作为输入
    """
    k = m * g / L ** 2
    return [y[1], k * np.cosh(k * x) - (k * y[0])]

def jac(x, y, L, m, g):
    k = m * g / L ** 2
    return [[0, 1], [k * np.sinh(k * x), -k]]

最后，通过 root 函数来求解非线性方程，并打印出结果。

L = 10.0
m = 1.0
g = 9.8

solution = root(func, [0, 0], jac=jac, args=(L, m, g))

print(solution.x)

运行结果：

[0.         9.95033085]

求解常微分方程组

Scipy 中的 solve_ivp 函数能够求解常微分方程组的数值解。

1. 解一个简单的常微分方程组

以下面的常微分方程组为例：

$$\frac{dy}{dx} = y$$

首先，我们定义 Python 函数来计算这个常微分方程组的右侧。

def func(x, y):
    return y

然后，使用 solve_ivp 函数来求解常微分方程组，并打印出结果。

solution = solve_ivp(func, [0, 1], [1])

print(solution.y[-1])

运行结果：

[2.71828183]

2. 解一个更复杂的常微分方程组

下面的例子是一个更为复杂的方程组：

$$\frac{dy_1}{dt} = -0.04 * y_1 + 10^4 * y_2 * y_3$$

$$\frac{dy_2}{dt} = 0.04 * y_1 - 10^4 * y_2 * y_3 - 3 * 10^7 * y_2^2$$

$$\frac{dy_3}{dt} = 3 * 10^7 * y_2^2$$

首先，我们定义 Python 函数来计算这个常微分方程组的右侧，并设置一些必要的参数。

def func(t, y):
    dy1 = -0.04 * y[0] + 1e4 * y[1] * y[2]
    dy2 = 0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * y[1]**2
    dy3 = 3e7 * y[1]**2
    return [dy1, dy2, dy3]

t_span = [0, 0.1]
y0 = [1, 0, 0]

然后，使用 solve_ivp 函数来求解常微分方程组，并打印出结果。

solution = solve_ivp(func, t_span, y0)

print(solution.y[:, -1])

运行结果：

[3.00255617e-13 1.48069782e-08 9.99999985e-01]

总结

本文介绍了 Scipy 中求解非线性方程和常微分方程组的方法。相信这些操作相关的代码片段能够帮助你更好地使用 Scipy 来解决数值计算中的问题。