# determining solution to PDE with Mathematica

Bug introduced in or after 10.3, persisting through 11.2.

It's been a while since I have touched separation of variable and not too keen to "get my hands dirty" if I may avoid it.

The given PDE is $$\frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2}$$

coupled with the initial and boundary conditions:

$$u_{x}\left ( 0,t \right )-u\left ( 0,t \right )=0$$

$$u\left ( 1,t \right )=0$$

$$u\left ( x,0 \right )=f\left ( x \right )$$

Direct computation with Mathematica deduce the solution to be

DSolve[{D[u[x, t], t] == D[u[x, t], {x, 2}],
D[u[0, t], x] - u[0, t] == 0, D[u[1, t]] == 0, u[x, 0] == f[x]},
u[x, t], {x, t}]


which is not in synced with my "pen and paper" solution.

Would someone be kind enough to verify that Mathematica is giving the right output?

• MMA has no information about f[x] , consequently it gives a solution with the restriction Integrate[f[x] Sin[\[Pi] x K[1]], {x, 0, 1}, Assumptions -> True]. A comparable restriction I would expect in your "pen and paper" (? to secret to show?) solution. – Ulrich Neumann Jan 19 '18 at 8:33
• @Nasser Typological error – Physkid Jan 19 '18 at 9:15
• @Nasser I am going to correct this. – Physkid Jan 19 '18 at 9:47
• Bug still present in M11.3 – user58955 Mar 23 '18 at 11:02

which is not in synced with my "pen and paper" solution.

But looking at it, I could not verify Mathematica solution.

But ignoring this for now, I Just looked at the eigenvalues it used. For simplicity, let $f(x)=x$. This does not affect the eigenvalues. We see

ClearAll[u, x, t]
pde = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == x;
bc = {D[u[0, t], x] - u[0, t] == 0, D[u[1, t]] == 0};
sol = First@DSolve[{pde, ic, bc}, u[x, t], {x, t}];
sol = u[x, t] /. sol /. {K[1] -> n, Infinity -> 20}


So Mathematica is saying the eigenvalues $\lambda_n= (n \pi)^2$ and the eigenfunctions are $\sin(\sqrt \lambda_n x)$. So Mathematica is saying the first few eigenvalues are

 Table[(n Pi)^2,{n,1,6}]//N


But by solving this by hand, the eigenvalues are the solution to

$$\tan(\sqrt \lambda) = \frac{1}{\sqrt \lambda}$$

And this does not have closed form solution. The first few eigenvalues are

   lam/.NSolve[Tan[Sqrt[lam]]==1/Sqrt[lam]&&0<lam<130,lam]


Also, I get the eignfunctions as combination of $\cos$ and $\sin$ and not just $\sin$. So I have no idea how Mathematica obtained this solution. ALso, plotting the numerical solutions vs. analytical, show they are different

(NDSolve gives warning initial and boundary conditions are inconsistent at end points. But this is becuase I picked $f(x)=x$ for initial conditions. This can be ignored for now, it does not affect the overall solution for $t>0$).

ClearAll[u, x, t]
pde = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == x;
bc = {D[u[0, t], x] - u[0, t] == 0, D[u[1, t]] == 0};
sol = First@DSolve[{pde, ic, bc}, u[x, t], {x, t}];
sol = u[x, t] /. sol /. {K[1] -> n, Infinity -> 20};
solN = NDSolve[{pde, ic, bc}, u, {x, 0, 1}, {t, 0, 1}];
Manipulate[
Plot[{Evaluate[Activate[sol] /. t -> t0],
Evaluate[u[x, t] /. solN /. t -> t0]}, {x, 0, 1},
PlotRange -> {Automatic, {0, 1}},
PlotLegends -> {"DSolve", "NDSolve"}],
{{t0, 0.001, "time"}, 0, 1, .01}
]


## Hand solution for first part

Solve \begin{align*} \frac{\partial u}{\partial t} & =\frac{\partial^{2}u}{\partial x^{2}}\\ \frac{\partial u\left( 0,t\right) }{\partial x}-u\left( 0,t\right) & =0\\ \frac{\partial u\left( L,t\right) }{\partial x} & =0\\ u\left( x,0\right) & =f\left( x\right) \end{align*}

Let $u\left( x,t\right) =X\left( x\right) T\left( t\right)$, which gives the $X\left( x\right)$ ODE as

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

Assuming $\lambda<0$, then solution is

$$X\left( x\right) =c_{1}\cosh\left( \sqrt{\lambda}x\right) +c_{2} \sinh\left( \sqrt{\lambda}x\right)$$

Then

