Python 那点事--我的一个有限差分python程序

最新推荐文章于 2024-07-05 12:38:48 发布

原创最新推荐文章于 2024-07-05 12:38:48 发布 · 置顶 · 1.7k 阅读

12 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#python #算法

本文介绍了一维波动方程的有限差分法在Python中的实现过程，对比了Matlab与Python在处理大型数组方面的优劣，详细阐述了控制方程、初始条件和边界条件，并展示了使用二维列表存储数据的示例代码。

1.很多同学刚从matlab转到python时有点无所适从，特别是对于matlab的数组处理和python的列表。

matlab转python是迟早的事情，一方面是matlab处理大型数组太消耗资源，一方面是matlab的绘图板渲染效果没有python好，还有一个原因是python的库越来越强大。

我的第一个python程序是一维波动方程有限差分实现，该程序用来模拟锚杆动测问题的数学物理模型，数学物理方程表述如下：

控制方程：

$\frac{\partial^2u }{\partial t^2}=V^2_{p}\frac{\partial ^2u}{\partial x^2}$

初始条件：

$u(x,t)|_{t=0}=0\\ \frac{\partial u}{\partial t}|_{t=0}=0$

边界条件：

$\frac{\partial u}{\partial x}|_x=0=-M\\ \frac{\partial u}{\partial x}|_x=l=0$

取空间步长 $h$ 和时间步长 $\tau$ ，向前差分格式表示，令 $r=\frac{v_{p}\tau }{h}$ ,则有5点显示差分格式：

$u^{n+1}_{j}=r^{2}{u^{n}_{j-1}+u^{n}_{j+1}}+2(1-r^2)u^{n}_{j}-u^{n-1}_{j}$

u = [[0 for i in range(N)] for i in range(J)]

u[:][0] = 0
u[:][1] = 0

#瞬时激振雷克子波响应
for i in range(0, N, 1):
    if(i < 100):
        fricker.append(10000*(1 - 2 * (math.pi * fm * (i-50) * 0.001) ** 2) * math.exp(-(math.pi * fm * (i-50) * 0.001) ** 2))
    else:
        fricker.append(0)
fricker = [fricker]

for n in range(1, N-1, 1):#N是时间序列
    for j in range(1, J-1, 1):#J是空间序列
        u[j][n+1] = r**2*(u[j-1][n] + u[j+1][n]) + 2*(1-r**2)*u[j][n] - u[j][n-1]
        u[J-1][n+1] = u[J-2][n+1]
        u[0][n+1] = u[1][n+1] + h*fricker[0][n+1]

用二维列表的格式存储数组，而不一定要转换成array格式。

子波、位移、加速度场：