python求解偏微分方程_Python数值计算----------求解简单的偏微分方程
很多物理現象的都可以用方程來描述,比如熱傳導與物質擴散可以用擴散方程來描述,流體的流動可以用NS方程描述等等。如果能夠將這些偏微分方程求解出來,就可以來對很多物理現象進行仿真,現在工程中的仿真軟件都是通過程序數值求解一類偏微分方程。今天我們嘗試求解一類偏微分方程,為了簡單起見,我們以一個簡單的平流方程為例,方程形式如下:平流方程
求解偏微分方程的數值解法非常多,這里也是采用一種較為直白的方法--------有限差分法,將上述的偏微分方程離散為差分方程,這里采用一種迎風格式的差分方法,首先我們可以看出這個平流方程由兩個維度,一個空間維度x,一個時間維度t,我們對這兩個維度進行離散:差分網格
差分方程
將差分方程簡化一下
通過上述的差分方程我們可以通過上一個時間步的結果推進得到下一個時間步的結果:
這里我們邊界條件是周期邊界如下圖,周期邊界的好處就是我們可以一直看到在我們的求解域內看到波的移動:
方程的初始條件我們給空間維度上一個指數函數的波形:
下面我們通過一段簡單(其實挺長的,但是很容易理解)的Python代碼對上述差分方程進行求解:
"""This code solves the advection equationU_t + vU_x = 0over the spatial domain of 0 <= x <= 1 that is discretizedinto 103 nodes, with dx=0.01, using the First-Order Upwind (FOU)scheme with Forward difference in Eq.for an initial profile of a Gaussian curve, defined byU(x,t) = exp(-200*(x-xc-v*t).^2)where xc=0.25 is the center of the curve at t=0.The periodic boundary conditions are applied either end of the domain.The velocity is v=1. The solution is iterated until t=1.5 seconds."""
import numpy as np
import matplotlib.pyplot as plt
import time
start =time.clock()
class UpwindMethod1:
def __init__(self, N, tmax):
self.N = N # number of nodes
self.tmax = tmax
self.xmin = 0
self.xmax = 1
self.dt = 0.009 # timestep
self.v = 1 # velocity
self.xc = 0.25
self.initializeDomain()
self.initializeU()
self.initializeParams()
def initializeDomain(self):
self.dx = (self.xmax - self.xmin)/self.N
self.x = np.arange(self.xmin-self.dx, self.xmax+(2*self.dx), self.dx)
def initializeU(self):
u0 = np.exp(-200*(self.x-self.xc)**2)
self.u = u0.copy()
self.unp1 = u0.copy()
def initializeParams(self):
self.nsteps = round(self.tmax/self.dt)
def solve_and_plot(self):
tc = 0
for i in range(self.nsteps):
plt.clf()
# The FOU scheme, Eq. (18.21)
for j in range(self.N+2):
self.unp1[j] = self.u[j] - (self.v*self.dt/(self.dx))*(self.u[j]-self.u[j-1])
self.u = self.unp1.copy()
# Periodic boundary conditions
self.u[0] = self.u[self.N+1]
self.u[self.N+2] = self.u[1]
plt.plot(self.x, self.u, 'bo-', label="First-order Upwind")
plt.axis((self.xmin-0.12, self.xmax+0.12, -0.2, 1.4))
plt.grid(True)
plt.xlabel("Distance (x)")
plt.ylabel("u")
plt.legend(loc=1, fontsize=12)
plt.suptitle("Time =%1.3f" % (tc+self.dt))
plt.pause(0.01)
tc += self.dt
def main():
sim = UpwindMethod1(100, 1.5)
sim.solve_and_plot()
plt.show()
if __name__ == "__main__":
main()
end = time.clock()
print('Running time:%sSeconds'%(end-start))
運行程序大家就可以看到一個波在來回移動。這里以一個簡單平流方程為例,然后構造了其迎風格式的差分方程,然后用我們可愛的python語言進行求解,這里再次說明了計算機語言是一個工具,有了這個工具我們可以用它來做我們想做的事,對我來言Python做科學計算非常的Nice。
總結
以上是生活随笔為你收集整理的python求解偏微分方程_Python数值计算----------求解简单的偏微分方程的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 作为一个软件实施工程师的发展方向
- 下一篇: ufldl学习笔记与编程作业:Multi