A general way to derive a weak form is to multiply a test function on both sides of the equation and then integrate them. The second step is to use some kind of divergence theorems to derive the weak solution such that the solution is some what not so smooth as in the strong form. For your question here, we can derive the weak form as follows:
Let the equation multiply a $\phi\in C_0^\infty$ on both sides and then integrate them, we get
\begin{equation}
\int_\Omega-\nabla\cdot (A(x)\nabla u)\phi dx=\int_\Omega\omega^2q(x)u\phi dx.
\end{equation}
Using the divergence theorem and the condition that $\phi$ vanishes on bdry, the LHS becoms
\begin{equation}
\int_\Omega (A(x)\nabla u)\cdot\nabla\phi dx.
\end{equation}
And we have the weak form
\begin{equation}
\int_\Omega (A(x)\nabla u)\cdot\nabla\phi dx=\int_\Omega\omega^2q(x)u\phi dx.
\end{equation}
for all $\phi\in C_0^\infty$.
Assuming that $\Omega $ is Lipschitz- continuous subset of $\mathbb{R}^n$ and making use of the identity I.2.17 and theorem 2.4 and 2.5 in [V. Girault and P. A. Raviart, Finite Element Methods for Navier-Stokes Equations. Theory and Algorithms, Springer Series in Computational Mathematics, 5 Springer, Berlin, 1986.]
Which state that\
- $\mathcal{D}(\bar{\Omega})^n$ is dense in $ H(\operatorname{div},\Omega)$\
and
- the map $\gamma: v \rightarrow v\cdot\eta|_{\Gamma_n}$ defined on $\mathcal{D}(\bar{\Omega})^n$ can be extended by continuity to a linear and continuous mapping from $H(\operatorname{div},\Omega)$ into $H^{-\frac{1}{2}}(\Gamma_n)$\
So, we have the following Green formula
$$(v,\nabla \phi)+(\operatorname{div} v,\phi)=(v\cdot \eta,\phi)_{\Gamma_n} \quad \forall v \in H(\operatorname{div},\Omega), \forall \phi \in H^1(\Omega)$$
Therefore, we have for every $v \in V_0$
\begin{align}
a(u,v)=\int_\Omega a(x)\nabla u(x) \cdot \nabla v(x)+\bigl (b(x) \cdot \nabla u(x) \bigr)v(x) \, \mathrm{dx} \\
= \int_\Omega -\operatorname{div} \big(a(x)\nabla u(x)\big)v(x)+\bigl (b(x) \cdot \nabla u(x) \bigr)v(x) \, \mathrm{dx}\\+\int_{\Gamma_n} a(x)\nabla u(x) \eta v(x) \mathrm{ds}
\end{align}
while $\eta$ is the outward unit normal vector.
Now, taking $v \in \mathcal{D}(\Omega)$ we have
$$ \int_\Omega -\operatorname{div} \big(a(x)\nabla u(x)\big)v(x)+\bigl (b(x) \cdot \nabla u(x) \bigr)v(x) \, \mathrm{dx}=\int_{\Omega} f(x)v(x)\, \mathrm{dx}$$
So,
$$ -\operatorname{div} \big(a(x)\nabla u(x)\big)+\bigl (b(x) \cdot \nabla u(x) \bigr)=f(x) \quad \text{a.e. in } \Omega$$
Which can be extended into all $\Omega$
Making use of this, we find
$$ a(v) \frac{\partial u}{\partial \eta}=g(x), \quad \forall x \in \Gamma_n$$
Finely, the equation is
$$
\left\{\begin{aligned}
-\operatorname{div} \big(a(x)\nabla u(x)\big)+\bigl (b(x) \cdot \nabla u(x) \bigr)=f(x) \quad \forall x \in \Omega\\
a(x) \frac{\partial u}{\partial \eta}=g(x), \quad \forall x \in \Gamma_n
\end{aligned}\right.
$$
The second one can be treated the same way to get
$$
\left\{\begin{aligned}
-\operatorname{div} \big(a(x)\nabla u(x)\big)-\bigl (b(x) \cdot \nabla u(x) \bigr)=f(x) \quad \forall x \in \Omega\\
a(x) \frac{\partial u}{\partial \eta}+b(x)u(x)\eta=g(x), \quad \forall x \in \Gamma_n
\end{aligned}\right.
$$
Best Answer
The weak form of a single PDE asserts an integral equality for all "test functions" in a suitable vector space.
If you retain the distinct test functions when summing several weak forms, so that we still quantify universally over them, then this summed-up form is equivalent to the system of weak forms because we could set all but one of the test functions to zero in order to recover a single weak form.