有限差分波动方程

1 投票
1 回答
40 浏览
提问于 2025-04-14 17:53

我正在尝试用有限差分法来解决波动方程(之后还想把结果做成动画)。但是,似乎每个时间步得到的结果都是一样的,导致动画看起来没有变化。请问这是什么问题呢?

Tension=100 #N
mu = .1 #

v= (Tension/mu)**0.5
L=1 #string length
Nx=200
SimulationTime=.2


x = linspace(0, L, Nx+1) # Mesh points in space
 
dx = x[1] - x[0]
Nt = 1000
dt = SimulationTime/Nt
t = linspace(0, Nt*dt, Nt+1)
V = v*dt/dx
V2 = V**2

Lambda=.1

u = zeros(Nx+1) # t
u_1 = zeros(Nx+1) # t-1dt
u_2 = zeros(Nx+1) # t-2dt

# function that can be used to describe initial conditions
#Note - I should be zero at the boundary and dI/dx should not have a very large value at the boundary.  
I = lambda x: 1e-3*x*(1-x)  
#I = lambda x: 1.e-3*np.sin(2*np.pi*x/Lambda) 
#I = lambda x: 1.e-3*signal.gausspulse(x-Nx*dx/2, fc=10)

for i in range(0,int(Nx)):
     u_1[i] = I(x[i]) #u_1 - solution at time n-1 - Here we put initial conditions for t=0; 


n=0  #First timepoint
for i in range(1, Nx):
     u[i] = u_1[i] -0.5*V2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])  


# Use variable user_data to store u  for each iteration over time.          
user_data = {}
user_data['x'] = x
user_data['u'] = []
#(n, u)


u_2[:] = u_1
u_1[:] = u
for n in range(1, Nt):
     #Compute the values of the whole string at a given timepoint
     for i in range(1, Nx):
          u[i] = -u_2[i]+2*u_1[i]+V2*(u_1[i+1]-2*u_1[i]+u_1[i-1])
     u[0]=0
     u[Nx]=0
     u_2[:] = u_1
     u_1[:] = u

     user_data['u'].append(u)

我尝试了很多不同的调整,但结果似乎还是没有变化。

这是动画的部分:

from matplotlib.animation import FuncAnimation


# Example data (list of lists)
data = user_data['u']
data.pop(0)


# Create figure and axis
fig, ax = plt.subplots()
x_data = np.arange(len(data[0]))

# Plot the initial data
line, = ax.plot(x_data, data[0], label='Series 1')

# Customize plot
ax.set_xlim(0, len(data[0]) - 1)
#ax.set_ylim(min(min(data)), max(max(data)))
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.legend()

# Animation update function
def update(frame):
    line.set_ydata(data[frame])
    return line,

# Create animation
animation = FuncAnimation(fig, update, frames=len(data), interval=1000, blit=True)

# Show the plot
plt.show()

from IPython.display import HTML
HTML(animation.to_jshtml())

1 个回答

1

你在每个时间点都在添加u的引用。所以每次添加的都是u[]的最新值!

你需要添加一个副本:

user_data['u'].append(u.copy())

撰写回答