Как правильно решить одномерное волновое уравнение, чтобы получить профиль смещения (периодическая проблема граничных условий)?

Question

Я пытаюсь решить 1D волновое уравнение для кучи с периодической БК (периодической нагрузкой).

Я почти уверен в своих формулах дискретизации. Единственное, в чем я не уверен, так это периодическое время до н.э. и время (t) там ==> sin(omega*t),

Когда я настраиваю его так, как оно есть сейчас, это дает мне странный профиль смещения. Однако, если я настрою это sin(omega*1) или же sin(omega*2)... и т. д., это напоминает синусоидальную волну, но в основном это означает, что sin(omega*t)т.е. sin(2*pi*f*t), равен 0, когда t является целочисленным значением..

Я кодирую все в Jupyter Notebook вместе с частью визуализации, но решение не близко к распространяющейся синусоиде.

Вот соответствующий код Python:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

def oned_wave(L, dz, T, p0, Ep, ro, f):

    """Solve u_tt=c^2*u_xx on (0,L)x(0,T]"""
    """Assume C = 1"""

    p = p0
    E = Ep
    ro = ro
    omega = 2*np.pi*f
    c = np.sqrt(E/ro)
    C = 1                                     # Courant number

    Nz = int(round(L/dz))
    z = np.linspace(0, L, Nz+1)               # Mesh points in space
    dt = dz/c                                 # Time step based on Courant Number
    Nt = int(round(T/dt))
    t = np.linspace(0, Nt*dt, Nt+1)           # Mesh points in time
    C2 = C**2                                 # Help variable in the scheme

    # Make sure dz and dt are compatible with z and t
    dz = z[1] - z[0]
    dt = t[1] - t[0]

    w = np.zeros(Nz+1)                        # Solution array at new time level
    w_n = np.zeros(Nz+1)                      # Solution at 1 time level back
    w_nm1 = np.zeros(Nz+1)                    # Solution at 2 time levels back

    # Set initial condition into w_n
    for i in range(0,Nz+1):
        w_n[i] = 0

    result_matrix = w_n[:]                    # Solution matrix where each row is displacement at given time step

    # Special formula for first time step
    for i in range(1, Nz):
        w[i] = 0.5*C2 * w_n[i-1] + (1 - C2) * w_n[i] + 0.5*C2 * w_n[i+1]

    # Set BC
    w[0] = (1 - C2) * w_n[i] + C2 * w_n[i+1] - C2*dz*((p*np.sin(omega*dt))/E) # this is where, I think, the mistake is: sin(omega*t)
    w[Nz] = 0

    result_matrix = np.vstack((result_matrix, w)) # Append a row to the solution matrix

    w_nm1[:] = w_n; w_n[:] = w                # Switch variables before next step

    for n in range(1, Nt):
        # Update all inner points at time t[n+1]
        for i in range(1, Nz):
            w[i] = - w_nm1[i] + C2 * w_n[i-1] + 2*(1 - C2) * w_n[i] + C2 * w_n[i+1]

        # Set BC
        w[0] = - w_nm1[i] + 2*(1 - C2) * w_n[i] + 2*C2 * w_n[i+1] - 2*dz*((p*np.sin(omega*(dt*n)))/E) # this is where, I think, the mistake is: sin(omega*t)
        w[Nz] = 0

        result_matrix = np.vstack((result_matrix, w)) # Append a row to the solution matrix

        w_nm1[:] = w_n; w_n[:] = w            # Switch variables before next step

    return result_matrix

Мой блокнот Jupyter.