$$X^{\prime}\left( x\right) =c_{1}\sqrt{\lambda}\sinh\left( \sqrt{\lambda }x\right) +c_{2}\sqrt{\lambda}\cosh\left( \sqrt{\lambda}x\right)$$

At $x=0$

\begin{align*} 0 & =X^{\prime}\left( 0\right) -X\left( 0\right) \\ & =c_{2}\sqrt{\lambda}-c_{1}\\ c_{1} & =c_{2}\sqrt{\lambda} \end{align*}

Second B.C. gives

\begin{align*} 0 & =X^{\prime}\left( L\right) \\ 0 & =c_{1}\sqrt{\lambda}\sinh\left( \sqrt{\lambda}L\right) \end{align*}

But this mean $c_{1}=0$ since $\sinh\left( \sqrt{\lambda}L\right)$ can't be zero. Therefore $c_{2}$ is also zero. Hence $\lambda<0$ not possible eigenvalue.

Assuming $\lambda=0$ then the solution

\begin{align*} X\left( x\right) & =c_{1}x+c_{2}\\ X^{\prime}\left( x\right) & =c_{1} \end{align*}

First B.C. gives

\begin{align*} 0 & =X^{\prime}\left( 0\right) -X\left( 0\right) \\ & =c_{1}-c_{2}\\ c_{1} & =c_{2} \end{align*}

Second B.C. gives

\begin{align*} 0 & =X^{\prime}\left( L\right) \\ & =c_{1} \end{align*}

Hence $c_{1}=0$ and $c_{2}=0$. Trivial solution. Therefore $\lambda=0$ is not an eigenvalue.

Assuming $\lambda>0$, then solution is

\begin{align*} X\left( x\right) & =c_{1}\cos\left( \sqrt{\lambda}x\right) +c_{2} \sin\left( \sqrt{\lambda}x\right) \\ X^{\prime}\left( x\right) & =-\sqrt{\lambda}c_{1}\sin\left( \sqrt {\lambda}x\right) +\sqrt{\lambda}c_{2}\cos\left( \sqrt{\lambda}x\right) \end{align*}

First B.C. gives

\begin{align} 0 & =X^{\prime}\left( 0\right) -X\left( 0\right) \nonumber\\ & =\sqrt{\lambda}c_{2}-c_{1}\nonumber\\ c_{1} & =\sqrt{\lambda}c_{2} \tag{1} \end{align}

Second B.C. gives

\begin{align} 0 & =X^{\prime}\left( L\right) \nonumber\\ & =-\sqrt{\lambda}c_{1}\sin\left( \sqrt{\lambda}L\right) +\sqrt{\lambda }c_{2}\cos\left( \sqrt{\lambda}L\right) \tag{2} \end{align}

Plugging (1) into (2) gives

\begin{align*} 0 & =-\sqrt{\lambda}\sqrt{\lambda}c_{2}\sin\left( \sqrt{\lambda}L\right) +\sqrt{\lambda}c_{2}\cos\left( \sqrt{\lambda}L\right) \\ 0 & =c_{2}\left( -\lambda\sin\left( \sqrt{\lambda}L\right) +\sqrt {\lambda}\cos\left( \sqrt{\lambda}L\right) \right) \end{align*}

For non-trivial solution, we want $-\lambda\sin\left( \sqrt{\lambda}L\right) +\sqrt{\lambda}\cos\left( \sqrt{\lambda}L\right) =0$ or

\begin{align*} -\lambda\tan\left( \sqrt{\lambda}L\right) +\sqrt{\lambda} & =0\\ \tan\left( \sqrt{\lambda}L\right) & =\frac{1}{\sqrt{\lambda}} \end{align*}

Eigenvalues $\lambda_{n}$ are given by solution to the above trig equation. There is no closed form solution. Hence eigenfunctions are

$$X_{n}\left( x\right) =c_{n}\left( \sqrt{\lambda_{n}}\cos\left( \sqrt{\lambda_{n}}x\right) +\sin\left( \sqrt{\lambda_{n}}x\right) \right)$$

Update

The eigenfuntion (solution to the space boundary value ODE) can be verified using Mathematica 11.2 as follows

ClearAll[y,x,lam];
DSolve[{y''[x]+ lam y[x]==0,y'[0]-y[0]==0,y'[1]==0},y[x],x]


The above shows that the eigenvalues are the solution to

    lam Sin[Sqrt[lam]]==Sqrt[lam] Cos[Sqrt[lam]]


Which agree with the above hand solution.

• There is an unfortunate error in your working. "at x=0:.....", Cosh(0) =1 – Physkid Jan 20 '18 at 3:57
• @Physkid thanks. Corrected. It does not affect the solution. Negative Eigenvalue still remain not valid. – Nasser Jan 20 '18 at 7:15