【数据建模】微分方程与动力系统

文章目录

  • 微分方程与动力系统
    • 1. 微分方程的理论基础
      • 1.1 函数、导数与微分
      • 1.2 一阶线性微分方程的解
      • 1.3 二阶常系数线性微分方程的解
    • 2. 使用python求解微分方程
      • 2.1 求解微分
      • 2.2 求解定积分
        • 2.2.1 quad函数求解
        • 2.2.2 梯型法则求解
    • 3. 使用Scipy和Sympy解微分方程
      • 3.1 使用sympy求解微分方程解析解
      • 3.2 求解微分方程的方法
        • 3.2.1 使用 SymPy 求解微分方程解析解
        • 2.2.2 使用 SciPy 求解微分方程数值解
          • 2.2.2.1使用 `odeint` 求解一阶微分方程
          • 2.2.2.2使用 `odeint` 求解高阶微分方程
          • 2.2.2.3 使用 `solve_ivp` 求解一阶微分方程
          • 2.2.2.4 `odeint`和`solve_ivp`方法对比
      • 2.3 偏微分方程的数值求解
        • 2.3.1 离散化区域
        • 2.3.2. 差分公式
        • 2.3.3 边界条件和初始条件
        • 2.3.4 构造代数方程组
        • 2.3.5 求解代数方程组
        • 2.3.6 后处理
      • 示例代码:一维热传导方程
    • 4.微分方程案例
    • 5. 差分方程案例
    • 6. 元胞自动机
      • 6.1 元胞自动机介绍
      • 6.2 元胞自动机的应用场景
      • 6.3 元胞自动机的优点
      • 示例:康威的“生命游戏”
    • 7. 数值计算方法与微分方程求解
      • 7.1 梯度下降法
      • 7.2 牛顿法
      • 7.3 Euler 法与 Runge-Kutta 法
        • 7.3.1 Euler 法
        • 7.3.2 Runge-Kutta 法
      • 7.4 Crank-Nicolson 法
    • 8.学习心得
    • 9. 对学习文档的一些疑问
      • 9.1 文档2.1.4中梯型求解代码是否有误

微分方程与动力系统

1. 微分方程的理论基础

微分方程在物理、工程、生物学等诸多领域有着广泛应用,它们用于描述许多自然现象和工程问题。掌握微分方程的基础理论对理解和解决实际问题至关重要。

1.1 函数、导数与微分

函数是数学中描述变量之间关系的基本概念。设 y = f ( x ) y = f(x) y=f(x) 是一个函数,表示变量 y y y 随变量 x x x 的变化关系。导数是函数变化率的度量,定义为:
f ′ ( x ) = lim ⁡ Δ x → 0 f ( x + Δ x ) − f ( x ) Δ x f'(x) = \lim_{\Delta x \to 0} \frac{f(x + \Delta x) - f(x)}{\Delta x} f(x)=Δx0limΔxf(x+Δx)f(x)

微分是导数在微小变化下的线性近似。对于函数 y = f ( x ) y = f(x) y=f(x),它的微分表示为:
d y = f ′ ( x ) d x dy = f'(x) dx dy=f(x)dx

示例:
考虑函数 y = x 2 y = x^2 y=x2,则其导数为 f ′ ( x ) = 2 x f'(x) = 2x f(x)=2x,对应的微分为:
d y = 2 x d x dy = 2x dx dy=2xdx

1.2 一阶线性微分方程的解

一阶线性微分方程的一般形式为:
d y d x + P ( x ) y = Q ( x ) \frac{dy}{dx} + P(x)y = Q(x) dxdy+P(x)y=Q(x)

它的解法包括求解齐次方程和非齐次方程。齐次方程的形式为:
d y d x + P ( x ) y = 0 \frac{dy}{dx} + P(x)y = 0 dxdy+P(x)y=0

其通解为:
y = C e − ∫ P ( x )   d x y = Ce^{-\int P(x) \, dx} y=CeP(x)dx
其中, C C C 是任意常数。

对于非齐次方程,可以通过引入积分因子 μ ( x ) = e ∫ P ( x )   d x \mu(x) = e^{\int P(x) \, dx} μ(x)=eP(x)dx 转化为:
μ ( x ) d y d x + μ ( x ) P ( x ) y = μ ( x ) Q ( x ) \mu(x) \frac{dy}{dx} + \mu(x) P(x) y = \mu(x) Q(x) μ(x)dxdy+μ(x)P(x)y=μ(x)Q(x)
d d x [ μ ( x ) y ] = μ ( x ) Q ( x ) \frac{d}{dx} [\mu(x) y] = \mu(x) Q(x) dxd[μ(x)y]=μ(x)Q(x)

积分得到通解:
y = 1 μ ( x ) ( ∫ μ ( x ) Q ( x )   d x + C ) y = \frac{1}{\mu(x)} \left( \int \mu(x) Q(x) \, dx + C \right) y=μ(x)1(μ(x)Q(x)dx+C)

示例:
解方程 d y d x + 2 y = e x \frac{dy}{dx} + 2y = e^x dxdy+2y=ex

  1. 求积分因子 μ ( x ) = e ∫ 2   d x = e 2 x \mu(x) = e^{\int 2 \, dx} = e^{2x} μ(x)=e2dx=e2x
  2. 变换方程: e 2 x d y d x + 2 e 2 x y = e x e 2 x e^{2x} \frac{dy}{dx} + 2 e^{2x} y = e^{x} e^{2x} e2xdxdy+2e2xy=exe2x
  3. 简化为: d d x [ e 2 x y ] = e 3 x \frac{d}{dx} [e^{2x} y] = e^{3x} dxd[e2xy]=e3x
  4. 积分得到: e 2 x y = ∫ e 3 x   d x = e 3 x 3 + C e^{2x} y = \int e^{3x} \, dx = \frac{e^{3x}}{3} + C e2xy=e3xdx=3e3x+C
  5. 最终解为: y = e x 3 + C e − 2 x y = \frac{e^x}{3} + Ce^{-2x} y=3ex+Ce2x

