python - 用 python 和 animate 求解一维波浪方程

Question

我正在尝试求解一维波动方程，并为数值计算解决方案和动画程序编写代码，将数据保存在文件中。我不知道如何修复错误并最终获得工作代码。

u_tt = a**2 * u_xx + f(x,t)

当输入附加函数和非零初始条件和边界条件时，程序有必要求解方程，图形可视化并将数据保存到文件中。

所以我附上了我的代码（Python 3.9）和错误消息：

import numpy as np
import math
import matplotlib.pyplot as plt
import os
import time
import glob


def sol(I, V, f, a, L, C, T, U_0, U_L, dt, user_func = None):
    """
    solver for wave equation
    u_tt = a**2*u_xx + f(x,t) (0,L) where u=0 for
    x=0,L, for t in (0,T].
    :param I:
    :param V:
    :param f:
    :param a:
    :param L:
    :param C:
    :param T:
    :param U_0:
    :param U_L:
    :param dt:
    :param user_func:
    :return:
    """

    nt = int(round(T / dt))
    t = np.linspace(0, nt * dt, nt + 1)  # array for time points
    dx = dt * a / float(C)
    nx = int(round(L / dx))
    x = np.linspace(0, L, nx + 1)  # array for coord points
    q = a ** 2
    C2 = (dt / dx) ** 2
    dt2 = dt * dt

    # --- checking f(x,t) ---
    if f is None or f == 0:
        f = lambda x, t: 0  

    # --- check the initial conds dU(x,0)/dt ---
    if V is None or V == 0:
        V = lambda x: 0

    # boundary conds
    if U_0 is not None:
        if isinstance(U_0, (float, int)) and U_0 == 0:
            U_0 = lambda t: 0
    if U_L is not None:
        if isinstance(U_L, (float, int)) and U_L == 0:
            U_L = lambda t: 0

    # ---  allocate memory  ---
    u = np.zeros(nx + 1)  
    u_n = np.zeros(nx + 1)  
    u_nm = np.zeros(nx + 1)  

    # --- valid indexing check ---
    Ix = range(0, nx + 1)
    It = range(0, nt + 1)

    # --- set the boundary conds ---
    for i in range(0, nx + 1):
        u_n[i] = I(x[i])

    if user_func is not None:
        user_func(u_n, x, t, 0)

    # --- finite difference step ---
    for i in Ix[1:-1]:
        u[i] = u_n[i] + dt * V(x[i]) + 0.5 * C2 * (0.5 * (q[i] + q[i + 1]) * (u_n[i + 1] - u_n[i]) -
               0.5 * (q[i] + q[i - 1]) * (u_n[i] - u_n[i - 1])) + 0.5 * dt2 * f(x[i], t[0])

    i = Ix[0]
    if U_0 is None:
        # set the boundary conds (x=0: i-1 -> i+1  u[i-1]=u[i+1]
        # where du/dn = 0, on x=L: i+1 -> i-1  u[i+1]=u[i-1])
        ip1 = i + 1
        im1 = ip1  # i-1 -> i+1
        u[i] = u_n[i] + dt * V(x[i]) + \
               0.5 * C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1])
                           * (u_n[i] - u_n[im1])) + 0.5 * dt2 * f(x[i], t[0])
    else:
        u[i] = U_0(dt)

    i = Ix[-1]
    if U_L is None:
        im1 = i - 1
        ip1 = im1  # i+1 -> i-1
        u[i] = u_n[i] + dt * V(x[i]) + \
               0.5 * C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1]) * (u_n[i] - u_n[im1])) + \
               0.5 * dt2 * f(x[i], t[0])
    else:
        u[i] = U_L(dt)

    if user_func is not None:
        user_func(u, x, t, 1)

    # update data
    u_nm, u_n, u = u_n, u, u_nm

    # --- time looping ---
    for n in It[1:-1]:
        # update all inner points
        for i in Ix[1:-1]:
            u[i] = - u_nm[i] + 2 * u_n[i] + \
                   C2 * (0.5 * (q[i] + q[i + 1]) * (u_n[i + 1] - u_n[i]) -
                         0.5 * (q[i] + q[i - 1]) * (u_n[i] - u_n[i - 1])) + dt2 * f(x[i], t[n])

        #  --- set boundary conds ---
        i = Ix[0]
        if U_0 is None:
            # set the boundary conds
            # x=0: i-1 -> i+1  u[i-1]=u[i+1] where du/dn=0
            # x=L: i+1 -> i-1  u[i+1]=u[i-1] where du/dn=0
            ip1 = i + 1
            im1 = ip1
            u[i] = - u_nm[i] + 2 * u_n[i] + \
                   C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1]) * (u_n[i] - u_n[im1])) + \
                   dt2 * f(x[i], t[n])
        else:
            u[i] = U_0(t[n + 1])

        i = Ix[-1]
        if U_L is None:
            im1 = i - 1
            ip1 = im1
            u[i] = - u_nm[i] + 2 * u_n[i] + \
                   C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1]) * (u_n[i] - u_n[im1])) + \
                   dt2 * f(x[i], t[n])
        else:
            u[i] = U_L(t[n + 1])

        if user_func is not None:
            if user_func(u, x, t, n + 1):
                break

        u_nm, u_n, u = u_n, u, u_nm

    return u, x, t



# --- here function for return functions ---
# return func(x)
def func(x):
    """
    :param x:
    :return:
    """

    return # expression

# start simulate and animate or visualisation and savin the data from file
def simulate(
        I, V, f, a, L, C, T, U_0, U_L, dt,  # params
        umin, umax,  # amplitude
        animate = True,  # animate or not?
        solver_func = sol,  # call the solver
        mode = 'plotter',  # mode: plotting the graphic or saving to file
):
   # code for visualization and simulate
   ...........
    # start simulate
    solver_func(I, V, f, a, L, C, T, U_0, U_L, dt, user_func)

    return 0

def task( ):
    '''
    test tasking for solver and my problem
    :return:
    '''

    I
    L = 1
    a = 1
    C = 0.85
    T = 1
    dt = 0.05
    U_0, U_L, V, f
    umax = 2
    umin = -umax

    simulate(I, V, f, a, L, C, T, U_0, U_L, dt, umax, umin, animate = True, solver_func = sol, mode = 'plotter',)


if __name__ == '__main__':
    task()

我得到同样的错误：

File "C:\\LR2-rep\wave_eq_1d.py", line 102, in sol
u[i] = u_n[i] + dt * V(x[i]) + 0.5 * C2 * (0.5 * (q[i] + q[i + 1]) * (u_n[i + 1] - u_n[i]) -

TypeError：“int”对象不可下标

我理解错误的含义，但我不明白如何修复它，并且将近两个星期我一直无法编写程序......我请求帮助解决这个问题！非常感谢您！

标签： pythonpde