The B-S PDE is defined the space positive half space $\mathbb{R}^{+}_{S}\times\mathbb{R}^{+}_{t}=\{(S,t): S\in(0,\infty),t\in(0,\infty)\}$ by the equation
$$(1)\;\;\;\;V_{t}+\frac{1}{2}\sigma^{2}S^{2}V_{SS}+rSV_{S}-rV=0.$$
For various pay-off functions $f=f(S,t)$ (where the $S$ in the argument is interpreted as $\{S_{t}\}_{t\in\mathbb{R}}$; in many cases it will just be $t=T$, the value of $S$ at expiration), boundary conditions may be instituted and the domain of definition becomes $(0,S_{\text{max}})\times(0,T)$. We shall discuss boundary conditions in more detail later; for now we will work on the positive plane.
For $\xi\in\mathbb{R}$ we have
$$-\frac{1}{2}\sigma^{2}\xi^{2}<\theta\xi^{2}$$
for some appropriate positive parameter $\theta$ depending only on $\sigma$. Thus, (1) is uniformly backward parabolic.
I. Converting the B-S PDE to the heat equation (systematic procedure)
We begin our analysis by performing a sequence of change of independent and dependent variables in order to reduce the equation to a more simple one. We first deal with the variable co-efficients, whose special structure suggests the change of variables $S\mapsto e^{-S}$, or using a new variable name, $x=\log S$. The differential operators then become
$$\frac{\partial}{\partial S}=\frac{\partial}{\partial x}\frac{\partial{x}}{\partial S}=\frac{1}{S}\frac{\partial}{\partial x}$$
and
$$\frac{\partial^{2}}{\partial S^{2}}=\frac{\partial}{\partial S}\left[\frac{1}{S}\frac{\partial}{\partial x}\right]=\frac{1}{S^{2}}\left(\frac{\partial}{\partial x^{2}}-\frac{\partial}{\partial x}\right).$$
Substituting these into (1) yields
$$V_{t}+\frac{1}{2}\sigma^{2}(V_{xx}-V_{x})+rV_{x}-rV=0,$$
or
$$V_{t}+\frac{1}{2}\sigma^{2}V_{xx}+(r-\frac{1}{2}\sigma^{2})V_{x}-rV=0.$$
The reaction term $rV$ will lead to an exponential increase in the solution, and so we assume $V$ has the form $e^{rt}u$ for some to be determined function $u$. This leads to the change of variables $V\mapsto e^{rt}u$ and the equation becomes
$$e^{rt}\left(ru+u_{t}+\frac{1}{2}\sigma^{2}u_{xx}+(r-\frac{1}{2}\sigma^{2})u_{x}-ru\right)=0,$$
or
$$u_{t}+\frac{1}{2}\sigma^{2}u_{xx}+(r-\frac{1}{2}\sigma^{2})u_{x}=0.$$
(If $r=r(t)$ is not constant, but a deterministic function of $t$ alone, then we can solve the ODE $w_{t}-r(t)w=0$ and make the change of variables $V=w(t)u$ in order to eliminate the source term; we will assume however that $r\equiv\text{const}$ for the remainder of this article.)
We now deal with the drift term $(r-\frac{1}{2}\sigma^{2})u_{x}$ by switching to the moving frame $x'=x-(r-\frac{1}{2}\sigma^{2})t$ (i.e. switching to the characteristic coordinate system). We also at this stage change the structure of the equation to forward parabolic by changing the direction of evolution to forward time through the change of variables $t\mapsto -t$, or $t'=-t.$ We compute as before the new differential operators
$$\frac{\partial}{\partial x}=\frac{\partial}{\partial x'}\frac{\partial x'}{\partial x}+\frac{\partial}{\partial t'}\frac{\partial t'}{\partial x}=\frac{\partial}{\partial x'},$$
$$\frac{\partial^{2}}{\partial x^{2}}=\frac{\partial^{2}}{\partial x'^{2}},$$
$$\frac{\partial}{\partial t}=\frac{\partial}{\partial x'}\frac{\partial x'}{\partial t}+\frac{\partial}{\partial t'}\frac{\partial t'}{\partial t}=-(r-\frac{1}{2}\sigma^{2})\frac{\partial}{\partial x'}-\frac{\partial}{\partial t'}.$$
Substitution then yields
$$-(r-\frac{1}{2}\sigma^{2})u_{x'}-u_{t'}+\frac{1}{2}\sigma^{2}u_{x'x'}+(r-\frac{1}{2}\sigma^{2})u_{x'}=0$$
or
$$(2)\;\;\;\;u_{t'}-ku_{x'x'}=0.$$
This is the heat equation, with $k=\frac{1}{2}\sigma^{2}.$ It therefore suffices to study the properties of (1) by studying (2).
II. Uniqueness and Stability
We shall restrict our attention to the domain $D=\{(x,t):0\leq t\leq T,0\leq x\leq \ell\}.$ This is consistent with domain used in practice with solving (1). We have the following theorem
Maximum Principle. If $u$ solves the heat equation $u_{t}=ku_{xx}$, then $u$ achieves its maximum on the boundary lines $t=0$, $x=0$, or $x=\ell$. Moreover, if $u$ achieves its maximum in the interior of $D$ or on the line $t=T$, then $u$ is constant.The first assertion is what is known as the weak maximum principle, and this is what we shall prove. The second assertion is called the strong maximum principle since it implies the weak version; however, its proof requires tools from the theory of harmonic functions (e.g. the mean-value property) and some non-trivial analysis (cf. Evans). Analogous statements hold for the minimum by considering $-u$.
Proof. Recall from calculus that if a function $u$ obtains its maximum at an interior point of a set $E$, then $D_{i}u=0$ and $D^{2}_{i}u\leq0$. Therefore we have $u_{t}=0$ and $u_{xx}\leq0$. The idea is to show that $u_{xx}<0$ in order to obtain a contradiction, and we will do this with a perturbation.
Let $\epsilon>0$, put $M=\sup_{(x,t)\in\partial D-\{(x,t):t=T\}}u$, and define $v:=u+\epsilon x^{2}.$ The definition of $v$ shows
$$(3)\;\;\;\;v_{t}-kv_{xx}=u_{t}-ku_{xx}-2\epsilon k=-2\epsilon k<0.$$
Applying the argument in the first paragraph shows that $v$ does not attain its maximum in $D^{\circ}$. Suppose now $u$ obtains its maximum at a point $(x_{0},T)$ on the line $t=T$. Then $v_{x}=0,v_{xx}\leq0$ as before and
$$v_{t}(x_{0},T)=\lim_{\delta\to0^{+}}\frac{v(x_{0},T)-v(x_{0},T-\delta)}{\delta}\geq0,$$
a contradiction. Since $D$ is compact, the continuous function $v$ must obtain its maximum on $\partial D-\{(x,t):t=T\}.$ Therefore $v\leq M+\epsilon\ell^{2}$ on $D$. By sending $\epsilon\to0$, we obtain the assertion of the weak maximum principle for $u$.
This theorem immediately settles the uniqueness of the boundary value problem for (2), for suppose $u$ and $v$ are both solutions to (2) with identical boundary data. Then $w:=u-v$ is also a solution by linearity has vanishes on the boundary $x=0,x=\ell,t=0$. Thus $w\equiv0$ in all of $D$ by combining the assertions of the maximum and minimum principles. This implies $u\equiv v$ and the claim follows.
Another important route to establishing uniqueness is through the concept of conservation of energy.
Consider the "energy" integral
$$E(t)=\int_{0}^{\ell}|u(x,t)|^{2}\;dx$$
where $u$ vanishes on the spatial boundary $0\leq x\leq\ell$ and at the intial time $t=0$. Differentiating with respect to $t$ yields
$$E'(t)=\frac{d}{dt}\int_{0}^{\ell}u^{2}\;dx=\int_{0}^{\ell}2uu_{t}\;dx.$$
where differentiation under the integral sign is valid, since as will be clear in the next section when we solve the heat equation, the mapping $x\mapsto 2u(x,t)u'(x,t)$ is absolutely integrable on $[0,\ell]$ for every $t\in[0,T)$. Substituting $u_{t}=ku_{xx}$ and integrating by parts yields
$$E'(t)=2k\int_{0}^{\ell}uu_{xx}\;dx=2kuu_{x}\Big|_{x=0}^{x=\ell}-2k\int_{0}^{\ell}u^{2}_{x}\;dx=-2k\int_{0}^{\ell}u_{x}^{2}\;dx\leq0.$$
However, $E(0)=0$, and therefore $E\equiv0$ for all $t$. Consequently, $u\equiv0$ and by considering the difference of two solutions with identical boundary data, we recover uniqueness (initially, we have uniqueness only almost everywhere, but since solutions to the heat equation are smooth, we recover pointwise uniqueness.)
We continue our efforts to establish well-posedness by showing that the PDE is stable with respect to the initial and boundary data.
III. Existence and the fundamental solution
There are by now multiple approaches to solving the heat equation recorded in the literature. One of the simplest makes use of the Fourier transform. We begin by extending our spatial domain to all of $\mathbb{R}$ and enforcing a decay condition so that $u\in L^{2}(dx).$ Our convention for the (one-dimensional) Fourier transform is
$$\mathcal{F}(u)(\xi):=\hat{u}(\xi)=(2\pi)^{-1/2}\int_{-\infty}^{\infty}u(x)e^{-ix\xi}\;dx$$
so that
$$||\hat{u}||_{2}=||u||_{2}\;\;\;\;(\text{Plancheral})$$
along with the associated inversion formula
$$u(x)=\int_{-\infty}^{\infty}\hat{u}(\xi)e^{ix\xi}\;d\xi.$$
The Fourier transform drops polynomial weights in exchange for anti-differentiation, so we apply it to the spatial derivative $x$ in order to get
$$\hat{u}_{t}=-k|\xi|^{2}\hat{u}.$$
Note that we have used the relation $\mathcal{F}(\partial_{t}u)=\partial_{t}\mathcal{F}u,$ which holds because of our condition on $u$.
Solving this ODE (in $t$) yields
$$\hat{u}=\hat{u}(0)e^{-k\xi^{2}t}.$$
A simple computation involving substituting the above result into the inversion formula and applying the convolution theorem then yields
$$u(x,t)=(4\pi kt)^{-\frac{1}{2}}\int_{-\infty}^{\infty}\phi(x)e^{-(x-y)^{2}/4kt}\;dy,$$
where $\phi(x)=u(x,0).$ This implies that our fundamental solution is
$$\Phi(x,t)=(4\pi kt)^{-\frac{1}{2}}e^{-x^{2}/4kt}.$$
IV. Boundary conditions and special solutions to the B-S PDE
Let us examine some typical boundary conditions in the Black-Scholes PDE.