2
$\begingroup$

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}]

enter image description here 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]
$\endgroup$
11
  • $\begingroup$ does this relate your question? $\endgroup$
    – Xminer
    Commented Mar 20, 2023 at 22:56
  • $\begingroup$ It can't solve the general heat pde because it is very hard to find symbolic solution to the eigenvalues to the spatial ode $X(x)$. Numerically it is possible to find eigenvalues, but not symbolically. Try DEigensystem on X[x] ode and you will see that it can't solve it. $\endgroup$
    – Nasser
    Commented Mar 20, 2023 at 23:07
  • $\begingroup$ @Nasser so how can I calculate the analytical solution? $\endgroup$ Commented Mar 20, 2023 at 23:26
  • $\begingroup$ @xzczd How could I involve integrate? $\endgroup$ Commented Mar 21, 2023 at 10:12
  • $\begingroup$ @xzczd I mean How could I write this code? I don't quite understand it. $\endgroup$ Commented Mar 21, 2023 at 10:48

2 Answers 2

1
$\begingroup$

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

enter image description here

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.

$\endgroup$
1
3
$\begingroup$

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]

Mathematica graphics

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)

enter image description here

Update to answer comment on plot

You could do the following

enter image description here

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]

Mathematica graphics

 makePlot[.3]

Mathematica graphics

etc...

$\endgroup$
17
  • $\begingroup$ Numerically I think it is correct what I have done! $\endgroup$ Commented Mar 20, 2023 at 23:43
  • $\begingroup$ @AthanasiosParaskevopoulos Yes, numerical approach is only practical way to solve this. I did not look/verify your numerical code, I assume it is correct. $\endgroup$
    – Nasser
    Commented Mar 20, 2023 at 23:49
  • $\begingroup$ Via Fourier transform what do I have to do? Do I have to write the solution like the Fourier series? $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}$ $\endgroup$ Commented Mar 20, 2023 at 23:52
  • 1
    $\begingroup$ @AthanasiosParaskevopoulos change? Nothing. The code I used is posted below the plot/Manipulate. It works as is on V 13.2.1. DO you mean you do not want to see the evolution of the solution for the times given and have specific time only? The code is there, feel free to change it as you like. $\endgroup$
    – Nasser
    Commented Mar 21, 2023 at 0:24
  • 1
    $\begingroup$ @AthanasiosParaskevopoulos sure. Will add that also. $\endgroup$
    – Nasser
    Commented Mar 21, 2023 at 0:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.