Weak Formulation of Helmholtz equation with a complex coefficient and complex source term

Question

I'm trying to solve

$$\nabla^2u + ku = f\text{ in } \Omega,$$ $$u=0 \text{ on }\partial\Omega,$$ where $\Omega$ is a unit square and $k,f$ have real and imaginary components as $k=k_r+k_ij$ and $f(x,y)=f_r(x,y) + j\cdot f_i(x,y)$, respectively.

My approach to solve this problem is to assume the solution has both a real and imaginary component as well: $u=u_r + ju_i$. Substituting this int the PDE above, we get:

$$\nabla^2u_r + j\nabla^2u_i+ (k_ru_r-k_iu_i)+j(k_iu_r+k_ru_i) = f_r+jf_i$$

Equating reals with reals and imaginaries with imaginaries, we obtain a system of two coupled PDE's:

$$\nabla^2u_r + (k_ru_r-k_iu_i) = f_r$$ $$\nabla^2u_i + (k_iu_r+k_ru_i) = f_i$$

Multipling each equation by a separate test function $v_r$, $v_i$, respectively, and integrating to get the weak forms, I get

$$-(\nabla u_r,\nabla v_r) + k_r(u_r,v_r)-k_i(u_i,v_r) = \int_\Omega f_rv_r$$ $$-(\nabla u_i,\nabla v_i) + k_i(u_r,v_i)+k_r(u_i,v_i) = \int_\Omega f_iv_i$$

Adding the two equations I get a single equation: $$-(\nabla u_r,\nabla v_r) + k_r(u_r,v_r)-k_i(u_i,v_r) -(\nabla u_i,\nabla v_i) + k_i(u_r,v_i)+k_r(u_i,v_i) = \int_\Omega f_rv_r+\int_\Omega f_iv_i$$

I then choose a finite dimensional space of piecewise polynomial functions V such that $u_r,v_r,u_i,v_i\in V$. That is, I choose all trial and test functions to come from the same space.

I've tried using Method of Manufactured Solutions to test my implementation in two cases. In case 1, I assume that the solution is real and of the form $u=sin(\pi x)sin(\pi y)$, which yields a source term $f = (-2\pi^2+2j)sin(\pi x)sin(\pi y)$. In case 2, I assume the solution has both real and imaginary parts and is of the form $u=sin(\pi x)sin(\pi y) + j sin(\pi x)sin(\pi y)$, which yields a source term $f= [(-2\pi^2-2) + j(-2\pi^2+2) ]sin(\pi x)sin(\pi y)$.

For case 1, the error decreases as I refine both the mesh size and the polynomial order, as expected. However, in case 2, I do not get convergence at all. Clearly, it could be a programming error, but I want to make sure that my approach is sound. In particular:

1. Is it ok to separate the complex equation into a system of real PDE's?

2. Am I wrong for multiplying real test functions into the separated equations? Should I have multiplied the complex equation by a complex test function of the form $v_r+j_vi$?

3. Is it ok to add the two weak forms together to come up with a single variational expression? I've seen this done when using mixed finite element spaces for the stokes equation, but I'm not entirely sure that this is ok for this problem.

Any help would be greatly appreciated.

Answer 1

Let $A$ and $B$ be real matrices, and $b$ and $c$ be real vectors. Then if we define $z$ as the complex solution to a complex system of equations, $$(A+jB)z = (b+jc)$$ then $z$ is equivalently given as $z=x+jy$, where $$\begin{bmatrix}A & -B \\ B & A\end{bmatrix}\begin{bmatrix}x \\ y\end{bmatrix} = \begin{bmatrix}b \\ c\end{bmatrix}$$ These simple facts extend to linear operator equations. That is, suppose that $A+jB$ were the complex Helmholtz differential operator that you are trying to invert. Then you can follow these identical steps to obtain two sets of real, coupled operator equations.

To answer your original equations:

1. Yes. This is simply the matrix representation of complex numbers extended to linear operators.
2. This step is fine, because complex testing / basis functions also have real representations, and can be split in the same way that the equations themselves are split.
3. This is where you've made a mistake. To give a finite-dimensional analogy, you have effectively collapsed the 2x2 system shown above (with $2n$ real equations and $2n$ real unknowns) into $(A+B)x+(A-B)y = b+c$, which is an underdetermined system with $n$ real equations and $2n$ real unknowns. You then made this problem well-defined by setting $y=0$. So you're really solving $(A+B)x = (b+c)$, and this is right only when $B=0$ and $c=0$.

If you remove that last step where you have added the two equations together, the system should converge for all complex cases.

Answer 2

Your step 2 is where I have doubts. From what I experienced, you could do go different approaches here, both of which I have seen implemented and give correct results.

Approach 1 (this is what you outlined): Split the original complex equation in a real and an imaginary part. The thought here is that both, the real and the imaginary part, have to equate to zero as to find a valid solution. Then multiply each of the two equations with a separate test-function.

Approach 2: Multiply the complex equation you want to solve by a complex valued test-function.

Example: Assume you want to solve the equation $R+jI=0$, where $R$ is the real and $I$ the imaginary part of the equation.

Approach 1 would yield (with test-functions $v_R, v_I$): \begin{align} v_R R &= 0 \\ v_I I &= 0 \end{align}

Approach 2 would yield the following: \begin{align} (v_R + jv_I)(R+jI) = 0 \\ v_R R - v_I I + j(v_IR + v_RI) = 0 \end{align} Assuming that each part, real and imaginary, has to equate to zero, this will yield: \begin{align} v_R R - v_I I = 0 \\ v_IR + v_RI = 0 \end{align} which is obviously not the same as approach 1. To me, it seems like approach 1 is equating the real and imaginary part of the original equation to zero, while approach 2 is the projection of the complex valued original equation into a complex valued test-function space, and then equation that to zero. However, I have no idea what the ramifications of both approaches are.