标签:python numpy numerical-methods differential-equations pde
我正在尝试编写一个python程序,使用具有二阶空间离散和周期边界条件的显式欧拉方法来求解一阶一维波动方程(传输方程).
我是python的新手,我使用numpy编写了这个程序,但我觉得我在某个地方犯了一个错误,因为波被扭曲了.一旦它离开左边界,它似乎变得扭曲,而不是简单地向左转换.我很确定这是一个编程错误,但有可能是舍入错误吗?我没有正确使用numpy吗?有关以更多Python风格的方式编写此程序的任何建议吗?谢谢!
PDE是
它是有限差分形式
解决
这是我尝试的:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# wave speed
c = 1
# spatial domain
xmin = 0
xmax = 1
n = 50 # num of grid points
# x grid of n points
X, dx = np.linspace(xmin,xmax,n,retstep=True)
# for CFL of 0.1
dt = 0.1*dx/c
# initial conditions
def initial_u(x):
return np.exp(-0.5*np.power(((x-0.5)/0.08), 2))
# each value of the U array contains the solution for all x values at each timestep
U = []
# explicit euler solution
def u(x, t):
if t == 0: # initial condition
return initial_u(x)
uvals = [] # u values for this time step
for j in range(len(x)):
if j == 0: # left boundary
uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][j+1]-U[t-1][n-1]))
elif j == n-1: # right boundary
uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][0]-U[t-1][j-1]))
else:
uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][j+1]-U[t-1][j-1]))
return uvals
# solve for 500 time steps
for t in range(500):
U.append(u(X, t))
# plot solution
plt.style.use('dark_background')
fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)
# animate the time data
k = 0
def animate(i):
global k
x = U[k]
k += 1
ax1.clear()
plt.plot(X,x,color='cyan')
plt.grid(True)
plt.ylim([-2,2])
plt.xlim([0,1])
anim = animation.FuncAnimation(fig,animate,frames=360,interval=20)
plt.show()
任何人都可以解释为什么这(波浪失真)正在发生?
解决方法:
您的实施是正确的.失真来自相对较大的空间步长dx.在其当前值0.2处,它与波的大小相当,这使得波在图上可见多边形.这些离散误差累积超过500步.这是我从plt.plot(X,U [-1])得到的代码:
这是我使用n = 100(时间和空间步长减半),运行范围(1000)中的t的解决方案以补偿较小的时间步长,并再次绘制plt.plot(X,U [-1 ]):
du / dx的对称差分近似具有与三阶导数成比例的阶数dx ** 3的误差.这些累积的方式很复杂,因为解决方案正在四处移动,但无论如何,如果dt与它成比例,则较小的dx会改善这些问题.
标签:python,numpy,numerical-methods,differential-equations,pde 来源: https://codeday.me/bug/20190627/1304625.html
本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享; 2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关; 3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关; 4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除; 5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。