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?

  • 1
    $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this meta Q&A helpful $\endgroup$
    – Michael E2
    Sep 27 '16 at 15:07
  • 1
    $\begingroup$ You have an error in bc2: It should be D[u[0,t], t]==1. In addition, Mathematica will not find an analytic solution for this situation. You should try NDSolve instead. $\endgroup$
    – march
    Sep 27 '16 at 15:47
  • 1
    $\begingroup$ @march I think he wants Derivative[1, 0][u][0, t] as the boundary condition. But you're right about not finding a solution to the PDE even after the correction. $\endgroup$
    – Jens
    Sep 27 '16 at 15:51
  • 1
    $\begingroup$ This doesn't appear to do what I want it will help to show what you obtained as answer, and what you expected. in version 11 your code above gives !Mathematica graphics $\endgroup$
    – Nasser
    Sep 27 '16 at 17:35
  • 1
    $\begingroup$ @Nasser Your solution is valid for bc2 = u[0, t] == 1, but not for the bc2 in the question. $\endgroup$
    – bbgodfrey
    Sep 27 '16 at 19:47

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] *)

enter image description here

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.

  • $\begingroup$ Can you explain why this method was necessary? What is the issue behind trying to use DSolve? $\endgroup$
    – Myridium
    Sep 28 '16 at 4:51
  • $\begingroup$ @Myridium My guess is DSolve is still not strong enough even after the improvement in v10.3. (Here's another example that I think DSolve should handle. ) $\endgroup$
    – xzczd
    Sep 28 '16 at 5:06

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

enter image description here

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.


\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:

A0 = 1; B0 = 25; L0 = 1; k = 1; f = 25; tmax = 1;
En = (2/L0) Assuming[Element[n, Integers] && n > 0,
    Integrate[ (f - (A0 (x - L0) + B0)) Cos[(n - 1/2) Pi/L0 x], {x, 0,
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}];
 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}]


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

Mathematica graphics

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

enter image description here


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.