很多物理現(xiàn)象的都可以用方程來(lái)描述，比如熱傳導(dǎo)與物質(zhì)擴(kuò)散可以用擴(kuò)散方程來(lái)描述，流體的流動(dòng)可以用NS方程描述等等。如果能夠?qū)⑦@些偏微分方程求解出來(lái)，就可以來(lái)對(duì)很多物理現(xiàn)象進(jìn)行仿真，現(xiàn)在工程中的仿真軟件都是通過(guò)程序數(shù)值求解一類偏微分方程。今天我們嘗試求解一類偏微分方程，為了簡(jiǎn)單起見(jiàn)，我們以一個(gè)簡(jiǎn)單的平流方程為例，方程形式如下：平流方程

求解偏微分方程的數(shù)值解法非常多，這里也是采用一種較為直白的方法--------有限差分法，將上述的偏微分方程離散為差分方程，這里采用一種迎風(fēng)格式的差分方法，首先我們可以看出這個(gè)平流方程由兩個(gè)維度，一個(gè)空間維度x，一個(gè)時(shí)間維度t，我們對(duì)這兩個(gè)維度進(jìn)行離散：差分網(wǎng)格

差分方程

將差分方程簡(jiǎn)化一下

通過(guò)上述的差分方程我們可以通過(guò)上一個(gè)時(shí)間步的結(jié)果推進(jìn)得到下一個(gè)時(shí)間步的結(jié)果：

這里我們邊界條件是周期邊界如下圖，周期邊界的好處就是我們可以一直看到在我們的求解域內(nèi)看到波的移動(dòng)：

方程的初始條件我們給空間維度上一個(gè)指數(shù)函數(shù)的波形：

下面我們通過(guò)一段簡(jiǎn)單(其實(shí)挺長(zhǎng)的，但是很容易理解)的Python代碼對(duì)上述差分方程進(jìn)行求解：

"""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))

運(yùn)行程序大家就可以看到一個(gè)波在來(lái)回移動(dòng)。這里以一個(gè)簡(jiǎn)單平流方程為例，然后構(gòu)造了其迎風(fēng)格式的差分方程，然后用我們可愛(ài)的python語(yǔ)言進(jìn)行求解，這里再次說(shuō)明了計(jì)算機(jī)語(yǔ)言是一個(gè)工具，有了這個(gè)工具我們可以用它來(lái)做我們想做的事，對(duì)我來(lái)言Python做科學(xué)計(jì)算非常的Nice。

總結(jié)

以上是生活随笔為你收集整理的python求解偏微分方程_Python数值计算----------求解简单的偏微分方程的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。