1.3 二阶常系数线性微分方程的解

二阶常系数线性微分方程的一般形式为:
a d 2 y d x 2 + b d y d x + c y = 0 a \frac{d^2y}{dx^2} + b \frac{dy}{dx} + c y = 0 adx2d2y+bdxdy+cy=0

其特征方程为:
a r 2 + b r + c = 0 ar^2 + br + c = 0 ar2+br+c=0

根据特征方程的根,方程的解有三种情况:

  1. 两个实根 r 1 r_1 r1 r 2 r_2 r2
    y = C 1 e r 1 x + C 2 e r 2 x y = C_1 e^{r_1 x} + C_2 e^{r_2 x} y=C1er1x+C2er2x

  2. 一个实根 r r r
    y = ( C 1 + C 2 x ) e r x y = (C_1 + C_2 x)e^{rx} y=(C1+C2x)erx

  3. 一对共轭复根 r = α ± β i r = \alpha \pm \beta i r=α±βi
    y = e α x ( C 1 cos ⁡ ( β x ) + C 2 sin ⁡ ( β x ) ) y = e^{\alpha x} (C_1 \cos(\beta x) + C_2 \sin(\beta x)) y=eαx(C1cos(βx)+C2sin(βx))

示例:
解方程 y ′ ′ − 3 y ′ + 2 y = 0 y'' - 3y' + 2y = 0 y′′3y+2y=0

  1. 写出特征方程: r 2 − 3 r + 2 = 0 r^2 - 3r + 2 = 0 r23r+2=0
  2. 求根: r 1 = 1 r_1 = 1 r1=1, r 2 = 2 r_2 = 2 r2=2
  3. 通解为: y = C 1 e x + C 2 e 2 x y = C_1 e^x + C_2 e^{2x} y=C1ex+C2e2x

通过对函数、导数与微分的理解,以及一阶和二阶常系数线性微分方程的解法,我们可以解决许多实际问题。在实际应用中,通常会涉及到初值条件和边值条件的微分方程,需要根据具体问题进一步求解。

2. 使用python求解微分方程

在Python中,我们可以使用Numpy和SciPy这两个库来进行函数的微分和积分计算。

2.1 求解微分

在此举例:求解函数f(x) = x^2在点x = 1处的导数:

import numpy as np
# 定义x的取值范围和步长
x = np.linspace(0, 2, 100)
y = x**2
# 计算导数
dydx = np.gradient(y, x)
# 在x=1处的导数值
derivative_at_1 = dydx[np.argmin(abs(x - 1))]
print(f'在x=1处的导数值是:{derivative_at_1}')

np.gradient 函数计算数组 y 的梯度,返回 y 对 x 的导数。这个函数采用有限差分方法来近似计算导数。

np.argmin(abs(x - 1)) 计算 x 数组中与 1 最接近的值的索引。通过这个索引,我们可以找到对应的导数值 dydx。

2.2 求解定积分

在此举例:计算函数f(x) = cos(2πx) * exp(-x) + 1.2在区间[0, 0.7]上的定积分。

import numpy as np
from scipy.integrate import quad
# 定义函数
def f(x):
    return np.cos(2 * np.pi * x) * np.exp(-x) + 1.2
2.2.1 quad函数求解

使用SciPy库中的quad函数求解定积分:

def test():
    # 计算定积分
    integral, error = quad(f, 0, 0.7)
    print(f'定积分的结果是:{integral}')
2.2.2 梯型法则求解

使用数值积分的方法来近似计算,在此以梯型法则为例编写程序:

def tixing():
    # h = x[1]-x[0];
    # 积分区间
    x0 = 0
    xn = 0.7

    # 步长
    n = 1000  # 分成1000个小区间
    h = (xn - x0) / n

    s = 0
    for i in range(n):
        xn = x0 + i * h
        xn1 = xn + h;
        yn = f(xn)
        yn1 = f(xn1)
        s0 = (yn + yn1) * h / 2
        s += s0
    print(s)

3. 使用Scipy和Sympy解微分方程

3.1 使用sympy求解微分方程解析解

背景:大多数微分方程是求不出解析解的,只能在不同取值条件下求一个数值解。

3.2 求解微分方程的方法

在解决微分方程时,我们可以使用解析方法或数值方法。解析方法能给出精确的解,数值方法则在解析解难以获得时提供近似解。下面我们将分别使用 SymPySciPy 来求解微分方程。

3.2.1 使用 SymPy 求解微分方程解析解

SymPy 是一个用于符号计算的 Python 库,能够求解各种类型的微分方程。

示例:求解一阶线性微分方程

考虑一阶线性微分方程:
d y d x + y = e x \frac{dy}{dx} + y = e^x dxdy+y=ex

使用 SymPy 求解:

import sympy as sp

# 定义符号变量
x = sp.symbols('x')
y = sp.Function('y')(x)

# 定义微分方程
ode = sp.Eq(y.diff(x) + y, sp.exp(x))

# 求解微分方程
solution = sp.dsolve(ode, y)
print(solution)

解释

  1. 导入 SymPy 库。
  2. 定义符号变量 x 和函数 y(x)
  3. 定义微分方程 (\frac{dy}{dx} + y = e^x)。
  4. 使用 dsolve 函数求解微分方程,得到解析解。

示例:求解二阶常系数线性微分方程

考虑二阶常系数线性微分方程:
y ′ ′ − 3 y ′ + 2 y = 0 y'' - 3y' + 2y = 0 y′′3y+2y=0

使用 SymPy 求解:

import sympy as sp

# 定义符号变量
x = sp.symbols('x')
y = sp.Function('y')(x)

# 定义微分方程
ode = sp.Eq(y.diff(x, x) - 3*y.diff(x) + 2*y, 0)

# 求解微分方程
solution = sp.dsolve(ode, y)
print(solution)

