5
$\begingroup$

I was trying to verify my hand solution for this PDE using Mathematica.

Is a trick to help Mathematica obtain solution to $u_{t}+u=u_{xx}$ with initial conditions $u(x,0)=f(x)$ and boundary conditions $u(0,t)=u(L,t)=0$ ?

I tried changing $L$ to $1$, and tried using explicit $f(x)$, but Mathematica still can't solve it. This is a basic separation of variables PDE which I think Mathematica 11.2 should have been able to solve it. So I thought may be some one can have a trick to help Mathematica with it.

Hand solution

Using separation of variables, let $u\left( x,t\right) =X\left( x\right) T\left( t\right) $. Substituting this back into the PDE gives \begin{align*} T^{\prime}X+TX & =X^{\prime\prime}T\\ \frac{T^{\prime}}{T}+1 & =\frac{X^{\prime\prime}}{X}=-\lambda \end{align*} Where the separation constant is some real value $-\lambda$. This gives the following two ODE's to solve \begin{align} T^{\prime}+(1+\lambda)T & =0\tag{1}\\ X^{\prime\prime}+\lambda X & =0\tag{2} \end{align} Starting with the spatial ODE in order to obtain the eigenvalues. The boundary conditions on the spatial ODE become \begin{align*} X\left( 0\right) & =0\\ X\left( 1\right) & =0 \end{align*}

The above boundary value ODE is standard one and its eigenvalues are $$ \lambda_{n}=\left( \frac{n\pi}{L}\right) ^{2}\qquad n=1,2,3,\cdots $$ The corresponding eigenfunctions are $$ X_{n}\left( x\right) =c_{n}\sin\left( \sqrt{\lambda_{n}}x\right) $$

The solution to the time ODE (1) is, using integrating factor method

$$ T\left( t\right) =e^{-\left( 1+\lambda_{n}\right) t} $$

Therefore \begin{align} u\left( x,t\right) & =\sum_{n=1}^{\infty}u_{n}\nonumber\\ & =\sum_{n=1}^{\infty}T_{n}\left( t\right) X_{n}\left( x\right) \nonumber\\ & =\sum_{n=1}^{\infty}c_{n}e^{-\left( 1+\lambda_{n}\right) t}\sin\left( \sqrt{\lambda_{n}}x\right) \tag{3} \end{align}

In order to determine $c_{n}$, the initial condition is now applied. At $t=0$, $u\left( x,0\right) =f\left( x\right) $ and the above becomes $$ f\left( x\right) =\sum_{n=1}^{\infty}c_{n}\sin\left( \sqrt{\lambda_{n} }x\right) $$ Multiplying both sides by $\sin\left( \sqrt{\lambda_{m}}x\right) $ and integrating over the domain of $f\left( x\right) $ gives $$ \int_{0}^{L}f\left( x\right) \sin\left( \sqrt{\lambda_{m}}x\right) dx=\int_{0}^{L}\sum_{n=1}^{\infty}c_{n}\sin\left( \sqrt{\lambda_{n}}x\right) \sin\left( \sqrt{\lambda_{m}}x\right) dx $$ Interchanging the order of summation and integrating gives $$ \int_{0}^{L}f\left( x\right) \sin\left( \sqrt{\lambda_{m}}x\right) dx=\sum_{n=1}^{\infty}c_{n}\int_{0}^{L}\sin\left( \sqrt{\lambda_{n}}x\right) \sin\left( \sqrt{\lambda_{m}}x\right) dx $$ By orthogonality of $sin$ functions, all terms in the right side vanish except when $n=m$, leading to \begin{align*} \int_{0}^{L}f\left( x\right) \sin\left( \sqrt{\lambda_{m}}x\right) dx & =c_{m}\int_{0}^{L}\sin^{2}\left( \sqrt{\lambda_{m}}x\right) dx\\ & =c_{m}\frac{L}{2} \end{align*} Therefore $$ c_{m}=\frac{2}{L}\int_{0}^{L}f\left( x\right) \sin\left( \sqrt{\lambda_{m} }x\right) dx $$ Since $m$ is an arbitrary number, replacing it back to $n$ \begin{equation} c_{n}=\frac{2}{L}\int_{0}^{L}f\left( x\right) \sin\left( \sqrt{\lambda_{n} }x\right) dx\qquad n=1,2,3,\cdots\tag{4} \end{equation}

