Analytic solution for 1D heat equation

Question

I would like to use Mathematica to solve a simple heat equation model analytically.

I have an insulated rod, it's 1 unit long.

• At x = 0, there is a Neumann boundary condition where the temperature gradient is fixed to be 1. This is to simulate constant heat flux.
• At x = 1, there is a Dirichlet boundary condition where the temperature is fixed at 25.

• The initial temperature of the rod is 25.

  heqn = D[u[x, t], t] == D[u[x, t], {x, 2}];
  ic = u[x, 0] == 25;
  bc1 = u[1,t] == 25;
  bc2 = D[u[0,t]]==1;
  DSolve[{heqn, ic, bc1, bc2},u[x,t], {x,t}]
  

This doesn't appear to do what I want, where have I gone wrong?

Answer 1

Here's a analytic solution based on LaplaceTransform, involving an InverseLaplaceTransform.

We first correct the trivial syntax mistake in bc2:

heqn = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == 25;
bc1 = u[1, t] == 25;
bc2 = D[u[x, t], x] == 1 /. x -> 0;

Then do a LaplaceTransform on the equation set and solve the obtained ODE set with DSolve:

tset = LaplaceTransform[{heqn, bc1, bc2}, t, s] /. Rule @@ ic /. 
  HoldPattern@LaplaceTransform[a_, __] :> a

tsol[x_, s_] = u[x, t] /. First@DSolve[tset, u[x, t], x]

Notice I've made a replacement HoldPattern@LaplaceTransform[a_, __] :> a because DSolve has difficulty in handling those LaplaceTransform[__]s, just keep in mind that the u[x, t] inside DSolve is actually LaplaceTransform[u[x, t], t, s].

OK, the final step is to do a InverseLaplaceTransform on tsol[x, s], but sadly, InverseLaplaceTransform gives back the input formula. Anyway, the goal of your question has been achieved, we've obtained an analytic solution, involving an integral. (Don't forget InverseLaplaceTransform is just an sometimes very hard to calculate integral. ) Here's the solution, in traditional math notation:

$$u(x,t)=\frac{1}{2 \pi i}\int _{\gamma -i \infty }^{\gamma +i \infty }\frac{ e^{-\sqrt{s} x} \left(25 \sqrt{s} e^{\sqrt{s} x}+25 \sqrt{s} e^{\sqrt{s} x+2 \sqrt{s}}+e^{2 \sqrt{s} x}-e^{2 \sqrt{s}}\right)}{\left(e^{2 \sqrt{s}}+1\right) s^{3/2}} e^{s t} d s$$

The solution can be used for evaluation, with the help of packages for numeric Laplace transform, I'll use this package as an example:

