# how to obtain solution for $u_t+u=u_{xx}$

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

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]


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]


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]


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;  May be next version of Mathematica can solve this? • Is this even general enough to be implemented? This only works for rectangular domains. Feb 9, 2018 at 22:33 • 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... Feb 9, 2018 at 23:49 • @VsevolodA. I am not sure I understand your comment. This is 1D problem. What exactly is the this when you say "this only works"? Feb 10, 2018 at 3:59 • @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. Feb 11, 2018 at 6:05 • @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. Feb 11, 2018 at 13:53

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]


% /. 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