____________________________________________________________________________________
From Mathematica symbolic PDE document.
Solve for \(V(S,t)\) where \(V\) is the price of the option as a function of stock price \(S\) and time \(t\). \(r\) is the risk-free interest rate, and \(\sigma \) is the volatility of the stock.
\[ \frac {\partial V}{\partial t} + \frac {1}{2} \sigma ^2 S^2 \frac {\partial ^2 V}{\partial S^2} = r V - r S \frac {\partial V}{\partial S} \]
With boundary condition \( V(S,T) = \max \{ S-k,0 \}\)
Reference https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_equation See the European call version at bottom of the page.
Mathematica ✓
ClearAll[u, t, x, sigma, k]; pde = D[u[x, t], t] == (1*sigma^2*D[u[x, t], {x, 2}])/2; ic = u[x, 0] == k*Exp[x - 1]*HeavisideTheta[x]; sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, ic}, u[x, t], {x, t}, Assumptions -> k > 0], 60*10]];
\[ \left \{\left \{u(x,t)\to \frac {1}{2} k e^{\frac {\sigma ^2 t}{2}+x-1} \left (\text {Erf}\left (\frac {\sigma ^2 t+x}{\sqrt {2} \sqrt {t} \left | \sigma \right | }\right )+1\right )\right \}\right \} \]
Maple ✓
x:='x';t:='t';u:='u';sigma:='sigma';k:='k'; pde := diff(u(x,t),t) = 1/2*sigma^2*diff(u(x,t),x$2); ic := u(x, 0) = k*exp(x - 1)*Heaviside(x); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde,ic],u(x,t)) assuming k>0),output='realtime'));
\[ u \left ( x,t \right ) =-k{{\rm e}^{-1}} \left ( i{\it invfourier} \left ( {\frac {{{\rm e}^{-1/2\,{s}^{2}{\sigma }^{2}t}}}{s+i}},s,x \right ) -{\it invfourier} \left ( {{\rm e}^{-1/2\,{s}^{2}{\sigma }^{2}t}}{\it fourier} \left ( {{\rm e}^{x}},x,s \right ) ,s,x \right ) \right ) \]
____________________________________________________________________________________
From Mathematica DSolve help pages.
Solve for \(V(t,s)\)
\[ \frac {\partial v}{\partial t} + \frac {1}{2} \sigma ^2 s^2 \frac {\partial ^2 v}{\partial s^2} +(r-q) s \frac {\partial v}{\partial s} - r v(t,s)=0 \]
With boundary condition \( v(T,s) = \psi (s)\)
Reference https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_equation
Mathematica ✓
ClearAll[t, s, v, sigma, psi]; pde = D[v[t, s], t] + (1*sigma^2*s^2*D[v[t, s], {s, 2}])/2 + (r - q)*s*D[v[t, s], s] - r*v[t, s] == 0; bc = v[T, s] == psi[s]; sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, bc}, v[t, s], {t, s}], 60*10]];
\[ \left \{\left \{v(t,s)\to \frac {e^{r (t-T)} \int _{-\infty }^{\infty } \psi \left (e^{K[1]}\right ) \exp \left (-\frac {\left (-K[1]+(T-t) \left (-q+r-\frac {\sigma ^2}{2}\right )+\log (s)\right )^2}{2 \sigma ^2 (T-t)}\right ) \, dK[1]}{\sqrt {2 \pi } \sqrt {\sigma ^2 (T-t)}}\right \}\right \} \]
Maple ✓
t:='t'; s:='s'; sigma:='sigma';v:='v';psi:='psi'; interface(showassumed=0); pde:=diff(v(t, s), t) +s^2*(diff(v(t, s), s, s))/(2*sigma^2)+(r-q)*s*(diff(v(t, s), s))-r*v(t, s) = 0; ic:=v(T, s) = psi(s); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde,ic],v(t,s))),output='realtime'));
\[ v \left ( t,s \right ) =\psi \left ( s \right ) +\sum _{n=1}^{\infty }{\frac { \left ( t-T \right ) ^{n} \left ( U\mapsto rU^{ \left ( n \right ) } \right ) \left ( \psi \left ( s \right ) \right ) }{n!}} \]