I've got an inhomogeneous PDE of the following form:
$$\alpha\partial^2_xu+\partial_tu=f$$
with $\alpha<0$ and a source term $f$. I descretise $u$ according to $u_{m,n}=u(m\Delta t,n\Delta x)$ ($f$ analogously) and take as the derivatives the finite-difference formulas
\begin{align}
(\partial_tu)_{m,n}&=\frac{u_{m+1,n}-u_{m,n}}{\Delta t}\quad\text{(foward difference)},\\\\
(\partial_x^2u)_{m,n}&=\frac{u_{m+1,n}-2u_{m,n}+u_{m-1,n}}{(\Delta x)^2}\quad\text{(two-fold central difference)}.
\end{align}
Now I would like to assess the numerical stability of this scheme. Therefore, I make the ansatz $u=\xi^m\operatorname{e}^{\mathrm{i}kn\Delta x}$ with temporal amplification factor $\xi$ and spatial mode $k$ (von Neumann stability analysis) and I yield the following equation for $\xi$:
$$\xi^{m+1}\operatorname{e}^{\mathrm{i}kn\Delta x}=\xi^m\operatorname{e}^{\mathrm{i}kn\Delta x}\bigg[1+\frac{4\alpha\Delta t}{(\Delta x)^2}\sin^2\bigg(\frac{k\Delta x}{2}\bigg)\bigg]+\Delta tf_{m,n}.$$
I would like to solve this for $\xi$ in order to assess the stability of my numerical scheme, but I don't know how to deal with the inhomogeneous summand of $\Delta tf$. (In the case of $f=0$ I could directly solve this equation for $\xi$, but the inhomogeneous summand spoils it.)
Best Answer
General Result
If the linear PDE is well-posed, it is enough to determine the stability of your numerical scheme by only studying the numerical discretization of the homogenous PDE.
The Details...
As long as the inhomogenous function $f$ isn't too nasty (i.e. the problem is well-posed), it suffices determine the numerical stability of your scheme by just working with the homogenous pde. To see why, let's consider the following linear inhomogenous differential equation defined on \begin{align} u_t = Du + f(x,t) \end{align} where $D$ is the linear differential operator and $f(x,t)$ is the inhomogenous forcing term. Let $u^n_j = u(j\Delta x ,n\Delta t)$ denote the discretization of $u$. Using the forward Euler method to discretize this equation in time, the discretized version of the PDE becomes $$\frac{u^{n+1}_j - u^n_j}{\Delta t} = Au_j^n + f^n_j$$ where $A$ is the matrix which approximates the linear differential operator $D$ to some desired level of accuracy. Furthermore, we can represent the grid functions above using Fourier analysis. We have $$u_j^n = \int_{-\frac{\pi}{\Delta x}}^{\frac{\pi}{\Delta x}} \hat{u}(\xi)e^{i\xi\Delta x j}d\xi, \quad f^n_j = \int_{-\frac{\pi}{\Delta x}}^{\frac{\pi}{\Delta x}}\hat{f}(\xi)e^{i\xi \Delta xj}d\xi.$$ Plugging in this expression into the discrete PDE above, we obtain $$\int_{-\frac{\pi}{\Delta x}}^{\frac{\pi}{\Delta x}} \left(\frac{\hat{u}^{n+1} - \hat{u}}{\Delta t} - g(\xi,\Delta x)\hat{u}^n - \hat{f}^n\right)e^{i\xi\Delta x j}d\xi = 0$$ where $g(\xi,\Delta x)$ is the amplification factor resulting from applying the matrix $A$ to the eigenvector $e^{i\xi\Delta x j}$. Since the entire integral is equal to zero, we must have that the integrand itself is zero and so $$\hat{u}^{n+1} = \hat{u}^n\left(1 + \Delta t g(\xi,\Delta x)\right) + \Delta t \hat{f}^n.$$ This implies that $$\hat{u}^n = \left(1 + \Delta t g(\xi,\Delta x)\right)^n\hat{u}^0 + \Delta t \sum_{k=0}^{n-1}\hat{f}^k(1+g(\xi,\Delta x))^{n-1-k}$$ Now, if we assume that $|1+\Delta tg(\xi,\Delta x)| \le 1$ where $|\cdot|$ is the complex modulus, we can multiply the above equation by its complex conjugate $\bar{\hat{u}}^n$, perform some somewhat tedious algebraic manipulation, and apply the triangle inequality to obtain that $$|\hat{u}|^2\le |\hat{u}^0|^2 + 16\Delta t|\hat{u}^0|\sum_{k=0}^{n-1}|\hat{f}^k| + \Delta t^2\sum_{k=0}^{n-1}|\hat{f}^k|^2$$ Integrating both sides of this equation from $\frac{-\pi}{\Delta x}$ to $\frac{\pi}{\Delta x}$ and applying the Cauchy-Schwarz inequality then gives $$||\hat{u}^n||_{2}^2 \le ||\hat{u}^0||_2^2 + 16\Delta t||\hat{u}^0||_2\sum_{k=0}^{n-1}||\hat{f}^k||_2 + \Delta t^2\sum_{k=0}^{n-1}||\hat{f}^k||_2^2$$ Therefore, thanks to Parseval's theorem, we may conclude that the solution is stable so long as $n\Delta t \sim \mathcal{O}(1)$, $f^n_j \in l_2$ ($f$ isn't too nasty) for each $n$ and our initial data $u^0_j\in l_2$. Since we only required that the modulus $|1+\Delta tg(\xi,\Delta x)|\le 1$, which is the same requirement needed for the homogenous problem, it is enough to deduce stability just by studying the numerical discretization of the homogenous problem. We applied the forward Euler time discretization to arrive at this conclusion, but you can easily apply this result to other time discretizations by applying slight modifications.