解释

  1. 导入 SymPy 库。
  2. 定义符号变量 x 和函数 y(x)
  3. 定义二阶微分方程 y ′ ′ − 3 y ′ + 2 y = 0 y'' - 3y' + 2y = 0 y′′3y+2y=0
  4. 使用 dsolve 函数求解微分方程,得到解析解。
2.2.2 使用 SciPy 求解微分方程数值解

SciPy 是一个用于科学和技术计算的 Python 库,它包含许多用于数值求解微分方程的工具。Python求解微分方程数值解可以使用scipy库中的integrate包。在这当中有两个重要的函数:odeintsolve_ivp。但本质上,从底层来讲求解微分方程数值解的核心原理都是Euler 法和Runge-Kutta 法。

2.2.2.1使用 odeint 求解一阶微分方程

odeint()函数需要至少三个变量,第一个是微分方程函数,第二个是微分方程初值,第三个是微分的自变量。

示例:使用 odeint 求解一阶微分方程

考虑一阶微分方程:
d y d x = − 2 y + 1 \frac{dy}{dx} = -2y + 1 dxdy=2y+1
初始条件 y ( 0 ) = 0 y(0) = 0 y(0)=0

使用 SciPy 求解:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# 定义函数
def model(y, x):
    dydx = -2*y + 1
    return dydx

# 初始条件
y0 = 0

# x 的取值范围
x = np.linspace(0, 5, 100)

# 求解微分方程
y = odeint(model, y0, x)

# 绘制结果
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of dy/dx = -2y + 1')
plt.show()

解释

  1. 导入必要的库:numpy, odeintmatplotlib
  2. 定义微分方程的函数 model,返回导数 dydx
  3. 设置初始条件 y0 = 0
  4. 定义 x 的取值范围。
  5. 使用 odeint 函数求解微分方程。
  6. 使用 matplotlib 绘制解的曲线。
2.2.2.2使用 odeint 求解高阶微分方程

使用 odeint 求解高阶微分方程
Python求解微分方程数值解的时候是无法直接求解高阶微分方程的,必须通过换元降次的方法实现低阶化,把一个高阶微分方程替换成若干个一阶微分方程组成的微分方程组才能求解。

odeint 函数可以用于求解高阶微分方程。为了使用 odeint 求解高阶微分方程,我们需要将其转换为一组一阶微分方程。这里是一个示例,展示如何使用 odeint 求解二阶微分方程:

考虑以下二阶微分方程:
y ′ ′ − 3 y ′ + 2 y = 0 y'' - 3y' + 2y = 0 y′′3y+2y=0
初始条件为 y ( 0 ) = 1 y(0) = 1 y(0)=1 y ′ ( 0 ) = 0 y'(0) = 0 y(0)=0

