Construct the solution as $u(x) = \sup\{v(x) : v \text{ is subharmonic with } v \leq g \text{ on } \partial\Omega\}$. This supremum over subharmonic functions yields the unique harmonic function with boundary values $g$, provided the domain and data are regular.

For Laplace equation, single-layer potential: $u(x) = \int_{\partial\Omega} \Phi(x-y)\mu(y)\,dS_y$ and double-layer potential: $v(x) = \int_{\partial\Omega} \frac{\partial\Phi}{\partial n_y}(x-y)\nu(y)\,dS_y$ provide integral representations. Solving for densities $\mu, \nu$ to match boundary conditions leads to boundary integral equations.

Minimize the Dirichlet energy: $J[u] = \frac{1}{2}\int_\Omega |\nabla u|^2\,dx - \int_\Omega fu\,dx$ over $u \in H^1_0(\Omega)$. The minimizer satisfies $-\Delta u = f$ (Euler-Lagrange equation). This direct method establishes existence via:

1. Show $J$ is bounded below and coercive
2. Take minimizing sequence
3. Use weak compactness in $H^1_0$
4. Show limit is the minimum (using lower semicontinuity)

For $-\Delta u + cu = f$, use Green's function: $u(x) = \int_\Omega G(x,y)f(y)\,dy + \int_{\partial\Omega} u(y)\frac{\partial G}{\partial n_y}(x,y)\,dS_y$ Determining $G$ for specific domains (disk, half-space, exterior of sphere) gives explicit formulas.

Discretize variational formulation using piecewise polynomial spaces $V_h \subset H^1_0$. Solve: $\int_\Omega \nabla u_h \cdot \nabla v_h = \int_\Omega fv_h$ for all $v_h \in V_h$. This reduces to sparse linear systems with convergence $\|u - u_h\|_{H^1} = O(h^k)$ for degree $k$ elements.