Why the Green's function is zero on the Dirichlet boundary surface? Given a differential equation
$$
Lu(x) = f(x),
$$
where $L$ is the differential operator, and $f$ is a given function of $x$, the general solution of $u(x)$ is
$$
u(x) = \int G(x, s) f(s)\, ds
$$
$G(x, s)$ is the Green's function of $L$ and satisfies
$$
LG = \delta(x - s)
$$
If we have Dirichlet boundary conditions like
$$
u(a) = \alpha,
$$
$\alpha$ is a given number, this differential problem is called the Dirichlet problem.
My question is why $G(a, s) = 0$ with given a boundary condition $u(a) = \alpha$?
 A: The short answer is that your problem is linear non-homogeneous. You have two sources of inhomogeneity: the boundary condition ($\alpha$) and the sources ($f(x)$). Instead of solving the entire problem in one go, you want to treat the inhomogeneities one at a time and use superposition to reconstruct the original solution.
This is simply linear algebra. Take for example a $d$ dimensional vector space $E$ and $e^1,…,e^d$ basis vectors of the dual space. Say you want to solve the system of equations for $x\in E$:
$$
\langle e^1,x\rangle =a^1\\
…\\
\langle e^d,x\rangle =a^d
$$
and $a^1,…,a^d$ known constants. As you can see, you could try to solve the equation directly or construct the dual basis $e_1,…,e_d$ defined by:
$$
\langle e^i,e_j\rangle=\delta^i_j
$$
i.e. you normalize only one inhomogeneity and make the other constraints homogeneous. Your solution can this be directly solved to:
$$
x=\sum_i a^ie_i
$$
Formally, you can apply the same reasoning to your problem. The solution to your problem is therefore:
$$
u(x)=v(x)+\int G(x,s)f(s)ds \tag{1}
$$
with:
$$
LG(x,s) = \delta(x-s), \quad G(a,s)=0 \\
Lv(x) = 0, \quad v(a)=\alpha
$$
I think that part of your confusion comes from the fact that you forgot the first term.
Hope this helps.
Examples
Applying the method when $E$ is finite dimensional shows that it’s a pedantic way of introducing the inverse matrix. Take $E=\mathbb R^2$, and:
$$
e^1 = (L^1_1, L^1_2) \\ 
e^2 = (L^2_1, L^2_2) \\
$$
then let:
$$
e_1=\begin{pmatrix}
G_1^1 \\
G_1^2
\end{pmatrix} \quad e_2=\begin{pmatrix}
G_2^1 \\
G_2^2
\end{pmatrix}
$$
The duality condition is equivalent to the definition of the inverse matrix:
$$
LG=I_2
$$
In particular if you want to solve the system:
$$
Lx=a
$$
then you’ll use:
$$
x=Ga
$$
A more advanced example would be $L=-\frac{d^2}{dx^2}$ giving the bulk equation:
$$
-\frac{d^2u}{dx^2}=f \tag{2}
$$
and Dirichlet boundary conditions at $x=0,1$:
$$
u(0)=\alpha_0 \\
u(1)=\alpha_1 \tag{3}
$$
You then have:
$$
G(x,s)= \min[(1-s)x,s(1-x)]\\
v(x) = \alpha_1x+\alpha_0(1-x)
$$
Hope this helps.
