1
$\begingroup$

I have implemented a finite difference solver for the 1d wave equation with variable wave speed:

$$ u_{tt} = c(x)u_{xx}, \hspace{10mm}c(x) = \dfrac{6 -x^2}{2} \hspace{5mm} $$ on $-2 \leq x \leq 2, t > 0$ with initial conditions: $$u(x,0) = \cos(\dfrac{\pi x}{0.2}), \hspace{5mm} u_t(x,0) = 0 $$ for $-0.1 \leq x \leq 0.1$.

I use second order finite differences in order to solve for each grid point in time: $$\frac{u_j^{n+1} - 2u_j^n + u_j^{n-1}}{k^2} = c(x_j)\frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{h^2}$$

Rearranging for $u_j^{n+1}$: $$u_{j}^{n+1} = 2u_j^n - u_j^{n-1} + c(x_j)q^2(u_{j+1}^n - 2u_j^n + u_{j-1}^n)$$

Where $q = k/h$. The boundary condition to be enforced is the open condition, so $qu_{x} = u_{t}$ at $x = -2$, and $qu_{x} = -u_{t}$ at $x = 2$. Using the first order finite differences, we can obtain an equation for the ghost points at the boundaries:

$$ q(u_{1}^{n} - u_{-1}^n) = u_0^{t+1} - u_0^{t-1} $$ at $x = -2$.

Next, we can substitute this into the second order equation ($c(x_j)$ disappears since it is 1 at $\pm 2$) in order to end up with our iterative formula for the boundary at $x=-2$:

$$u_{0}^{n+1} = \dfrac{2u_0^n + (q - 1) u_0^{n-1} + 2r^2(u_{1}^n - u_0^n)}{q + 1}$$

The same can be done for the boundary at $x = 2$.

My python implementation is shown below.

def one_dimensional_wave_solver(x_steps, t_steps):
    # discrete points in space
    x = np.linspace(-2, 2, x_steps + 1) 
    # distance of step
    dx = x[1] - x[0]
    # discrete points in time
    t = np.linspace(0, 10, t_steps + 1)
    # time steps 
    dt = t[1] - t[0]

    # define cur_u and previous u
    cur_u = np.zeros(x_steps + 1)
    prev_u = np.zeros((x_steps + 1, t_steps + 2))

    # Define mesh constant
    c = 1
    mesh = (c * dt) / dx

    # Initial conditions 
    for i in range(0, x_steps + 1):
        if (-0.2 <= x[i] and x[i] <= 0.2):
            prev_u[i, 0] = np.cos((np.pi * x[i]) / (0.2))
        else:
            prev_u[i, 0] = 0

    # Iterate through time steps
    for t in range(0, t_steps + 1):
        # For the first time step
        if t == 0:
            for i in range(1, x_steps):
                cur_u[i] = prev_u[i, t] + (cx(x[i]) * (mesh ** 2) / 2) * (prev_u[i - 1, t] - 2 * prev_u[i, t] + prev_u[i + 1, t])
            # Boundary conditions
            cur_u[0] = prev_u[0, t] +  mesh * (prev_u[1, t] - prev_u[0, t])
            cur_u[-1] = prev_u[-1, t] - mesh * (prev_u[-1, t] - prev_u[-2, t])
            
        else:
            # equation for computing u
            for i in range(1, x_steps):
                cur_u[i] = 2 * prev_u[i, t] - prev_u[i, t - 1] + cx(x[i]) * (mesh ** 2) * (prev_u[i - 1, t] - 2 * prev_u[i, t] + prev_u[i + 1, t])
            # Boundary conditions
            cur_u[0] = (2 * prev_u[0, t] + (mesh - 1) * prev_u[0, t - 1] + (2 * mesh ** 2) * (prev_u[1, t] - prev_u[0, t])) / (mesh + 1)
            cur_u[-1] = (2 * prev_u[-1, t] + (mesh - 1) * prev_u[-1, t - 1] + 2 * (mesh ** 2) * (prev_u[-2, t] - prev_u[-1, t])) / (mesh + 1)

        # Swap variables
        prev_u[:, t + 1] = cur_u.copy()

    # Return final u value
    return prev_u[:, -1], prev_u

Using this allows me to get an almost fully transparent boundary. But there is still some slight reflection. Is there an issue with the first order equations for the boundaries?

enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ That doesn't look like a "slight" reflection. It's clearly visible. Without having examined the question in detail, I vote for error. (When I've done similar things years ago, I had really slight reflections of the order $10^{-6}$ or something like that). $\endgroup$
    – davidhigh
    Mar 26 at 14:09
  • 1
    $\begingroup$ 1. The initial conditions you give do not match the implemented ones, also you do not state the initial conditions outside the interval $[-0.1, 0.1]$ or $[-0.2, 0.2]$. 2. Why do you have only a plain $c$ in your wave equation, and not $u_{tt} = c{^\color{red} 2}(x) u_{xx}$? 3. In what sense is the B.C. you give the open condition? How about the simple outflow $u_{n+1} = u_n, u_{-1} = u_{0}$? 4. The image you provide looks quite underresolved, did you try smaller time-steps? 5. Related to 4: I am not sure that a CFL numbe rof 1 works here, since the max. of $c$ can be 3. $\endgroup$
    – Dan Doe
    Mar 31 at 9:43

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy