Analytical Solution in Generalized Heat Equation

Question

Given the following problem $$u_{t}-\frac{v}{2}\cdot u_{xx}+\frac{v}{2}\cdot x^2 \cdot u(x)=0 $$ $$u(x,0)=f(x)$$

Where $ f(x)=y_1(x)+0.2y_4(x)+0.01y_6 (x) $

$y_n$ are the Eigenfunctions that defined as $ y_n(x)=\exp\left(\displaystyle{\frac{-x^2}{2}}\right)H_n(x) $

By separation we assume $$u(x,t)=X(x)T(x)$$

I have found that the general solution of $T$ is $T(x)=c_1e^{-\lambda v t}$ and the general solution of $X$ is $ X_n(x)=\exp\left(\displaystyle{\frac{-x^2}{2}}\right)H_n(x) $

So the general solution to our initial problem is $u_n(x,t)=\sum_{n=1}^{\infty}A_n \exp\left(\displaystyle{\frac{-x^2}{2}}\right)H_n(x) c_1e^{-\lambda v t}$

Using the I.C we got $$ u(x,0)=\sum_{n=1}^{\infty}A_n \exp\left(\displaystyle{\frac{-x^2}{2}}\right)H_n(x)=f(x) $$

I found the numerical solution with the following code.

eigen[n_] := n + 1/2
Her[n_] := HermiteH[n, x]
Y[n_] := Exp[-x^2/2]*Her[n]
Y[1]
Y[4]
Y[6]
f[x_] = Y[1] + 0.2 Y[4] + 0.01 Y[6]
IC = u[x, 0] == f[x]
PDE = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[x, t]\)\) - 1/2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(x, x\)]\(u[x, t]\)\) + 
   1/2*x^2*u[x, t] == 0
BCs = {u[-21, t] == 0, u[21, t] == 0}
sol = First[NDSolve[{PDE, BCs, IC}, u, {t, 0, 10}, {x, -21, 21}]]
Plot3D[u[x, t] /. sol, {t, 0, 10}, {x, -21, 21}]

Update

Could anyone help me with the code?

(* Define the Eigenfunctions *)
y[n_, x_] := Exp[-x^2/2] HermiteH[n, x]
(* Define the constant *)
v = 1;
(* Find the general solution for T *)
T[x_, t_] := c1 Exp[-v \[Lambda] t]
(* Find the general solution for X *)
X[n_, x_] := y[n, x]
(* Combine to get the general solution for u *)
u[n_, x_, t_] := X[n, x] T[x, t]
(* Define the initial condition *)
f[x_] := y[1, x] + 0.2 y[4, x] + 0.01 y[6, x]
(* Solve for the coefficients *)
An[n_] := (Integrate[f[x] X[n, x], {x, -\[Infinity], \[Infinity]}]) / (Integrate[X[n, x]^2, {x, -\[Infinity], \[Infinity]}])
(* Define the analytical solution *)
u[x_, t_] := Sum[An[n] X[n, x] T[x, t], {n, 1, \[Infinity]}]
(* Plot the solution for t=0, 0.1, 0.2, 0.3 *)
Plot[{u[x, 0], u[x, 0.1], u[x, 0.2], u[x, 0.3]},{x, -5, 5},PlotRange -> Al]

Answer 1

I solved it!

The issue was that BC must be $\pm \infty$ and the correct normalization needed to be used. Similar to what was done for your generalized wave pde question Difficulty to Solve the Generalized Wave Equation

The analytical solution to the general heat pde in 1D

\begin{align*} u_{t}+\frac{v}{2}x^{2}u & =\frac{v}{2}u_{xx}\\ u\left( x,0\right) & =f\left( x\right) =e^{-\frac{x^{2}}{2}}H_{1}\left( x\right) +\frac{1}{5}e^{-\frac{x^{2}}{2}}H_{4}\left( x\right) +\frac {1}{100}e^{-\frac{x^{2}}{2}}H_{6}\left( x\right) \\ u\left( -\infty,t\right) & =0\\ u\left( +\infty,t\right) & =0 \end{align*}

is the following

$$ u(x,t)=e^{-\frac{x^{2}}{2}}\left( e^{-\frac{v}{2}3t}H_{1}\left( x\right) +\frac{1}{5}e^{-\frac{v}{2}9t}H_{4}\left( x\right) +\frac{1}{100}% e^{-\frac{v}{2}13t}H_{6}\left( x\right) \right) $$

Animation

Code for the above

ClearAll[x, t];
v = 1;
sol[x_, t_] := Exp[-x^2/2]*(Exp[-v/2*3*t]*HermiteH[1, x] + 
   1/5*Exp[-v/2*9*t]*HermiteH[4, x] + 1/100*Exp[-v/2*13*t]*HermiteH[6, x])

Y[n_] := Exp[-x^2/2]*HermiteH[n, x]
ic    = u[x, 0] == Y[1] + 0.2 Y[4] + 0.01 Y[6]
pde   = D[u[x, t], t] + v/2*x^2*u[x, t] == v/2*D[u[x, t], {x, 2}];
L     = 10; (*for NDSolve only*)
bc    = {u[-L, t] == 0, u[L, t] == 0}
nsol  = NDSolveValue[{pde, ic, bc}, u, {x, -L, L}, {t, 0, 3}];

Manipulate[
 Grid[{{Plot[sol[x0, t0], {x0, -10, 10}, PerformanceGoal -> "Quality",
      PlotRange -> {Automatic, {-4, 4}}, ImageSize -> 300, 
     PlotLabel -> "Analytical"]},
   {Plot[Evaluate[nsol[x0, t0]], {x0, -10, 10}, 
     PerformanceGoal -> "Quality", PlotRange -> {Automatic, {-4, 4}}, 
     ImageSize -> 300, PlotLabel -> "Numerical"]
    }}],
 {{t0, 0, "time"}, 0, 3, .01, Appearance -> "Labeled", 
  ContinuousAction -> True}, TrackedSymbols :> {t0}]

Hand solution

\begin{align*} u_{t}+\frac{v}{2}x^{2}u & =\frac{v}{2}u_{xx}\\ u\left( x,0\right) & =f\left( x\right) =e^{-\frac{x^{2}}{2}}H_{1}\left( x\right) +\frac{1}{5}e^{-\frac{x^{2}}{2}}H_{4}\left( x\right) +\frac {1}{100}e^{-\frac{x^{2}}{2}}H_{6}\left( x\right) \\ u\left( -\infty,t\right) & =0\\ u\left( +\infty,t\right) & =0 \end{align*} Let $u=X\left( x\right) T\left( t\right) $, then the pde becomes \begin{align*} T^{\prime}X+\frac{v}{2}x^{2}XT & =\frac{v}{2}X^{\prime\prime}T\\ \frac{T^{\prime}}{T}+\frac{v}{2}x^{2} & =\frac{v}{2}\frac{X^{\prime\prime}% }{X}\\ \frac{2}{v}\frac{T^{\prime}}{T} & =\frac{X^{\prime\prime}}{X}-x^{2}% \end{align*} Hence \begin{align*} \frac{2}{v}\frac{T^{\prime}}{T} & =-\lambda\\ T^{\prime}+\frac{v}{2}\lambda T & =0 \end{align*} The solution to the above is $$ T=Ae^{-\frac{v}{2}\lambda t}% $$ And \begin{align*} \frac{X^{\prime\prime}}{X}-x^{2} & =-\lambda\\ -X^{\prime\prime}+x^{2}X & =\lambda X\\ X\left( -\infty\right) & =0\\ X\left( \infty\right) & =0 \end{align*} The eigenvalues are (see below) $\lambda_{n}=1+2n$ for $n=0,1,2,3,\cdots$. Hence $\lambda_{n}=\left\{ 1,3,5,7,\cdots\right\} $ and corresponding eigenfunction are \begin{align*} X_{n} & =e^{-\frac{x^{2}}{2}}H_{n}\left( x\right) \qquad n=0,1,2,3\cdots\\ \lambda_{n} & =\left\{ 1,3,5,7,\cdots\right\} \end{align*} Hence the pde solution is linear combination of the solutions $X_{n}T_{n}$ or \begin{align} u & =\sum_{n=0}^{\infty}T_{n}X_{n}\tag{1}\\ & =\sum_{n=0}^{\infty}A_{n}e^{-\frac{v}{2}\lambda_{n}t}e^{-\frac{x^{2}}{2}% }H_{n}\left( x\right) \nonumber \end{align} At $t=0$ the above becomes, since $u\left( x,0\right) =f\left( x\right) =e^{-\frac{x^{2}}{2}}H_{1}\left( x\right) +\frac{1}{5}e^{-\frac{x^{2}}{2}% }H_{4}\left( x\right) +\frac{1}{100}e^{-\frac{x^{2}}{2}}H_{6}\left( x\right) $% \begin{align*} e^{-\frac{x^{2}}{2}}H_{1}\left( x\right) +\frac{1}{5}e^{-\frac{x^{2}}{2}% }H_{4}\left( x\right) +\frac{1}{100}e^{-\frac{x^{2}}{2}}H_{6}\left( x\right) & =\sum_{n=0}^{\infty}A_{n}e^{-\frac{x^{2}}{2}}H_{n}\left( x\right) \\ e^{-\frac{x^{2}}{2}}\left( H_{1}\left( x\right) +\frac{1}{5}H_{4}\left( x\right) +\frac{1}{100}H_{6}\left( x\right) \right) & =\sum _{n=0}^{\infty}A_{n}e^{-\frac{x^{2}}{2}}H_{n}\left( x\right) \\ e^{-x^{2}}\left( H_{1}\left( x\right) +\frac{1}{5}H_{4}\left( x\right) +\frac{1}{100}H_{6}\left( x\right) \right) & =\sum_{n=0}^{\infty} A_{n}e^{-x^{2}}H_{n}\left( x\right) \end{align*}

multiplying both sides by arbitrary $H_{m}\left( x\right) $

$$ e^{-x^{2}}H_{m}\left( x\right) \left( H_{1}\left( x\right) +\frac{1}% {5}H_{4}\left( x\right) +\frac{1}{100}H_{6}\left( x\right) \right) =\sum_{n=0}^{\infty}A_{n}e^{-x^{2}}H_{n}\left( x\right) H_{m}\left( x\right) $$

Integrating, and moving the integral into the sum on the RHS gives $$ \int_{-\infty}^{\infty}e^{-x^{2}}H_{m}\left( x\right) H_{1}\left( x\right) dx+\frac{1}{5}\int_{-\infty}^{\infty}e^{-x^{2}}H_{m}\left( x\right) H_{4}\left( x\right) dx+\frac{1}{100}\int_{-\infty}^{\infty}e^{-x^{2}}% H_{m}\left( x\right) H_{6}\left( x\right) dx=\sum_{n=0}^{\infty}A_{n}% \int_{-\infty}^{\infty}e^{-x^{2}}H_{n}\left( x\right) H_{m}\left( x\right) dx $$ By orthogonality $\int_{-\infty}^{\infty}e^{-x^{2}}H_{n}\left( x\right) H_{m}\left( x\right) dx=\sqrt{\pi}2^{m}m!\delta_{nm}$ (see Wikipedia) hence the above becomes

\begin{align*} \int_{-\infty}^{\infty}e^{-x^{2}}H_{m}\left( x\right) H_{1}\left( x\right) dx+\frac{1}{5}\int_{-\infty}^{\infty}e^{-x^{2}}H_{m}\left( x\right) H_{4}\left( x\right) dx+\frac{1}{100}\int_{-\infty}^{\infty}e^{-x^{2}}% H_{m}\left( x\right) H_{6}\left( x\right) dx & =A_{m}\sqrt{\pi}2^{m}m!\\ \sqrt{\pi}2\delta_{n,1}+\frac{1}{5}\left( \sqrt{\pi}2^{4}4!\right) \delta_{n,4}+\frac{1}{100}\left( \sqrt{\pi}2^{6}6!\right) \delta_{n,6} & =A_{n}\sqrt{\pi}2^{n}n! \end{align*}

Hence \begin{align*} A_{n} & =\frac{1}{\sqrt{\pi}2^{n}n!}\left( \sqrt{\pi}2\delta_{n,1}+\frac {1}{5}\left( \sqrt{\pi}2^{4}4!\right) \delta_{n,4}+\frac{1}{100}\left( \sqrt{\pi}2^{6}6!\right) \delta_{n,6}\right) \end{align*}

Therefore the solution (1) becomes \begin{align} u\left( x,t\right) & =\sum_{n=0}^{\infty}A_{n}e^{-\frac{v}{2}\left( 2n+1\right) t}e^{-\frac{x^{2}}{2}}H_{n}\left( x\right) \tag{2}\\ & =A_{1}e^{-\frac{v}{2}\left( 2+1\right) t}e^{-\frac{x^{2}}{2}}H_{1}\left( x\right) +A_{4}e^{-\frac{v}{2}\left( 8+1\right) t}e^{-\frac{x^{2}}{2}} H_{4}\left( x\right) +A_{6}e^{-\frac{v}{2}\left( 12+1\right) t} e^{-\frac{x^{2}}{2}}H_{6}\left( x\right) \\ & =e^{-\frac{v}{2}3t}e^{-\frac{x^{2}}{2}}H_{1}\left( x\right) +\frac{1} {5}e^{-\frac{v}{2}9t}e^{-\frac{x^{2}}{2}}H_{4}\left( x\right) +\frac{1} {100}e^{-\frac{v}{2}13t}e^{-\frac{x^{2}}{2}}H_{6}\left( x\right) \end{align}

Therefore $$ \boxed{ u(x,t)=e^{-\frac{x^{2}}{2}}\left( e^{-\frac{v}{2}3t}H_{1}\left( x\right) +\frac{1}{5}e^{-\frac{v}{2}9t}H_{4}\left( x\right) +\frac{1}{100} e^{-\frac{v}{2}13t}H_{6}\left( x\right) \right) } $$

Attempt with DSolve

ClearAll[x, t, n];
Y[n_] := Exp[-x^2/2]*HermiteH[n, x]
v = 1;
ic = u[x, 0] == Y[1] + 0.2 Y[4] + 0.01 Y[6]
pde = D[u[x, t], t] + v/2*x^2*u[x, t] == v/2*D[u[x, t], {x, 2}];
bc = {u[-Infinity, t] == 0, u[Infinity, t] == 0}
sol = DSolve[{pde, ic, bc}, u[x, t], {x, t}]

No solution.

Answer 2

just to add more to the comment, to help show why this is hard to solve analytically using separation of variables.

To solve using separation of variables we must be able to find the eigenvalues of the BVP ode (the spatial ode) which results. But in this case, since the solution to the ode is very complicated and involves special functions, it can not find the eigenvalues.

\begin{align*} u_{t}+x^{2}u & =u_{xx}\\ u\left( x,0\right) & =f\left( x\right) \\ u\left( 0,t\right) & =0\\ u\left( L,t\right) & =0 \end{align*}

Let $u=X\left( t\right) T\left( t\right) $ then the pde becomes

$$ XT^{\prime}+x^{2}XT-X^{\prime\prime}T=0 $$

Dividing by $XT\neq0$ gives