But $\sqrt{\lambda_{n}}=\frac{n\pi}{L}$, therefore $$ c_{n}=\frac{2}{L}\int_{0}^{L}f\left( x\right) \sin\left( \frac{n\pi} {L}x\right) dx\qquad n=1,2,3,\cdots $$

The above shows that $c_{n}$ is the Fourier sine series of $f\left( x\right) $. Since $f\left( x\right) $ is not given, explicit solution for $c_{n}$ can not be found. Therefore the final solution is

\begin{align*} u\left( x,t\right) & =\sum_{n=1}^{\infty}c_{n}e^{-\left( 1+\lambda _{n}\right) t}\sin\left( \sqrt{\lambda_{n}}x\right) \\ & =\sum_{n=1}^{\infty}\left( \frac{2}{L}\int_{0}^{L}f\left( x\right) \sin\left( \sqrt{\lambda_{n}}x\right) dx\right) e^{-\left( 1+\lambda _{n}\right) t}\sin\left( \sqrt{\lambda_{n}}x\right) \end{align*}

With $\lambda_{n}=\left( \frac{n\pi}{L}\right) ^{2}$.

Mathematica attempts

ClearAll[x, t, u, f];
pde = D[u[x, t], t] + u[x, t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == f[x];
bc = {u[0, t] == 0, u[L, t] == 0};
DSolve[{pde, ic, bc}, u[x, t], x, t]

Mathematica graphics

ClearAll[x, t, u, f];
pde = D[u[x, t], t] + u[x, t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == x (1 - x);
bc = {u[0, t] == 0, u[L, t] == 0};
DSolve[{pde, ic, bc}, u[x, t], x, t]

Mathematica graphics

ClearAll[x, t, u, f];
pde = D[u[x, t], t] + u[x, t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == x (1 - x);
bc = {u[0, t] == 0, u[1, t] == 0};
DSolve[{pde, ic, bc}, u[x, t], x, t]

Mathematica graphics

Maple 2017.3 is able to solve this

pde:=diff(u(x,t),t)+u(x,t)=diff(u(x,t),x$2);
ic:=u(x,0)=f(x);
bc:=u(0,t)=0,u(L,t)=0;
pdsolve({pde,ic,bc},u(x,t)) assuming L>0;

Mathematica graphics

May be next version of Mathematica can solve this?

$\endgroup$
  • 2
    $\begingroup$ Is this even general enough to be implemented? This only works for rectangular domains. $\endgroup$ – Vsevolod A. Feb 9 '18 at 22:33
  • $\begingroup$ The maple result seems to be hard-coded. You could in principle do the same thing in MMA, but I'm not sure how useful this would be... $\endgroup$ – AccidentalFourierTransform Feb 9 '18 at 23:49
  • $\begingroup$ @VsevolodA. I am not sure I understand your comment. This is 1D problem. What exactly is the this when you say "this only works"? $\endgroup$ – Nasser Feb 10 '18 at 3:59
  • $\begingroup$ @Nasser this is 2D domain you are solving PDE in. And "this" is a spectral method, which can also be applied to 3D Possion equation for example (see FFT poisson solver), but only if boundary is rechtangular. $\endgroup$ – Vsevolod A. Feb 11 '18 at 6:05
  • $\begingroup$ @VsevolodA. You say "this is 2D domain". I really thought it was 1D, since domain is $0<x<L$. May be we are using different definition of 1D and 2D. $\endgroup$ – Nasser Feb 11 '18 at 13:53
6
$\begingroup$

Since DSolve doesn't work well at the moment, you may try finite Fourier transform:

(* Definition of finiteFourierSinTransform and finiteFourierCosTransform
   is not included in this post, please find it in the link above. *)
pde = D[u[x, t], t] + u[x, t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == f[x];
bc = {u[0, t] == 0, u[L, t] == 0};

Format@finiteFourierSinTransform[f_, __] := Subscript[\[ScriptCapitalF], s][f]
Format@finiteFourierCosTransform[f_, __] := Subscript[\[ScriptCapitalF], c][f]

finiteFourierSinTransform[{pde, ic}, {x, 0, L}, n]

Mathematica graphics

% /. Rule @@@ bc
tset = % /. HoldPattern@finiteFourierSinTransform[f_ /; ! FreeQ[f, u], __] :> f;

tsol = DSolve[tset, u[x, t], t][[1, 1, -1]];

sol = inverseFiniteFourierSinTransform[tsol, n, {x, 0, L}] // transformToIntegrate

Mathematica graphics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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