plotlst = Table[
   ListLinePlot[(Compile[{{x, _Real, 1}}, #] &@FT[tsol[x, #1] &, t])@Range[0, 1, 1/50], 
    DataRange -> {0, 1}, PlotRange -> {24, 25.1}], {t, 10^-3, 1, 1/50}];

plotlst // ListAnimate
(* Export["a.gif", plotlst] *)

Finally I'd like to say something about the failure of direct usage of DSolve. Personally, I think it's just because DSolve is still not strong enough even after the improvement in v10.3. One may think, it's the inconsistent i.c. and b.c. that stops DSolve from solving the problem, but it's not true, because even if we change the bc2 to

bc2 = D[u[x, t], x] == 1 - Exp[-1000 t] /. x -> 0;

DSolve still returns the input.

Answer 2

Here is a full analytical solution derived by hand calculation

$$ u\left( x,t\right) =x+24+\sum_{n=1}^{\infty}\frac{8}{\left( 1-2n\right) ^{2}\pi^{2}}\cos\left( \left( n-\frac{1}{2}\right) \pi x\right) e^{-\left( \left( n-\frac{1}{2}\right) \pi\right) ^{2}t} $$

And compared to Mathematica's above solution by xzczd result, and they agree.

DSolve does not seem to like the non-homogenous Neumann boundary conditions in this problem. The idea to solve this analytically is to break the problem in two stages. Solve for steady state first, which make the PDE an ODE, but now the ODE with nonhomogeneous B.C. is easily solved. Then use the solution at steady state $u_{E}$ to find the difference solution $v\left( x,t\right) =u\left( x,t\right) -u_{E}$ but now this has a homogeneous B.C. which is easily solved by separation of variables method.

Problem:

\begin{align*} \frac{\partial u\left( x,t\right) }{\partial t} & =k\frac{\partial ^{2}u\left( x,t\right) }{\partial x^{2}}\\ \frac{\partial u}{\partial x}\left( 0,t\right) & =A\\ u\left( L,t\right) & =B\\ u\left( x,0\right) & =f\left( x\right) \end{align*} In this problem $L=1,A=1,B=25,f\left( x\right) =25,k=1$ (thermal diffusivity). But the problem is solved in symbols and only at the very end these numerical values are used in animation.

Assuming the solution at steady state is $u_{E}\left( x\right) $, hence $u_{E}\left( x\right) $ is the solution to $\frac{\partial^{2}u\left( x,t\right) }{\partial x^{2}}=0$, or since now it is time independent, it becomes $$ u_{E}^{\prime\prime}=0 $$ With \begin{align*} u_{E}^{\prime}\left( 0\right) & =A\\ u_{E}\left( L\right) & =B \end{align*} This has the solution $$ u_{E}=c_{1}x+c_{2} $$ After applying the boundary conditions, the solution becomes \begin{equation} u_{E}=A\left( x-L\right) +B\tag{1} \end{equation} The difference solution is \begin{equation} v\left( x,t\right) =u\left( x,t\right) -u_{E}\left( x\right) \tag{2} \end{equation} With $v\left( x,t\right) $ now satisfying \begin{align} \frac{\partial v\left( x,t\right) }{\partial t} & =k\frac{\partial ^{2}v\left( x,t\right) }{\partial x^{2}}\tag{3}\\ \frac{\partial v}{\partial x}\left( 0,t\right) & =0\nonumber\\ v\left( L,t\right) & =0\nonumber\\ v\left( x,0\right) & =f\left( x\right) -u_{E}\nonumber \end{align} PDE (3) is solved now by separation of variables. Let $v=X\left( x\right) T\left( t\right) $. Substituting this into (3) gives $$ \frac{1}{k}\frac{T^{\prime}}{T}=\frac{X^{\prime\prime}}{X} $$ These must be equal to same constant, say $-\lambda$, giving two ODE's to solve, with corresponding B.C. \begin{align*} X^{\prime\prime}+\lambda X & =0\\ X^{\prime}\left( 0\right) & =0\\ X\left( L\right) & =0 \end{align*} And $$ T^{\prime}=k\lambda T $$ The solution to the space $X\left( x\right) $ ODE above is solved.

Case $\lambda<0$

The solution in this case is $X=c_{1}\cosh\left( \sqrt{\lambda}x\right) +c_{2}\sinh\left( \sqrt{\lambda}x\right) $. And $X^{\prime}=\sqrt{\lambda }c_{1}\sinh\left( \sqrt{\lambda}x\right) +\sqrt{\lambda}c_{2}\cosh\left( \sqrt{\lambda}x\right) .$ Applying first B.C. gives $0=\sqrt{\lambda}c_{2}$. Hence $c_{2}=0$. The solution becomes $X\left( x\right) =c_{1}\cosh\left( \sqrt{\lambda}x\right) $. Applying second B.C. gives $0=c_{2}\cosh\left( \sqrt{\lambda}L\right) $. But $\cosh$ is never zero, hence $c_{2}=0$ and solution is trivial. Therefore $\lambda<0$ is not an eigenvalue.

Case $\lambda=0$

ODE\ becomes $X^{\prime\prime}=0$. Solution is $X\left( x\right) =c_{1}x+c_{2}$. And $X^{\prime}=c_{1}$. First B.C. gives $0=c_{1}$. Solution becomes $X=c_{2}$. Applying second B.C. gives $0=c_{2}$. Hence $c_{2}=0$ and trivial solution. Therefore $\lambda=0$ is not an eigenvalue.

Case $\lambda>0$

Solution is $X=c_{1}\cos\left( \sqrt{\lambda}x\right) +c_{2}\sin\left( \sqrt{\lambda}x\right) $. And $X^{\prime}=-\sqrt{\lambda}c_{1}\sin\left( \sqrt{\lambda}x\right) +\sqrt{\lambda}c_{2}\cos\left( \sqrt{\lambda }x\right) $. First B.C. gives $0=-\sqrt{\lambda}c_{2}$ or $c_{2}=0$. Solution becomes $X=c_{1}\cos\left( \sqrt{\lambda}x\right) $. Second B.C. gives $0=c_{2}\cos\left( \sqrt{\lambda}L\right) $. Non trivial solution implies $\cos\left( \sqrt{\lambda}L\right) =0$ or \begin{align*} \sqrt{\lambda_{n}} & =\left( n-\frac{1}{2}\right) \frac{\pi}{L}\qquad n=1,2,3,\cdots\\ \lambda_{n} & =\left( \left( n-\frac{1}{2}\right) \frac{\pi}{L}\right) ^{2}\qquad n=1,2,3,\cdots \end{align*} Giving the eigenfunctions $$ X\left( x\right) =\sum_{n=1}^{\infty}C_{n}\cos\left( \sqrt{\lambda_{n} }x\right) $$ The corresponding time solution for same set of eigenvalues is $$ T\left( t\right) =\sum_{n=1}^{\infty}D_{n}e^{-\lambda_{n}kt} $$ Therefore \begin{align} v\left( x,t\right) & =XT\nonumber\\ & =\sum_{n=1}^{\infty}E_{n}\cos\left( \sqrt{\lambda_{n}}x\right) e^{-\lambda_{n}kt}\tag{4} \end{align} Where $E_{n}=C_{n}D_{n}$. Now $E_{n}$ is found from initial conditions. From (3) $v\left( x,0\right) =f\left( x\right) -u_{E}$, therefore $$ f\left( x\right) -u_{E}=\sum_{n=1}^{\infty}E_{n}\cos\left( \left( n-\frac{1}{2}\right) \frac{\pi}{L}x\right) $$ Multiplying both sides by $\cos\left( \left( m-\frac{1}{2}\right) \frac {\pi}{L}x\right) $ and integrating gives $$ \int_{0}^{L}\left( f\left( x\right) -u_{E}\right) \cos\left( \left( m-\frac{1}{2}\right) \frac{\pi}{L}x\right) dx=\int_{0}^{L}\cos\left( \left( m-\frac{1}{2}\right) \frac{\pi}{L}x\right) \sum_{n=1}^{\infty} E_{n}\cos\left( \left( n-\frac{1}{2}\right) \frac{\pi}{L}x\right) dx $$ Exchanging order of integration and summation $$ \int_{0}^{L}\left( f\left( x\right) -u_{E}\right) \cos\left( \left( m-\frac{1}{2}\right) \frac{\pi}{L}x\right) dx=\sum_{n=1}^{\infty}E_{n} \int_{0}^{L}\cos\left( \left( m-\frac{1}{2}\right) \frac{\pi}{L}x\right) \cos\left( \left( n-\frac{1}{2}\right) \frac{\pi}{L}x\right) dx $$ By orthogonality, $\int_{0}^{L}\cos\left( \left( m-\frac{1}{2}\right) \frac{\pi}{L}x\right) \cos\left( \left( n-\frac{1}{2}\right) \frac{\pi} {L}x\right) dx=\frac{L}{2}$ for $m=n$ and zero otherwise, and the above simplifies to $$ \int_{0}^{L}\left( f\left( x\right) -u_{E}\right) \cos\left( \sqrt{\lambda_{n}}x\right) dx=\frac{L}{2}E_{m} $$ Solving for $B_{m}$ and renaming to $n$ gives \begin{align*} E_{n} & =\frac{2}{L}\int_{0}^{L}\left( f\left( x\right) -u_{E}\right) \cos\left( \sqrt{\lambda_{n}}x\right) dx\\ & =\frac{2}{L}\int_{0}^{L}\left( f\left( x\right) -\left( A\left( x-L\right) +B\right) \right) \cos\left( \sqrt{\lambda_{n}}x\right) dx \end{align*} Hence the solution (4) becomes $$ v\left( x,t\right) =\sum_{n=1}^{\infty}E_{n}\cos\left( \sqrt{\lambda_{n} }x\right) e^{-\lambda_{n}kt} $$ Therefore from $u\left( x,t\right) =u_{E}\left( x\right) +v\left( x,t\right) $ the final solution is \begin{align} u\left( x,t\right) & =u_{E}\left( x\right) +\sum_{n=1}^{\infty}E_{n} \cos\left( \sqrt{\lambda_{n}}x\right) e^{-\lambda_{n}kt}\tag{5}\\ & =A\left( x-L\right) +B+\sum_{n=1}^{\infty}E_{n}\cos\left( \sqrt {\lambda_{n}}x\right) e^{-\lambda_{n}kt}\nonumber \end{align} Now for $A=1,B=25,f\left( x\right) =25$,$k=1,L=1$ then \begin{align*} E_{n} & =\frac{2}{L}\int_{0}^{L}\left( f\left( x\right) -\left( A\left( x-L\right) +B\right) \right) \cos\left( \left( n-\frac{1}{2}\right) \frac{\pi}{L}x\right) dx\\ & =2\int_{0}^{1}\left( 1-\left( \left( x-1\right) +25\right) \right) \cos\left( \left( n-\frac{1}{2}\right) \frac{\pi}{L}x\right) dx\\ & =\frac{8}{\left( 1-2n\right) ^{2}\pi^{2}} \end{align*} Hence the solution is $$ u\left( x,t\right) =x+24+\sum_{n=1}^{\infty}\frac{8}{\left( 1-2n\right) ^{2}\pi^{2}}\cos\left( \left( n-\frac{1}{2}\right) \pi x\right) e^{-\left( \left( n-\frac{1}{2}\right) \pi\right) ^{2}t} $$

Here is the Manipulate code:

ClearAll[u]
A0 = 1; B0 = 25; L0 = 1; k = 1; f = 25; tmax = 1;
En = (2/L0) Assuming[Element[n, Integers] && n > 0,
   Simplify[
    Integrate[ (f - (A0 (x - L0) + B0)) Cos[(n - 1/2) Pi/L0 x], {x, 0,
       L0}]]]
uME[x_, t_, k_, A0_, B0_, L0_] := 
  A0 (x - L0) + B0 + 
   Sum[En  Cos[(n - 1/2) Pi/L0 x] Exp[-((n - 1/2) Pi/L0)^2 k t], {n, 
     1, 50}];
Manipulate[
 Plot[uME[x, t, k, A0, B0, L0], {x, 0, L0}, AxesOrigin -> {0, 0}, 
  PlotRange -> {All, {23.8, 25.4}}, Frame -> True, 
  FrameLabel -> {{"u(t)", None}, {"x", Row[{"time ", t}]}}, 
  GridLines -> Automatic, GridLinesStyle -> LightGray, 
  PlotStyle -> Red, BaseStyle -> 14],
 {t, 0, 1, .01}]

Update

I wanted to add this before but had no time. This below is the analytical solution to this problem using DSolve. By breaking the problem as above, into steady state and transient parts. This is temporary solution for solving such problems with DSolve until it can add support for non homogeneous Neumann boundary conditions (version 12 may be?), which seems the reason why it can not solve it in one shot.

A0 = 1; B0 = 25; L0 = 1; k = 1; f = 25; tmax = 1;
(*solve steady state*)
uE = u[x] /. 
   First@DSolve[{u''[x] == 0, u[L0] == B0, u'[0] == A0}, u[x], x];

(*solve trasnient state*)
ic = v[x, 0] == f - uE;
bc1 = v[L0, t] == 0;
bc2 = Derivative[1, 0][v][0, t] == 0;
sol = v[x, t] /. 
   First@DSolve[{ D[v[x, t], t] == k D[v[x, t], {x, 2}], ic, bc1, 
      bc2}, v[x, t], {x, t}];
sol = sol /. {K[1] -> n, Infinity -> 40}; (*40 terms seems enough*)

(*add solutions*)
sol = sol + uE

Here is animation of the above

Animate[Plot[Activate[sol] /. t -> i, {x, 0, L0}, 
  PlotRange -> {All, {23.8, 25.4}}, Frame -> True, 
  FrameLabel -> {{"u(t)", None}, {"x", Row[{"time ", i}]}}, 
  GridLines -> Automatic, GridLinesStyle -> LightGray, 
  PlotStyle -> Red, BaseStyle -> 14], {i, 0, 1, .01}]

Answer 3

In version 13 on Windows 10

heqn = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == 25;bc1 = u[1, t] == 25;
bc2 = D[u[x, t], x] == 1 /. x -> 0;
DSolve[{heqn, ic, bc1, bc2}, u[x, t], {x, t}]

$$\left\{\left\{u(x,t)\to \underset{K[1]=1}{\overset{\infty }{\sum }}\frac{8 e^{-\frac{1}{4} \pi ^2 t (2 K[1]-1)^2} \cos \left(\frac{1}{2} \pi x (2 K[1]-1)\right)}{\pi ^2 (1-2 K[1])^2}+x+24\right\}\right\} $$

{{u[x, t] -> 24 + x + Inactive[Sum][(8 E^(-(1/4) π^2 t (-1 + 2 K[1])^2) 
        Cos[1/2 π x (-1 + 2 K[1])])/(π^2 (1 - 2 K[1])^2), {K[1], 1, ∞}]}}