Python光學仿真數值分析求解波動方程繪制波包變化圖

Posted on by

波動方程數值解

波動方程是三大物理方程之一，也就是弦振動方程，其特點是時間與空間均為二階偏導數。其自由空間解便是我們熟知的三角函數形式，也可以寫成自然虛指數形式。

一般來說，既然有瞭精確的解析解，那也就沒必要再去做不精確的數值模擬，但數值模擬的好處有兩個，一是避免無窮小，從而在思維上更加直觀；二是頗具啟發性，對於一些解析無解的情況也有一定的處理能力。

對此，我們首先考慮一維波動方程

在這裡插入圖片描述 

import numpy as np
import matplotlib.pyplot as plt
def set_y0(x,k,L):
    y = np.zeros_like(x)
    y[x<L] = np.sin(k*x[x<L])*np.sin(np.abs(x[x<L]*np.pi/L))
    return y

if __name__ == "__main__":
    x = np.linspace(0,10,1000)
    k = np.pi*2/1.064
    L = 5
    y = set_y0(x,k,L)

    plt.plot(x,y)
    plt.show()

其形狀為

在這裡插入圖片描述

現考慮讓這個光波在 [ 0 , L ] 范圍內往返傳播，在此采用Dirichlet邊界條件，取

在這裡插入圖片描述

至此，我們得到瞭光場的所有信息，原則上可以預測這個波包的所有行為，其迭代過程為

def wave1d(x,t,k,L):
    dx = x[1]-x[0]
    dt = t[1]-t[0]
    d2 = (dt/dx)**2
    y = np.zeros([len(t),len(x)])
    y[0,:] = set_y0(x,k,L)
    y[1,:] = set_y0(x-dt,k,L)
    for n in range(2,len(t)):
        y[n] = 2*y[n-1] - y[n-2] - d2*2*y[n-1]
        y[n,1:] += d2*y[n-1,:-1]
        y[n,:-1] += d2*y[n-1,1:]
        #邊界條件
        y[n,0] = 0
        y[n,-1] = 0   
    return y

由於 y y y是隨時間變化的參量，現有的matplotlib.pyplot已經無法滿足我們繪制動態圖片的需求，所以引入animation來進行繪制，其代碼為

import matplotlib.animation as animation
#輸入時間，自變量，因變量，圖題標記
def drawGif(t,x,ys,mark="time="):
    tAxis = np.linspace(0,len(t)-1,100).astype(int)
    fig = plt.figure()
    ax = fig.add_subplot(111,xlim=(0,10),ylim=(-1.5,1.5))
    ax.grid()
    line, = ax.plot([],[],lw=0.2)
    time_text = ax.text(0.1,0.9,'',transform=ax.transAxes)
    def init():
        line.set_data([],[])
        time_text.set_text("")
        return line, time_text   
    def animate(i):
        y = ys[i]
        line.set_data(x,y)
        time_text.set_text(mark+str(t[i]))
        return line, time_text
    # 動態圖繪制命令
    # 輸入分別為畫圖窗口，動畫函數，動畫函數輸入變量，延時，初始函數
    ani = animation.FuncAnimation(fig, animate, tAxis,
        interval=200, init_func=init)
    #通過imagemagick引擎來保存gif
    ani.save('wave.gif',writer='imagemagick')
    plt.show()
if __name__ == "__main__":
    x = np.linspace(0,10,1000)
    t = np.linspace(0,12,2041)
    k = np.pi*2/1.064
    L = 5
    y = wave1d(x,t,k,L)
    drawGif(t,x,y)

得到結果為

在這裡插入圖片描述

這個圖雖然很符合我們的預期，但有些物理過程並不清晰，我們不妨把初始波包設置為隻有一個波峰的孤波

def set_y0(x,k,L):
    y = np.zeros_like(x)
    y[x<L] = np.sin(np.abs(x[x<L]*np.pi/L))
    return y

其圖像為

在這裡插入圖片描述

我們可以清晰地看到，正弦波通過腔壁後，其震動方向發生瞭變化，此即半波損失。

以上就是Python光學仿真數值分析求解波動方程繪制波包變化圖的詳細內容，更多關於Python數值分析求解波動方程繪制波包變化圖的資料請關註WalkonNet其它相關文章！

推薦閱讀：