步骤:

  1. 将二阶方程转换为一阶方程组:
    y 1 = y y_1 = y y1=y y 2 = y ′ y_2 = y' y2=y。则原方程可以转换为以下系统:
    { y 1 ′ = y 2 y 2 ′ = 3 y 2 − 2 y 1 \begin{cases} y_1' = y_2 \\ y_2' = 3y_2 - 2y_1 \end{cases} {y1=y2y2=3y22y1

  2. 定义初始条件:
    y 1 ( 0 ) = 1 y_1(0) = 1 y1(0)=1 y 2 ( 0 ) = 0 y_2(0) = 0 y2(0)=0

  3. 使用 odeint 进行求解。

下面是具体的代码示例:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# 定义微分方程组
def model(y, t):
    y1, y2 = y
    dy1dt = y2
    dy2dt = 3*y2 - 2*y1
    return [dy1dt, dy2dt]

# 初始条件
y0 = [1, 0]

# 时间点
t = np.linspace(0, 5, 100)

# 求解微分方程组
solution = odeint(model, y0, t)

# 提取解
y1 = solution[:, 0]
y2 = solution[:, 1]

# 绘制结果
plt.plot(t, y1, label='y(t)')
plt.plot(t, y2, label="y'(t)")
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.title("Solution of the differential equation y'' - 3y' + 2y = 0")
plt.show()

代码解释:

  1. 导入库:导入 numpyodeintmatplotlib 库。
  2. 定义微分方程组:将二阶微分方程转化为两个一阶微分方程,并定义为函数 model
  3. 初始条件:设置初始条件 y 1 ( 0 ) = 1 y_1(0) = 1 y1(0)=1 y 2 ( 0 ) = 0 y_2(0) = 0 y2(0)=0
  4. 时间点:定义时间范围 [0, 5],并生成 100 个等间距点。
  5. 求解微分方程组:使用 odeint 函数求解,并返回解数组。
  6. 提取解:从解数组中提取 y 1 ( t ) y_1(t) y1(t) y 2 ( t ) y_2(t) y2(t)
  7. 绘制结果:使用 matplotlib 绘制 y 1 ( t ) y_1(t) y1(t) y 2 ( t ) y_2(t) y2(t) 随时间变化的曲线。

通过这些步骤,我们成功地使用 odeint 求解了一个高阶微分方程,并将结果进行了可视化。这个方法可以扩展到更高阶的微分方程,只需将其转换为对应的一阶微分方程组即可。

2.2.2.3 使用 solve_ivp 求解一阶微分方程

Solve_ivp函数的用法与odeint非常类似,只不过比odeint多了两个参数。一个是t_span参数,表示自变量的取值范围;另一个是method参数,可以选择多种不同的数值求解算法。常见的内置方法包括RK45, RK23, DOP853, Radau, BDF等多种方法,通常使用RK45多一些。
示例:使用 solve_ivp 求解一阶微分方程

考虑相同的一阶微分方程:
d y d x = − 2 y + 1 \frac{dy}{dx} = -2y + 1 dxdy=2y+1
初始条件 y ( 0 ) = 0 y(0) = 0 y(0)=0

使用 solve_ivp 求解:

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 定义函数
def model(x, y):
    dydx = -2*y + 1
    return dydx

# 初始条件
y0 = [0]

# x 的取值范围
x_span = (0, 5)
x = np.linspace(0, 5, 100)

# 求解微分方程
sol = solve_ivp(model, x_span, y0, t_eval=x)

# 绘制结果
plt.plot(sol.t, sol.y[0])
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of dy/dx = -2y + 1')
plt.show()

解释

  1. 导入必要的库:numpy, solve_ivpmatplotlib
  2. 定义微分方程的函数 model,返回导数 dydx
  3. 设置初始条件 y0 = [0]
  4. 定义 x 的取值范围 x_span 和用于评估的点 x
  5. 使用 solve_ivp 函数求解微分方程。
  6. 使用 matplotlib 绘制解的曲线。

由于没有设置method参数,这里默认solve_ivp使用RK45(4-5阶 Runge-Kutta 法)方法进行求解。

以上示例展示了如何使用 SymPy 求解微分方程的解析解,以及如何使用 SciPy 求解微分方程的数值解。通过这些方法,我们可以处理各种类型的微分方程,满足不同应用场景的需求。

2.2.2.4 odeintsolve_ivp方法对比

上面学习了如何调用odeint和solve_ivp解微分方程,接下来对两种方法进行对比总结:

  • odeint分析

odeintSciPy 库中经典的 ODE 求解器,基于 LSODA 算法,能够自动在非刚性和刚性方法之间切换。以下是一些特点:

优点:

  1. 简单易用:适合快速求解常见的 ODE 问题。
  2. 自动方法选择:根据问题的性质自动选择合适的算法(非刚性或刚性)。
  3. 广泛使用:在很多早期的应用和文献中使用广泛,社区支持强大。

缺点:

  1. 有限的功能:无法处理事件(events)或不连续点,不支持更多高级功能。
  2. 未来维护不确定:虽然目前仍在使用,但未来可能逐渐被 solve_ivp 所取代。
  • solve_ivp分析

solve_ivpSciPy 库中较新的 ODE 求解器,提供了更多的灵活性和功能,支持多种算法。以下是一些特点:

优点:

  1. 多种算法选择:支持多种求解算法(如 RK45, RK23, BDF, LSODA 等),用户可以根据问题选择最合适的算法。
  2. 事件处理:支持事件(events)功能,可以处理不连续点和终止条件等高级功能。
  3. 灵活性高:提供更多配置选项,适合复杂问题的求解。

缺点:

  1. 稍微复杂:相对于 odeint,使用稍微复杂一些,特别是涉及到高级功能时。
  2. 较新:虽然功能更强大,但一些传统用户可能更习惯于使用 odeint
  • 对比总结
特点odeintsolve_ivp
易用性较简单相对复杂
算法选择自动选择非刚性或刚性算法用户可选择多种算法
高级功能不支持事件处理支持事件处理、更多配置选项
维护和未来前景经典方法,可能逐渐被 solve_ivp 取代较新方法,未来维护和发展前景更好
社区支持使用广泛,社区支持强大支持不断增加,功能更强大

2.3 偏微分方程的数值求解

有限差分法是一种数值方法,通过使用离散点上的差分来逼近微分,进而求解偏微分方程。

以下是使用有限差分法求解偏微分方程的一般步骤:

2.3.1 离散化区域

将求解区域离散化,即将连续的空间和时间划分为有限个离散点。例如,在二维空间上,网格点可以表示为 ( x i , y j ) (x_i, y_j) (xi,yj),时间点可以表示为 t k t_k tk

示例
对于一个在二维空间和时间上的偏微分方程,我们可以将区域划分为一个 m × n m \times n m×n 的网格:

x i = x min + i Δ x , y j = y min + j Δ y , t k = t min + k Δ t x_i = x_{\text{min}} + i \Delta x, \quad y_j = y_{\text{min}} + j \Delta y, \quad t_k = t_{\text{min}} + k \Delta t xi=xmin+iΔx,yj=ymin+jΔy,tk=tmin+kΔt

其中, Δ x \Delta x Δx Δ y \Delta y Δy Δ t \Delta t Δt 分别是空间和时间步长,(i, j, k) 是整数索引。

2.3.2. 差分公式

将偏微分方程中的微分项用差分公式代替。常见的差分公式包括前向差分、后向差分和中心差分。

示例:
对于一维热传导方程:

∂ u ∂ t = α ∂ 2 u ∂ x 2 \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} tu=αx22u

可以使用前向差分来逼近时间导数,中心差分来逼近空间导数:

u i k + 1 − u i k Δ t = α u i + 1 k − 2 u i k + u i − 1 k ( Δ x ) 2 \frac{u_i^{k+1} - u_i^k}{\Delta t} = \alpha \frac{u_{i+1}^k - 2u_i^k + u_{i-1}^k}{(\Delta x)^2} Δtuik+1uik=α(Δx)2ui+1k2uik+ui1k

2.3.3 边界条件和初始条件

在离散化网格上设置边界条件和初始条件。

示例:
对于热传导方程,边界条件可以是定值边界条件(Dirichlet)或绝热边界条件(Neumann)。初始条件是 t = 0 t = 0 t=0 时刻的温度分布:

u i 0 = f ( x i ) u_i^0 = f(x_i) ui0=f(xi)

2.3.4 构造代数方程组

根据差分公式,将偏微分方程转化为代数方程组。在每个离散点上,写出相应的代数方程。

2.3.5 求解代数方程组

使用数值方法(如高斯消元法、LU分解法或迭代法)求解代数方程组,得到离散点上的解。

2.3.6 后处理

对得到的离散点上的数值解进行后处理,如插值、绘图、误差分析等。

示例代码:一维热传导方程

以下是使用有限差分法求解一维热传导方程的Python代码示例:

import numpy as np
import matplotlib.pyplot as plt

# 参数设置
L = 1.0      # 杆的长度
T = 0.1      # 总时间
alpha = 0.01 # 热传导系数
Nx = 50      # 空间离散点数
Nt = 1000    # 时间离散点数

dx = L / (Nx - 1)
dt = T / Nt
x = np.linspace(0, L, Nx)
t = np.linspace(0, T, Nt)

# 初始条件 u(x,0) = sin(pi*x)
u = np.sin(np.pi * x)

# 时间步进
for k in range(Nt):
    u_new = u.copy()
    for i in range(1, Nx - 1):
        u_new[i] = u[i] + alpha * dt / dx**2 * (u[i+1] - 2*u[i] + u[i-1])
    u = u_new

# 绘图
plt.plot(x, u, label='Numerical Solution')
plt.xlabel('x')
plt.ylabel('u')
plt.title('1D Heat Conduction')
plt.legend()
plt.show()

4.微分方程案例

5. 差分方程案例

6. 元胞自动机

6.1 元胞自动机介绍

元胞自动机(Cellular Automaton,CA)是一种离散的数学模型,由具有有限状态的元胞和规则组成,用于模拟复杂系统中的局部互动和全局演化。元胞自动机的基本结构包括:

  1. 网格:通常是二维的矩形网格,但也可以是其他形状。
  2. 元胞:每个网格点代表一个元胞,具有有限的状态集。
  3. 邻居:每个元胞有一个邻居集,定义了元胞之间的互动范围。常见的邻居类型包括Moore邻居(八个方向)和Von Neumann邻居(四个方向)。
  4. 规则:元胞状态根据特定的局部规则进行更新,规则通常依赖于元胞当前状态及其邻居的状态。

一个经典的元胞自动机示例是康威的“生命游戏”(Conway’s Game of Life),其中元胞的状态是生或死,通过简单的规则模拟生物群体的演化。

6.2 元胞自动机的应用场景

元胞自动机因其简单的结构和强大的模拟能力,在多个领域得到了广泛应用:

  1. 生物学:用于模拟细胞的生长和分裂、生物群体的演化等。
  2. 物理学:模拟晶体生长、流体动力学、扩散过程等。
  3. 计算科学:用于研究计算过程、开发新型计算模型,如量子计算。
  4. 社会科学:模拟城市发展、交通流量、疫情传播等社会现象。
  5. 生态学:用于模拟生态系统的动态变化和物种之间的互动。
  6. 图像处理:应用于图像分割、边缘检测等领域。

6.3 元胞自动机的优点

元胞自动机具有以下优点:

  1. 简单性:模型结构简单,规则易于理解和实现。
  2. 并行性:元胞自动机的计算可以天然地并行化,适合高性能计算和分布式计算。
  3. 局部性:规则基于局部信息,不需要全局知识,适用于大规模系统的模拟。
  4. 灵活性:可以通过改变规则和邻居结构来适应不同的应用场景。
  5. 自组织行为:能够模拟复杂系统中的自组织现象,如模式形成、集群行为等。

示例:康威的“生命游戏”

“生命游戏”是最著名的元胞自动机之一,用于展示元胞自动机的基本原理和应用,生命游戏(Game of Life)是由英国数学家约翰·康威(John Horton Conway)于1970年发明的一种元胞自动机。它是一个零玩家游戏,意味着游戏的演化是由初始状态决定的,不需要玩家的进一步输入。生命游戏的规则非常简单,但它能够产生极其复杂的行为,被誉为“元胞自动机的代表作”。生命游戏的规则如下:

  • (Exposure)当前细胞为存活状态时,当周围的存活细胞低于2个时(不包含2个),该细胞变成死亡状态
  • (Survive)当前细胞为存活状态时,当周围有2个或3个存活细胞时,该细胞保持原样
  • (Overcrowding)当前细胞为存活状态时,当周围有超过3个存活细胞时,该细胞变成死亡状态
  • (Reproduction)当前细胞为死亡状态时,当周围有3个存活细胞时,该细胞变成存活状态

Python 实现:

import numpy as np
import matplotlib.pyplot as plt

def initialize_grid(size):
    """初始化网格,随机填充0和1"""
    return np.random.choice([0, 1], size*size).reshape(size, size)

def count_neighbors(grid, x, y):
    """计算元胞(x, y)的活邻居数"""
    return np.sum(grid[max(0, x-1):min(grid.shape[0], x+2), max(0, y-1):min(grid.shape[1], y+2)]) - grid[x, y]

def update_grid(grid):
    """根据生命游戏规则更新网格"""
    new_grid = grid.copy()
    for i in range(grid.shape[0]):
        for j in range(grid.shape[1]):
            alive_neighbors = count_neighbors(grid, i, j)
            if grid[i, j] == 1:
                # Exposure: 少于2个邻居,细胞死亡
                if alive_neighbors < 2:
                    new_grid[i, j] = 0
                # Survive: 2或3个邻居,细胞保持存活
                elif alive_neighbors in [2, 3]:
                    new_grid[i, j] = 1
                # Overcrowding: 超过3个邻居,细胞死亡
                else:
                    new_grid[i, j] = 0
            else:
                # Reproduction: 恰好3个邻居,细胞复活
                if alive_neighbors == 3:
                    new_grid[i, j] = 1
    return new_grid

# 设置网格大小和步数
size = 50
steps = 100

# 初始化网格
grid = initialize_grid(size)

# 动画展示
fig, ax = plt.subplots()
img = ax.imshow(grid, interpolation='nearest', cmap='binary')

def animate(frame):
    global grid
    grid = update_grid(grid)
    img.set_data(grid)
    return img,

ani = plt.FuncAnimation(fig, animate, frames=steps, interval=100, blit=True)
plt.show()

7. 数值计算方法与微分方程求解

在Python中,数值计算方法主要依赖于一些专门的库,如NumPy、SciPy和SymPy等。这些库提供了丰富的数学函数和算法,用于处理线性代数、微积分、优化问题等。例如,SciPy中的scipy.optimize模块可以用于求解函数的极值,scipy.integrate模块可以用于数值积分,而scipy.linalg模块则提供了线性代数的相关功能。通过这些工具,我们可以在Python中有效地进行数值计算和求解微分方程。
程序库函数求解极大缩短我们的计算耗时,使得模型求解不再头疼,但我们不应只成为调包侠,还应了解程序的底层数值计算方法,下面是我们日常建模中常用到的微分方程求解方法。

7.1 梯度下降法

梯度下降法是一种用于寻找函数最小值的优化算法,尤其在机器学习中广泛应用。它的基本思想是沿着函数梯度的反方向迭代,逐步逼近最小值。

基本步骤

  1. 初始化参数 θ \theta θ
  2. 计算目标函数 J ( θ ) J(\theta) J(θ) 对参数的梯度 ∇ J ( θ ) \nabla J(\theta) J(θ)
  3. 更新参数: θ = θ − α ∇ J ( θ ) \theta = \theta - \alpha \nabla J(\theta) θ=θαJ(θ),其中 α \alpha α 是学习率。
  4. 重复步骤 2 和 3 直到收敛。

公式
θ t + 1 = θ t − α ∇ J ( θ t ) \theta_{t+1} = \theta_t - \alpha \nabla J(\theta_t) θt+1=θtαJ(θt)

用梯度下降发求解二次函数获取极值是的解
目标函数是一个简单的二次函数:

J ( θ ) = ( θ − 3 ) 2 J(\theta) = (\theta - 3)^2 J(θ)=(θ3)2

这是一个抛物线形函数,其获取极小值时的解为 θ = 3 \theta = 3 θ=3
示例代码

import numpy as np
import matplotlib.pyplot as plt

# 定义目标函数及其梯度
def J(theta):
    return (theta - 3)**2

def gradient_J(theta):
    return 2 * (theta - 3)

# 梯度下降法
def gradient_descent(initial_theta, learning_rate, iterations):
    theta = initial_theta
    history = []  # 用于保存每次迭代的参数和目标函数值
    for i in range(iterations):
        grad = gradient_J(theta)
        theta = theta - learning_rate * grad
        history.append((theta, J(theta)))
        print(f'Iteration {i+1}: theta = {theta}, J(theta) = {J(theta)}')
    return theta, history

# 参数设置
initial_theta = 0.0
learning_rate = 0.1
iterations = 100

# 求解
optimal_theta, history = gradient_descent(initial_theta, learning_rate, iterations)
print(f'Optimal theta: {optimal_theta}')

# 绘制收敛过程
thetas, costs = zip(*history)
plt.plot(thetas, costs, '-o')
plt.xlabel('Theta')
plt.ylabel('J(Theta)')
plt.title('Convergence of Gradient Descent')
plt.show()

在这里插入图片描述
从图中可以看出, T h e t a Theta Theta随着步长的增加, J ( T h e t a ) J(Theta) J(Theta)不断接近最小值0。
循环100步时候, T h e t a Theta Theta=2.9999999993888893

7.2 牛顿法

牛顿法是一种用于求解非线性方程的迭代方法,通过泰勒展开近似来寻找函数的根,具有较快的收敛速度。

基本步骤

  1. 初始化值 x 0 x_0 x0
  2. 计算函数 f ( x ) f(x) f(x) 和导数 f ′ ( x ) f'(x) f(x)
  3. 更新值: x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} xn+1=xnf(xn)f(xn)
  4. 重复步骤 2 和 3 直到收敛。

