Python解微分方程组

发表于 更新于

近期在解决一组微分方程的问题，所以需要使用Python 编程。首先介绍一下编程的相同原理：

任何程序都有一些已经编译好的，完成通用功能的部分，也就是说当你解决问题时要实现某些功能，这些功能已经内置于现有语言中，不需要作者再独立编写。比如在LaTeX中在导言引入宏包, 这些宏包定义了一些命令，不需要作者再自己定义了。Python中引入一些库, 这里的库和LaTeX中的宏包本质上是同一个东西，就是一些已经实现了的功能模块，将其加入到所写的程序中即可。对于高级玩家，还可以添加自己定义的一些特殊需求的模块，几乎所有的语言都是这个模式。

对于Python， 其常用的一些库，可以参考文章： 38个常用Python库：数值计算、可视化、机器学习等8大领域都有了, 直接上我写的一组代码：

Time Cristal
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
#! /usr/bin/env python3
# vim:fenc=utf-8
#
# Copyright © 2023 feng <feng@archlinux>
#
# Distributed under terms of the MIT license.

import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from sympy.abc import t

def fun_new(y, t, a, b):
    f, df, g, dg = y
    # my formular
    dydt = [ df , (-a+b)*f- t*g , dg, (-a-b)*g-t*f ]
    # article formular
    # dydt = [ df , (-a+t)*f+b*g , dg, (-a-t)*g+b*f ]
    return dydt

# 初始条件
y0=[0,-0.36012, 1, 0]
# t=np.linspace(0,5,101)
t_left=np.arange(0,-6,-0.001)
sol_left = odeint(fun_new , y0, t_left, args=(5/(2**(1/3)),2**(2/3)))
t_right=np.arange(0,6,0.001)
sol_right= odeint(fun_new , y0, t_right, args=(5/(2**(1/3)),2**(2/3)))

plt.plot(t_left,sol_left[:,0], color="r", label="$\phi_1$")
plt.plot(t_left,sol_left[:,2], color="b", label="$\phi_2$")
plt.plot(t_right,sol_right[:,0], color="r")
plt.plot(t_right,sol_right[:,2], color="b")
plt.axhline(0,color="gray") # 画过原点的水平线
plt.axvline(0,color="gray") # 画过原点的竖直线
plt.text(0.5, 1.5, r"$(-\frac{d^2}{dy^2}-\frac{\epsilon}{2^{1/3}}+2^{2/3})\phi_1(y)-y\phi_2(y)=0$", fontsize=20)
plt.text(0.5, 1.1, r"$(-\frac{d^2}{dy^2}-\frac{\epsilon}{2^{1/3}}-2^{2/3})\phi_2(y)-y\phi_1(y)=0$", fontsize=20)
plt.title("Zhen-Hua Feng FIG.11 $\epsilon=5$ $\phi_1=0$, $\phi_1'=-0.36012$, $\phi_2=1$, $\phi_2'=0$ ", fontproperties="SimHei", fontsize=20)
plt.xlabel("y")
plt.ylabel("$\phi_1$ and $\phi_2$")
plt.legend()
plt.show()

这里需要解释的是关于高阶微分的定义，代码如下

高阶微分的定义
1
2
3
4
5
def fun_new(y, t, a, b):
    f, df, g, dg = y
    # my formular
    dydt = [ df , (-a+b)*f- t*g , dg, (-a-b)*g-t*f ]
    return dydt
  • 第2行，f是一个函数，df是函数的一次微分， g 是另一个函数， dg 是函数的另一个微分
  • 第4行，使用一个数组存存数据，第1列为, df 表示函数 f 的微分，第2列是 df 的微分（也就是 f 的二阶导数）， 第3列是 g 的微分， 第4列是 dg 的微分（也就是 g 的二阶导数）
  • 第5行，一个 1X 4 数组，为plt绘图提供数据。
根据初值生成数据
1
2
3
y0=[0,-0.36012, 1, 0]
t_left=np.arange(0,-6,-0.001)
sol_left = odeint(fun_new , y0, t_left, args=(5/(2**(1/3)),2**(2/3)))
  • 第1行，定义初值，当然这里的初值应当根据物理情景提出。
  • 第2行，使用arange命令，生成自变量，步长取为-0,001, 这个会向前取值，同理步长取为正，则会向后生成自变量离散值。
  • 第3行， 使用odeint 代入自变量及参数后生成由fun_new构成的数组，这个数组和自变量的离散数一样多，然后再使用plt据这些数据绘图。