I have the following set of non-dimensional equations that I am trying to solve using the finite difference method:
$$ u_{rr} + {1 \over r} u_r+u_{zz}-{1\over c^2}u_{tt}-f=0 \tag 1$$ $$ f_{rr} + {1 \over r} f_r+f_{zz}-{c^2}f_{tt}=0 \tag 2$$ $$ V_{rr} + {1 \over r} V_r - {1\over r^2}V+\big(2c^4-1\big)V_{zz}-c^2V_{tt}=0 \tag 3$$
- $u(r,z,t)$ is the longitudinal displacement
- $V(r,z,t)$ is a longitudinal derivative of radial displacement
- $f(r,z,t)$ is a function defined by equation $(1)$
- $c$ is the square root of the S-wave and P-wave speed ratio
Using the central difference approximation scheme, it is not hard to calculate the interior grid point values. The problem is, I have radial boundary conditions that cause my algorithm to be unstable. Because of that, my numerical solution does not converge.
At $r=r_{max}$ I use the following boundary condition to eliminate the ghost grid point of the longitudinal displacement:
$$ u_r=r_{max} \Bigg({1\over 2 - c^4}f + u_{zz} \Bigg) \tag 4 $$
Also, to eliminate the ghost grid point of the longitudinal derivative of radial displacement, I use the boundary condition:
$$ V_r = \big(2c^4-1\big)\bigg({1 \over r_{max} }V + u_{zz}\bigg) \tag 5 $$
As for the ghost grid point of $f$, I use the combination of the above boundary condition and the following boundary condition to eliminate it:
$$ f_r = \bigg(1 - {1 \over c^4}\bigg) \bigg(V_{rr} + {1 \over r_{max} }V_r - {1 \over r_{max}^2 }V - V_{zz} \bigg) \tag 6$$
No matter how small the time step, or how high the spatial step is, these boundary conditions cause instability. Below is a graph of longitudinal displacement from which this can be seen clearly in the form of oscillations:
I am trying to eliminate this instability, but I do not know how. My questions are:
$1$. Can I approximate these boundary conditions (or at least one of them) in some way so that the right-hand side of them does not contain derivatives anymore? For example, swap them with some polynomials, Fourier series, or other types of functions. Note that I know the boundary values at $z = 0$ and $z = 1$, as well as the initial time conditions.
$2$. Is there any other approach I can take to solve this boundary value problem?