Python——一维波动方程求解

1、一维波动方程简介

     以线的横振动方程作为一维波动方程的例子，其运动学方程满足：

 

     在上式中为线密度；T为张力；代表强迫力的大小。当没有外力，且线的密度均匀等于，我们可以取，将一维波动方程化简为：

2、一维波动方程的差分解法

     在给定的初始条件：

 

     以及给定的边界条件：

     我们可以使用差分的方法将一维波动方程写为两种差分格式：

     第一种差分格式：

     第二种差分格式：   

   两种差分收敛条件：

 

2、问题

     用一种差分格式计算波动方程混合问题：

1. 取α=1，h=0.2，计算k=1，2，3，4层的值；
2. 取α=1，h=0.05，分别用两种差分格式计算k=1，2⋯，20层的近似值。

 3、程序

#python
import numpy as np
import math as ma
from matplotlib import pyplot as plt
#问题1算法
def n1(a,h,N,v):
    def f1(x,h,t):
        y=ma.sin(x*h*ma.pi)+x*h*t*(1-x*h)
        return y
    def f2(x,h):
        y=ma.sin(x*h*ma.pi)
        return y
    def f3(x1,x2,x3):
        y=x1+x2-x3
        return y
    t=a*h/v
    y_t=np.zeros((N-1,N))
    for i in range(N-1):
        if i==0:
            for j in range(N-1):
                y_t[i,j]=f2(j,h)
        if i==1:
            for j in range(N-1):
                y_t[i,j]=f1(j,h,t)
        if i>1:
            for j in range(1,N-1):
                y_t[i,j]=f3(y_t[i-1,j+1],y_t[i-1,j-1],y_t[i-2,j])
    return y_t
#问题2算法
def n2(a,h,N,v):
    def f1(x,h,t):
        y=1/2*(ma.sin((x+1)*h*ma.pi)+ma.sin((x-1)*h*ma.pi)+x*h*t*(1-x*h))
        return y
    def f2(x,h):
        y=ma.sin(x*h*ma.pi)
        return y
    def f3(x1,x2,x3):
        y=x1+x2-x3
        return y
    t=a*h/v;
    y_t=np.zeros((N-1,N));
    for i in range(N-1):
        if i==0:
            for j in range(N-1):
                y_t[i,j]=f2(j,h);
        if i==1:
            for j in range(1,N-1):
                y_t[i,j]=f1(j,h,t)
        if i>1:
            for j in range(1,N-2):
                y_t[i,j]=f3(y_t[i-1,j+1],y_t[i-1,j-1],y_t[i-2,j])
    return y_t
#问题1答案
a1=1;h1=0.2;N1=6;v1=1;
y1=n1(a1,h1,N1,v1)
#问题2答案
a2=1;h2=0.05;v2=1;N2=21
y2=n1(a2,h2,N2,v2)
a3=1;h3=0.05;v3=1;N3=21
y3=n2(a3,h3,N3,v3);
#结果展示
x1=[];t1=[]
for n in range(0,6):
    x1.append(n)
for m in range(0, 5):
    t1.append(m*0.05)
x1=np.array(x1)
t1=np.array(t1)
x,y=np.meshgrid(x1,t1)
z=y1
d3=plt.axes(projection='3d')
d3.plot_surface(x,y,z,cmap='rainbow')
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
d3.set_zlabel('Z(y方向位移)/m')
d3.set_xlabel('X(位置)/m')
d3.set_ylabel('T(时间)/s')
d3.set_title('一维波动方程')
plt.show()
x2=[];t2=[];t_m2=a2*h2/v2
for n in range(0,N2):
    x2.append(n)
for m in range(0,N2-1):
    t2.append(m*t_m2)
x2=np.array(x2)
t2=np.array(t2)
x,y=np.meshgrid(x2,t2)
z=y2
d3=plt.axes(projection='3d')
d3.plot_surface(x,y,z,cmap='rainbow')
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
d3.set_zlabel('Z(y方向位移)/m')
d3.set_xlabel('X(位置)/m')
d3.set_ylabel('T(时间)/s')
d3.set_title('一维波动方程')
plt.show()
x3=[];t3=[];t_m3=a3*h2/v2
for n in range(0,N3):
    x3.append(n)
for m in range(0,N3-1):
    t3.append(m*t_m3)
x3=np.array(x3)
t3=np.array(t3)
x,y=np.meshgrid(x3,t3)
z=y3
d3=plt.axes(projection='3d')
d3.plot_surface(x,y,z,cmap='rainbow')
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
d3.set_zlabel('Z(y方向位移)/m')
d3.set_xlabel('X(位置)/m')
d3.set_ylabel('T(时间)/s')
d3.set_title('一维波动方程')
plt.show()

4、结果

Fig1.问题1

Fig2.问题2-1

Fig3.问题2-2

Fig4.N=40

Fig5.N=60