公式
x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} xn+1=xnf(xn)f(xn)

示例代码

import numpy as np

# 定义函数及其导数
def f(x):
    return x**2 - 2

def df(x):
    return 2 * x

# 牛顿法
def newton_method(initial_x, iterations):
    x = initial_x
    for _ in range(iterations):
        x = x - f(x) / df(x)
    return x

# 参数设置
initial_x = 1.0
iterations = 10

# 求解
root = newton_method(initial_x, iterations)
print(f'Root: {root}')

7.3 Euler 法与 Runge-Kutta 法

Euler 法和 Runge-Kutta 法是数值求解常微分方程的基本方法。Euler 法简单但误差较大,而 Runge-Kutta 法通过更高阶的近似提高了精度。

7.3.1 Euler 法

公式
y n + 1 = y n + h f ( t n , y n ) y_{n+1} = y_n + h f(t_n, y_n) yn+1=yn+hf(tn,yn)

示例代码

import numpy as np
import matplotlib.pyplot as plt

# 定义微分方程
def f(t, y):
    return t - y

# Euler 法
def euler_method(f, t0, y0, h, n):
    t = np.zeros(n)
    y = np.zeros(n)
    t[0], y[0] = t0, y0
    for i in range(1, n):
        t[i] = t[i-1] + h
        y[i] = y[i-1] + h * f(t[i-1], y[i-1])
    return t, y

