Conceptual difference between strong and weak formulations
If we say a solution is weak/strong/classical/viscous, the following aspects are concerned (or more):
How we obtain the solution.
The regularity of the solution (how smooth this solution is, integrability, differentiability).
The solution satisfies the equation in what sense.
Weak solution:
- We can obtain the solution by Ritz-Galerkin formulation: find the minimizer of the following quadratic functional in an appropriate Hilbert space, $$ \mathcal{F}(u) = \frac{1}{2}\int_{\Omega} |\nabla u|^2 - \int_{\Omega} fu. $$
Smoothness depends on the right side data. If the $f\in H^{-1}$, then $u\in H^1$. If $f\in L^2$, then $u\in H^2_{loc}$. Moreover if $\Omega$ is $C^{1,1}$, we have an $H^2$-solution $u$ globally.
The solution satisfies the equation in distribution sense (see following explanation).
Why "weak":
The term "weak" normally refers to the 2 and 3: The solution $u$ is only in $H^1$ in the most general setting, this means that $u$ is the only differentiable once, notice $-\Delta$ has second partial derivative in it. The strong solution, however, indeed have twice differentiability, normally if we say $u$ is a strong solution, we mean that $u$ has $W^{2,p}$-regularity (Please refer to Gilbarg and Trudinger). The solution satisfies the equation only in the "weak" formulation $$ \int_{\Omega} \nabla u \cdot \nabla v \, dx = \int_{\Omega} fv \, dx \quad \forall v \in V, \tag{1} $$ where $V$ is certain Sobolev space.
Two ways to get this weak form: first is to write what condition the minimizer of $\mathcal{F}(u)$ must satisfy: if $u$ is a minimizer, then $$ \lim_{\epsilon \to 0}\frac{d}{d\epsilon} \mathcal{F}(u+\epsilon v) =0 $$ and the weak form of Euler-Lagrange equation is (1).
Another is multiplying the original equation by a test function then integration by parts. The intuition behind this should be Riesz representation theorem (at least to me it makes sense), we have: $$ \langle (-\Delta)u,v\rangle = l_u(v) = (u,v)_{V}, $$ from the differential operator $-\Delta$ $\to$ linear functional $l_u$ $\to$ representation using inner product $(\cdot ,\cdot)_V$. The inner product $(\cdot ,\cdot)_V$ on this Hilbert space $V$ is the left hand side of (1), if we make the test function space have zero boundary condition (We can use Poincaré inequality to prove the equivalence of the standard $H^1$-inner product). If you have taken any numerical PDE course in finite element, the professor would introduce Lax-Milgram theorem, and Lax-Milgram relies on Riesz.
Why weak form is useful in finite element method:
Short answer: Weak form is very handy in that it helps us formulate a linear equation system which can be solved by computer!
Long answer: The essential of Galerkin type approach is that we are exploiting the fact that the infinite dimensional Hilbert space has a set of basis $\{\phi_i\}_{i=1}^{\infty}$, if we can expand the $u$ in this basis: $$ u = \sum_{i=1}^{\infty} u_i\phi_i, $$ where $u_n$ is a number, pluggin back to (1), and let the test function $v$ run through all $\phi_j$ (same function, different subscript): $$ \int_{\Omega} \nabla (\sum_{i=1}^{\infty} u_i\phi_i) \cdot \nabla \phi_j \, dx =\sum_{i=1}^{\infty} u_i \int_{\Omega} \nabla \phi_i \cdot \nabla \phi_j \, dx = \int_{\Omega} f\phi_j \, dx \quad \forall j =1,2,\ldots. \tag{2} $$ We have an infinite dimensional linear equation system: $$ AU = F, $$ where $A_{ji} = \displaystyle\int_{\Omega} \nabla \phi_i \cdot \nabla \phi_j \, dx$, $U_i = u_i$, and $F_j = \displaystyle\int_{\Omega} f\phi_j\, dx$.
Finite element method essentially choose a finite dimensional subspace $V_h\subset V$ (may not be a subspace, please google Discontinuous Galerkin method), so that we approximate the solution in this finite dimensional subspace $V_h$! The summation in (2) does not have an infinite upper limit any more, instead there are finitely many $\phi_i$ and $v$ runs from $\phi_1$ to $\phi_N$, so that the linear system generated is still $AU = F$, but this time, it only has $N$ equations, and we can use computer to solve it.
First of all, it's fundamental for Finite Element Method because without that formulation, the correspondent numerical method to solve it would be more like Finite Difference.
When we consider a weak formulation of a PDE we are deliberately searching for solutions with less regularity conditions then the classical form imposes. We are trying to include in the class of solution for our PDE the candidates that almost satisfy the equation except for example for having a discontinuity on the derivative or a dirac jump. It wold be pretty interesting to have a a definition of solution and consequently a correspondent weak theory of derivation which allows us to include this cases in our set of solutions (they have frequently a physical interest). That is why pass to the weak form where derivatives a taken in the sense of the distributions (also called weak sense).
You can see a distribution as a extension of the definition of functions. See definition.
Another examples of issu that the distribution theory allows to treat is when the non homogeneous term and/or limit conditions are not regular. I give you a example of a real ODE to be solved in the distribution sense.
Consider $y' +2xy = \delta_0$, where $\delta_0$ is a distribution (delta de dirac) the classical theory don't covers this equation but it has a solution of the form $(H+c)e^{x^2}$( where $H$ is the Heavside function)
In weak form the equation is writen: $$\forall \phi \in \mathcal C_c (\mathbb R), \ \int_{\mathbb R} [ 2xy(x)\phi(x)-y(x)\phi '(x) ]~dx =\phi(0)$$
Apologies for reviving an old question, but I feel the answers don't sufficiently appeal to physical intuition for what "weak form" is accomplishing.
Consider systems in the real world where at an interface there can be abrupt changes to the material properties. These real physical systems can violate smoothness constraints for a PDE. You could say that the classical PDE for a system is too strict when looking at a real, arbitrary, physical system.
One simple example is the temperature diffusion equation:
$$q(x)=-\partial_x T(x)$$
where $q$ is the heat flux for a temperature profile $T$. The conservation of heat flux is:
$$\partial_xq(x)=0$$
which includes the second derivative of $T$. At a boundary where two materials have different coefficients of thermal conductivity, the first derivative of $T$ becomes discontinuous and the second derivative can't be found numerically.
We turn the equation into an integral equation to make our solution less strict. We could integrate over the entire domain of $x$ to satisfy the conservation equation, but this isn't strict enough to match with reality (since in reality the integral equation must equal zero at every point, not just simply zero).
Instead we split the integration region up into segments.
$$\int_0^1\partial_xq(x)dx=0 \rightarrow \int_0^{.01}\partial_xq(x)dx=0, \int_{.01}^{.02}\partial_xq(x)dx=0, ...\int_{.99}^1\partial_xq(x)dx=0 $$
Splitting this into more elements locks down the solution with more resolution. Essentially, the number of elements you have is now defining the strictness of your solution. In summary, the weak form allows us to represent and solve for real physical systems that aren't reliably represented but the strictness of the original PDE.
For more details, I would recommend COMSOL's blog series, beginning with The Strength of the Weak Form, where they describe the weak form and how it's implemented in their software.