0
$\begingroup$

I'm working with simulating both the heat and wave equation in 2D in a Python code. When simulating the heat equation, I learned about the CFL which I used to get a numerical stable solution.

I found CFL for the 2D wave equation as well, but my solution keeps being numerically unstable no matter how I change the time increments $dt$. I'm presenting the code here:

# Libraries
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib as mpl



def f(X,Y): # f(x,y,0)
    return X**2+Y**2

def g(X,Y): # g(x,y,0)
    return X

def laplacian(arr, row,col, dx ):
    return (arr[row + 1, col]
                            + arr[row - 1, col] + arr[row, col + 1]
                            + arr[row, col - 1]- 4*arr[row,col])/(dx**2)

def simulation(rect,hs,BC, c, frames, eq, eps = 1e-10):

    """ Simulation of the heat equation on a
    2D rectangular grid using the Finite Difference Method """

    X, Y = np.meshgrid(np.linspace(rect[0],rect[1],hs[0]),
                       np.linspace(rect[2],rect[3],hs[1])) # 2D meshgrid

    # Initial Conditions
    Z_init = f(X,Y) # u(x,y,0) = f(x,y,0)
    Z_dot_init = g(X,Y) # u_t(x,y,0) = g(x,y,0)

    zmax = max(Z_init.max(), BC[0], BC[1], BC[2], BC[3])
    zmin = min(Z_init.min(), BC[0], BC[1], BC[2], BC[3])
    zs = [Z_init]

    # Boundary Conditions
    if eq == "Heat":
        Z_init[0] = np.ones(len(Z_init[0])) * BC[0]
        Z_init[-1] = np.ones(len(Z_init[-1]))* BC[1]
        Z_init[:,0] =  np.ones(len(Z_init[:,0]))* BC[2]
        Z_init[:, -1] = np.ones(len(Z_init[:, -1])) * BC[3]

    # Figure settings
    fig = plt.figure()
    ax = plt.axes(projection='3d')
    ax.axes.set_xlim3d(rect[0] + eps, rect[1] - eps)
    ax.axes.set_ylim3d(rect[2] + eps, rect[3] - eps)
    ax.axes.set_zlim3d(zmin, zmax)
    plt.rcParams['mathtext.fontset'] = 'stix'
    plt.rcParams['font.family'] = 'STIXGeneral'

    # Plots for lines constituting the rectangle
    ax.plot([rect[0], rect[1]], [rect[2], rect[2]], [BC[0], BC[0]], color='black', linewidth=2)
    ax.plot([rect[1], rect[1]], [rect[2], rect[3]], [BC[3], BC[3]], color='black', linewidth=2)
    ax.plot([rect[0], rect[0]], [rect[2], rect[3]], [BC[2], BC[2]], color='black', linewidth=2)
    ax.plot([rect[0], rect[1]], [rect[3], rect[3]], [BC[1], BC[1]], color='black', linewidth=2)

    ax.plot([rect[0], rect[0]], [rect[2], rect[2]], [BC[0], BC[2]], color='black', linewidth=2)
    ax.plot([rect[1], rect[1]], [rect[2], rect[2]], [BC[0], BC[3]], color='black', linewidth=2)
    ax.plot([rect[0], rect[0]], [rect[3], rect[3]], [BC[1], BC[2]], color='black', linewidth=2)
    ax.plot([rect[1], rect[1]], [rect[3], rect[3]], [BC[1], BC[3]], color='black', linewidth=2)
    ax.set_title('Temperature development in a \n'
                 'rectangular room', fontsize = 18, fontname = 'STIXGeneral')

    # Infitesimals
    dx = (rect[1] - rect[0]) / ((hs[0] - 1))
    if eq == 'Heat':
        dt = dx**2/(10*c**2)
    else:
        dt = dx/(10*c)

    surf = ax.plot_surface(X, Y, Z_init, alpha=0.7, cmap=plt.cm.jet, vmin=zmin, vmax=zmax)
    cbar = fig.colorbar(surf)
    cbar.ax.set_ylabel('Temperature in ' + r'$^\circ\mathrm{C}$', rotation=270,fontsize = 14, labelpad=20)
    ax.set_axis_off()


    # Finite Difference Method
    if eq == 'Wave':

        for row in range(1, hs[0] - 1):
            for col in range(1, hs[1] - 1):
                Z_init[row,col] = Z_init[row,col] - 2 * laplacian(Z_dot_init,row,col, dx) + 1/2 * c ** 2 * dt ** 2 / (
                    dx) ** 2 * (Z_init[row+1, col] - 4 * Z_init[row, col]
                                + Z_init[row-1, col] + Z_init[row, col+1]
                                + Z_init[row, col-1])

        zs.append(Z_init)

        for iteration in range(2,frames):
            surf.remove()
            for row in range(1, hs[0]-1):
                for col in range(1,hs[1]-1):
                    Z_init[row,col] += 2 * zs[iteration - 1][row, col] - zs[iteration - 2][row, col] + c ** 2 * dt ** 2 / (
                            dx**2) * (zs[iteration - 1][row + 1, col] - 4 * zs[iteration - 1][row, col]
                                        + zs[iteration - 1][row - 1, col] + zs[iteration - 1][row, col + 1]
                                        + zs[iteration - 1][row, col - 1])

            zs.append(Z_init)
            surf = ax.plot_surface(X, Y, zs[iteration-2], alpha = 0.7, cmap =plt.cm.jet, vmin = zmin, vmax = zmax)
            plt.pause(dt * 10)

    if eq == 'Heat':

        for iteration in range(frames):

            surf.remove()
            for row in range(1, hs[0]-1):
                for col in range(1,hs[1]-1):
                    Z_init[row,col] = Z_init[row,col] + c**2 * dt / (dx)**2  *(Z_init[row+1,col] - 4*Z_init[row,col]
                                                   + Z_init[row-1,col] + Z_init[row,col+1]
                                                   + Z_init[row,col-1])

            # Plotting surface with slight pause
            surf = ax.plot_surface(X, Y, Z_init, alpha = 0.7, cmap =plt.cm.jet, vmin = zmin, vmax = zmax)
            plt.pause(dt)