# 参数设置
t0, y0 = 0, 1
h = 0.1
n = 100

# 求解
t, y = euler_method(f, t0, y0, h, n)

# 绘制结果
plt.plot(t, y, label='Euler method')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.show()
7.3.2 Runge-Kutta 法

最常用的是四阶 Runge-Kutta 方法(RK4),具有较高的精度。

公式
k 1 = h f ( t n , y n ) k 2 = h f ( t n + h 2 , y n + k 1 2 ) k 3 = h f ( t n + h 2 , y n + k 2 2 ) k 4 = h f ( t n + h , y n + k 3 ) y n + 1 = y n + 1 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) \begin{aligned} k_1 &= h f(t_n, y_n) \\ k_2 &= h f(t_n + \frac{h}{2}, y_n + \frac{k_1}{2}) \\ k_3 &= h f(t_n + \frac{h}{2}, y_n + \frac{k_2}{2}) \\ k_4 &= h f(t_n + h, y_n + k_3) \\ y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned} k1k2k3k4yn+1=hf(tn,yn)=hf(tn+2h,yn+2k1)=hf(tn+2h,yn+2k2)=hf(tn+h,yn+k3)=yn+61(k1+2k2+2k3+k4)

示例代码

import numpy as np
import matplotlib.pyplot as plt

# 定义微分方程
def f(t, y):
    return t - y