\begin{align*} \frac{T^{\prime}}{T}+x^{2}-\frac{X^{\prime\prime}}{X} & =0\\ \frac{T^{\prime}}{T} & =\frac{X^{\prime\prime}}{X}-x^{2} \end{align*}

Hence both are equal to some constant, say $-\lambda$. Therefore

\begin{align*} \frac{T^{\prime}}{T} & =-\lambda\\ T^{\prime}+\lambda T & =0 \end{align*}

Which has the solution

$$ T\left( t\right) =c_{1}e^{-\lambda t} $$

The $X\left( x\right) $ ode is now used to determine $\lambda\,$.

\begin{align*} \frac{X^{\prime\prime}}{X}-x^{2} & =-\lambda\\ X\left( 0\right) & =0\\ X\left( L\right) & =0 \end{align*}

And here we see the problem. It is not possible to determine $\lambda_{n}$ values (eigenvalues) which satisfies the BC given and give non zero solution for the eigenfunction $X_{n}\left( x\right) $. The solution to this ode involves ParabolicCylinderD and unresolved integrals. So it is a hopeless case to try to tackle this analytically. DEigensystem can't do it. But numerically it can find eigenvalues and corresponsing eigenfunctions.

 NDEigensystem[{X''[x] - x^2*X[x]}, X[x], {x, 0, 10}, 3]

There might be other approaches to solving this analytically (may be via Fourier transform) but I do not see myself a way via separation of variables method.

Update To answer comment about Transform methods for PDE, I think the best book I found on this is

Transform Methods for Solving Partial Differential Equations (Symbolic and Numeric Computation Series)

Update to answer comment on plot

You could do the following

eigen[n_] := n + 1/2
Her[n_] := HermiteH[n, x]
Y[n_] := Exp[-x^2/2]*Her[n]
Y[1]
Y[4]
Y[6]
f[x_] = Y[1] + 0.2 Y[4] + 0.01 Y[6]
IC = u[x, 0] == f[x]
PDE = \!\(
\*SubscriptBox[\(∂\), \(t\)]\(u[x, t]\)\) - 1/2*\!\(
\*SubscriptBox[\(∂\), \(x, x\)]\(u[x, t]\)\) + 
      1/2*x^2*u[x, t] == 0
BCs = {u[-21, t] == 0, u[21, t] == 0}
sol = First[NDSolve[{PDE, BCs, IC}, u, {t, 0, 10}, {x, -21, 21}]]
Manipulate[
 Plot[u[x, t] /. sol /. t -> t0, {x, -21, 21}, 
  PlotRange -> {Automatic, {-4, 4}}, GridLines -> Automatic, 
  GridLinesStyle -> LightGray, PlotStyle -> Red]
 ,
 {{t0, 0, "time"}, 0, 1, .01, Appearance -> "Labeled"}
 ]

Update to make separate snap shots

eigen[n_] := n + 1/2
Her[n_] := HermiteH[n, x]
Y[n_] := Exp[-x^2/2]*Her[n]
Y[1]
Y[4]
Y[6]
f[x_] = Y[1] + 0.2 Y[4] + 0.01 Y[6]
IC = u[x, 0] == f[x]
PDE = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[x, t]\)\) - 1/2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(x, x\)]\(u[x, t]\)\) + 
         1/2*x^2*u[x, t] == 0
BCs = {u[-21, t] == 0, u[21, t] == 0}
sol = First[NDSolve[{PDE, BCs, IC}, u, {t, 0, 10}, {x, -21, 21}]]

makePlot[t0_?NumericQ] := Module[{},
  Plot[u[x, t] /. sol /. t -> t0, {x, -21, 21}, 
     PlotRange -> {Automatic, {-4, 4}}, GridLines -> Automatic, 
     GridLinesStyle -> LightGray, PlotStyle -> Red, 
   PlotLabel -> Row[{"Solution at t=", t0}]]
  ]

And now

 makePlot[0]

 makePlot[.3]

etc...