simulation(rect = [-1,1,-1,1], hs = [100,100], BC = [0,0,0,0], c = 3, frames = 1000, eq = 'Wave')

I'm adding an image of the user instructions too:

enter image description here

Either way, I've tried searching for how to make my code numerically stable, but with no results. I tried changing the time steps to be even lower, I tried adding more elements into the grid, tried changing the maximum velcity, tried changing the initial conditions. I even tried removing the boundary conditions to only apply for the heat equation.

If you have any ideas on how to help me solve this, it'd be greatly appreciated.

Improve this question
$\endgroup$
9
  • $\begingroup$ What are your boundary conditions? How are you computing your Courant number? $\endgroup$
    – nicoguaro ♦
    Jan 29 at 21:07
  • $\begingroup$ @nicoguaro My boundary conditions are given by the array BC, in this case, all boundary conditions are 0 meaning that the wave always is 0 at the boundaries. I compute my Courant number as $dt = K * dx/(v_{max})$ and choosing $K<1$ and $v_{max} = c$. See "infitesimals" in my code. $\endgroup$
    – Tanamas
    Jan 30 at 8:10
  • 3
    $\begingroup$ Big error: You overwrite the old values with the new values and then use the new values to update later new values. I'm not sure if that is sufficient to destabilize the solution process, but it will lead to a wrong dynamic. $\endgroup$
    – Lutz Lehmann
    Jan 30 at 11:22
  • $\begingroup$ @LutzLehmann Wow thank you really! I have a question though regarding this, if I for instance make a list with an array that I later on change in a loop will this change go through to the array although it's "inside" the array in Python? Either way, I managed to solve it my just creating a temporary array. I appreciate your help :) $\endgroup$
    – Tanamas
    Jan 30 at 13:59
  • 1
    $\begingroup$ Another variant that gives not exactly the same method but is correct of the same order is to organize two sweeps in a checkerboard, red-blue, leapfrog pattern. So in a first sweep you update all red/white fields and in a second all blue/black fields. This gives something like a combination of a forward and a backward Euler step. $\endgroup$
    – Lutz Lehmann
    Jan 30 at 15:20

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.