# Runge-Kutta 法
def runge_kutta_method(f, t0, y0, h, n):
    t = np.zeros(n)
    y = np.zeros(n)
    t[0], y[0] = t0, y0
    for i in range(1, n):
        t[i] = t[i-1] + h
        k1 = h * f(t[i-1], y[i-1])
        k2 = h * f(t[i-1] + h/2, y[i-1] + k1/2)
        k3 = h * f(t[i-1] + h/2, y[i-1] + k2/2)
        k4 = h * f(t[i-1] + h, y[i-1] + k3)
        y[i] = y[i-1] + (k1 + 2*k2 + 2*k3 + k4) / 6
    return t, y

# 参数设置
t0, y0 = 0, 1
h = 0.1
n = 100

# 求解
t, y = runge_kutta_method(f, t0, y0, h, n)

# 绘制结果
plt.plot(t, y, label='Runge-Kutta method')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.show()

7.4 Crank-Nicolson 法

Crank-Nicolson 法是一种用于求解偏微分方程(PDE)的隐式数值方法,具有较高的精度和稳定性。它是 Euler 法和中心差分法的组合,适用于时间和空间的二阶精度计算。

公式
y n + 1 − y n h = 1 2 ( f ( t n , y n ) + f ( t n + 1 , y n + 1 ) ) \frac{y_{n+1} - y_n}{h} = \frac{1}{2} \left( f(t_n, y_n) + f(t_{n+1}, y_{n+1}) \right) hyn+1yn=21(f(tn,yn)+f(tn+1,yn+1))

示例代码

import numpy as np
import matplotlib.pyplot as plt

# 定义微分方程
def f(t, y):
    return t - y

# Crank-Nicolson 法
def crank_nicolson_method(f, t0, y0, h, n):
    t = np.zeros(n)
    y = np.zeros(n)
    t[0], y[0] = t0, y0
    for i in range(1, n):
        t[i] = t[i-1] + h
        y_pred = y[i-1] + h * f(t[i-1], y[i-1])  # 预测
        y[i] = y[i-1] + (h / 2) * (f(t[i-1], y[i-1]) + f(t[i], y_pred))  # 校正
    return t, y

# 参数设置
t0, y0 = 0, 1
h = 0.1
n = 100

# 求解
t, y = crank_nicolson_method(f, t0, y0, h, n)

# 绘制结果
plt.plot(t, y, label='Crank-Nicolson method')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.show()

8.学习心得

本次学习了微分方程的基础概念,了解到公式推导和python编程求解两种方法求解微分方程,通过应用案例掌握了微分方程和差分方程题目的解题思路和编程技巧。

微分方程的求解有两种情况:求解解析解或者给定特定条件求解精确解(部分微分方程是没有解析解的)

  • 通过SymPy 求解微分方程的解析式解;
  • 通过Scipy的odeintsolve_ivp求精确解,其中solve_ivpSciPy 库中较新的 ODE 求解器,提供了更多的灵活性和功能,支持多种算法。如果遇到高阶微分方程,要通过变元替换高阶项,使方程变成一阶方程再进行代码求解。

对于偏微分方程,python没有提供专门的解析包,需要更加实际情况编写代码求解。核心思想是将连续离散化处理
为解决差分方程模型往往难以有效地模拟大多数离散模型问题,进而引入了元胞自动机模型

9. 对学习文档的一些疑问

9.1 文档2.1.4中梯型求解代码是否有误

在这里插入图片描述
此处的梯型法求解的代码有误,for循环中 X n Xn Xn初值应当为0,而不是0.7,以下是修改后的代码,仅供参考

def tixing():
    # 积分区间
    x0 = 0
    xn = 0.7

    # 步长
    n = 1000  # 分成1000个小区间
    h = (xn - x0) / n
    s = 0

    xn = 0 #此处先令xn为0,随循环逐渐递进
    for i in range(n):
        xn1 = xn + h;
        yn = f(xn)
        yn1 = f(xn1)
        s0 = (yn + yn1) * h / 2
        s += s0
        xn = xn1
    print(s)

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/748503.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

记录一个80端口被占用导致软件打不开的问题

今天有个电脑&#xff0c;安装完我们的软件后&#xff0c;在浏览器上面打不开。但是我看虚拟机里面的配置啥的都很正常&#xff0c;我感觉不是软件挂了&#xff0c;应该是系统哪里的配置出了问题&#xff0c;导致软件打不开。 跟做软件的联系了&#xff0c;他让我直接访问虚拟机…

酷开系统丨开启家庭智能教育让学习成为一种乐趣

在数字化时代&#xff0c;孩子们接触的信息日益增多&#xff0c;而酷开系统洞察到了家长对孩子成长环境的关切。酷开系统&#xff0c;作为家庭娱乐与教育的融合平台&#xff0c;不仅注重提供丰富的教育资源&#xff0c;更致力于创造一个健康、有益的学习和娱乐环境。 在酷开系…

OPenFast软件中的NRELOffshrBsline5MW_Onshore_ServoDyn.dat文件详解

我先简单概括一下&#xff0c;后续我再详细总结&#xff1a;文件“NRELOffshrBsline5MW_Onshore_ServoDyn.dat”是用于NREL 5.0 MW基准风力发电机的ServoDyn模块的输入文件。它定义了仿真控制、变桨控制、发电机和扭矩控制、偏航控制以及输出设置等各种参数。以下是主要内容的总…

Uniapp的使用

为什么要使用uniapp uniapp 可以进行多端开发&#xff0c;uniapp 在设计的时候就拥有许多兼容性代码&#xff0c;可以兼容很多的平台 如 支付宝小程序 html页面 微信小程序等&#xff0c;注重开发效率而不是运行效率时 &#xff0c;就可以考虑一下 uniapp 当然也可以去…

mysql窗口函数选择每个窗口的第一条数据

需求 假设我们有一个名为sales的表&#xff0c;我们想要按product分组&#xff0c;并为每个产品选择销售额最高的那一天 sales表 首先给每个产品按照销售量进行排名 SELECT product,sale_date,sales_amount,ROW_NUMBER() OVER (PARTITION BY product ORDER BY sales_amount …

Langchain-Chatchat 搭建知识库教程,老奶奶都能看懂的教程

本文将在 google 实验室中使用 Langchain-Chatchat 搭建一个知识库&#xff0c;还可以进行聊天等功能。 由于是在 google 实验室上面跑代码&#xff0c;所以本地电脑什么配置都无所谓&#xff01; 效果图 运行起来后可以上传各种文档文件到知识库。 Langchain-Chatchat 是什…

sheng的学习笔记-hadoop,MapReduce,yarn,hdfs框架原理

目录 搭建hadoop 下载hadoop JAVA 下载bin windows 改环境变量 将winutils.exe和hadoop.dll放到C:\Windows\System32下&#xff0c;然后重启 修改配置 vim core-site.xml vim hdfs-site.xml hadoop-env.sh mapred-site.xml yarn-site.xml 格式化命令 启动集群 …

浅谈 MySQL 复制架构

Author&#xff1a;Arsen Date&#xff1a;2024/06/26 目录 前言一、参数设置1.1 slave_exec_mode1.2 max_allowed_packet1.3 binlog-do-db1.4 binlog-ignore-db1.5 replicate-ignore-db1.6 replicate-ignore-table1.7 replicate-wild-ignore-table1.8 slave_compressed_protoc…

【React】变量 useState

开发需要&#xff0c;随便学学react。上手第一天&#xff0c;感觉这个JS语法很怪&#xff0c;没有什么逻辑性&#xff0c;比较抽象。随便写写笔记。 跟着网上找的项目写写感觉这个项目还不错&#xff1a; 分享给码友 https://zh-hans.react.dev/learn/tutorial-tic-tac-toe 参…

最新MDYS14源码影视视频网站模板/苹果CMS系统/附搭建教程

最新MDYS14源码影视视频网站模板/苹果CMS系统/附搭建教程 基本介绍&#xff1a; 1、后台增加自定义参数&#xff0c;对应会员升级页面&#xff0c;以及积分充值 2、视频&#xff0c;演员&#xff0c;专题&#xff0c;收藏&#xff0c;会员系统模块齐全&#xff0c;支持子分类…

Taro + vue3 中微信小程序中实现拉起支付

在前端开发中 H5 的拉起支付和微信小程序的拉起支付 是不太一样的 现在分享一下微信小程序中的拉起支付 逻辑都在后端 我是用的Taro 框架 其实就是一个Api Taro 文档

【vue3】【vant】 移动端古诗词句子发布收藏项目

更多项目点击&#x1f446;&#x1f446;&#x1f446;完整项目成品专栏 【vue3】【vant】 移动端古诗词句子发布收藏项目 获取源码方式项目说明&#xff1a;其中功能包括素材包含&#xff1a;项目运行环境运行截图 获取源码方式 加Q群&#xff1a;632562109项目说明&#xf…

2024智能驾驶兴趣人群研究报告

来源&#xff1a;百分点舆情中心 近期历史回顾&#xff1a; 劳动力效能提升指引白皮书》人效研究院.pdf 【标准】企业ESG管理体系(T-CERDS 5—2023).pdf 【实用标准】GB_T 43868-2024 电化学储能电站启动验收规程.pdf 【实用模板】用户侧新型储能项目管理流程图及备案资料清单…

二叉树——另一颗树的子树

目录 1&#xff1a;题目分析及思路 2&#xff1a;代码实现和分析 1&#xff1a;代码 2&#xff1a;分析 1&#xff1a;题目分析及思路 给我们两棵二叉树&#xff0c;分别是 root 和 subRoot 。检验 root 中是否包含和 subRoot 具有相同结构和节点值的子树。如果存在&…

Planned independent reguirement can only be maintained via the network

背景&#xff1a;用户上线ps系统&#xff0c;物料用策略70跑需求 但是因为通用料被改了策略&#xff0c;改成其他的了&#xff0c;影响到计划独立需求了。 如果用户不需要了哪个料就会把数量改为0&#xff0c;或者直接删掉物料。之前建议是改成0&#xff0c;这样还有个记录在…

Pandas中的数据转换[细节]

今天我们看一下Pandas中的数据转换&#xff0c;话不多说直接开始&#x1f387; 目录 一、⭐️apply函数应用 apply是一个自由度很高的函数 对于Series&#xff0c;它可以迭代每一列的值操作&#xff1a; 二、⭐️矢量化字符串 为什么要用str属性 替换和分割 提取子串 …

如何ubuntu安装wine/deep-wine运行exe程序(包括安装QQ/微信/钉钉)

1.失败的方法&#xff1a; ubuntu22.04尝试下面这个链接方法没有成功&#xff0c; ubuntu22.04安装wine9.0_ubuntu 22.04 wine-CSDN博客 上面链接里面也提供了wine官方方法&#xff0c;链接如下&#xff1a;https://wiki.winehq.org/Ubuntu_zhcn 但是运行最后一个命令时候报…

HTTP-02

常用HTTP状态码是怎么分类的 常用的HTTP状态码是按照以下几个分类进行的&#xff1a; 1xx 信息类状态码&#xff1a;表示请求已被接收&#xff0c;需要进一步处理。 2xx 成功类状态码&#xff1a;表示请求已成功被服务器接收、理解和处理。 3xx 重定向类状态码&#xff1a;表…

你的生日是星期几?HTML+JavaScript帮你列出来

0 源起 上周末&#xff0c;大宝发现今年自己的生日不是周末&#xff0c;这样就不好约同学和好友一起开生日Party了&#xff0c;很是郁闷。一直嘀咕自己哪年的生日才是周末。 于是我用JavaScript写了一个小程序来帮她测算了未来100年中每年的生日分别是星期几。 1 设计交互界面…

RabbitMQ的Direct交换机

Direct交换机 BindingKey 在Fanout模式中&#xff0c;一条消息&#xff0c;会被所有订阅的队列都消费。但是&#xff0c;在某些场景下&#xff0c;我们希望不同的消息被不同的队列消费。这时就要用到Direct类型的Exchange。 在Direct模型下&#xff1a; 队列与交换机的绑定&a…