# Solve initial-boundary value problem of heat equation without NDSolve or DSolve

I want to solve the following pde

pde = D[u[t, x], t] - D[u[t, x], x, x] == 0


with the boundary conditions defined as

bc = {u[t, 0] == Sin[t], u[t, 1] == 0, u[0, x] == 0}


However, I just learned the separation of variables in class and this kind of pde is not discussed in my textbook. I am wondering how I can solve it without using DSolve or NDSolve?

• One would expect DSolve[{pde, bc}, u, {x, t}] to provide a symbolic solution, but it returns unevaluated. NDSolveValue[{pde, bc}, u, {x, 0, 1}, {t, 0, 4}] gives a numerical solution. Dec 19, 2016 at 0:09
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Dec 19, 2016 at 2:51
• "I just learned the separation of variables in class and this kind of pde is not discussed in my textbook. " Well, then I should say the introduction for separation of variables in your textbook is too brief. Dec 19, 2016 at 3:17

Here is analytical solution using eigenfunction expansion method.

Clear[x, t, a, b, an, sol, L0, k, qn, uEquib]
pde = Hold[D[u[t, x], t] - D[u[t, x], x, x] == 0];

(* find equilibrium solution, from setting D[u[t,x],t]=0*)
uEquib = a x + b;
uEquib = uEquib /. First@Solve[{(uEquib /. x -> 0) == Sin[t],
(uEquib /. x -> L0 ) == 0}, {a, b}] (*use known eigenfunction for homogenouse heat PDE*)
eigFunction = Sin[n Pi x /L0];
pde = ReleaseHold[pde /. u[t, x] -> (an[t] eigFunction + uEquib)] (*replace the new source by eigenfunction expansion*)
pde = pde /. (D[uEquib, t] -> qn eigFunction) (*find qn by orthogonality*)
qn = 2/L0 Assuming[Element[n, Integers], Integrate[ D[uEquib, t] eigFunction, {x, 0, L0}]] (*initial conditions, since v=u-r and u at t=0 is zero*)
ic = uEquib /. t -> 0

(*solve for an[t]*)
c = an[t] /. First@DSolve[{pde, an == ic}, an[t], t] L0 = 1; (*length*)
sol[t_, x_, k_] := Evaluate[Sum[ c eigFunction + uEquib , {n, 1, k}]] Plot the solution for 15 terms

Plot3D[sol[t, x, 15], {t, 0, 2 Pi}, {x, 0, L0}, AxesLabel -> {t, x, u}] Complete code in one block

Clear[x, t, a, b, an, sol, L0, k, qn, uEquib]
pde = Hold[D[u[t, x], t] - D[u[t, x], x, x] == 0];
(* equilibrium solution, from setting D[u[t,x],t]=0*)
uEquib = a x + b;
uEquib = uEquib /. First@Solve[{(uEquib /. x -> 0) == Sin[t], (uEquib /. x -> L0 ) == 0}, {a, b}];
eigFunction = Sin[n Pi x /L0];
pde = ReleaseHold[pde /. u[t, x] -> (an[t] eigFunction + uEquib)];
(*replace the new source by eigenfunction expansion*)
pde = pde /. (D[uEquib, t] -> qn eigFunction) ;
(*find qn by orthogonality*)
qn = 2/L0 Assuming[Element[n, Integers],Integrate[ D[uEquib, t] eigFunction, {x, 0, L0}]];
(*initial conditions, since v=u-r and u at t=0 is zero*)
ic = uEquib /. t -> 0;
(*solve for an[t]*)
c = an[t] /. First@DSolve[{pde, an == ic}, an[t], t];
L0 = 1;
sol[t_, x_, k_] := Evaluate[Sum[ c eigFunction + uEquib , {n, 1, k}]]
Plot3D[sol[t, x, 15], {t, 0, 2 Pi}, {x, 0, L0}, AxesLabel -> {t, x, u}]


Animation:

    Animate[Plot[sol[t, x, 10], {x, 0, L0},
PlotRange -> {{0, L0}, {-10, 10}}, Frame -> True,
FrameLabel -> {{"u[x,t]", None}, {"x",
Row[{"time = ", NumberForm[N@t, {3, 2}], " seconds"}]}}], {t, 0,
20, 1/10}] Here's a solution based on Laplace transform. First transform the equation and b.c.s and substitute the i.c. into them:

pde = D[u[t, x], t] - D[u[t, x], x, x] == 0;
ic = {u[0, x] == 0};
bc = {u[t, 0] == Sin[t], u[t, 1] == 0};

teqn = LaplaceTransform[{pde, bc}, t, s] /. Rule @@@ ic /.
HoldPattern@LaplaceTransform[a_, __] :> a /. u -> (U[#2] &) Here U[x] denotes Laplace transform of u[t, x], I make this change because DSolve can't handle expression involving LaplaceTransform directly.

The next step is to solve the ODE and obtain the transformed solution:

tsol = DSolve[teqn, U[x], x][[1, 1, -1]]


Oops… I still use DSolve for solving the ODE, but I guess what you mean in the question is "not to solve the PDE with DSolve or NDSolve directly", right?

The last step is to do the inverse Laplace transform, but InverseLaplaceTransform can't handle tsol directly, luckily this can be circumvented by expanding tsol with Fourier sine series first:

easyFourierSinCoefficient[expr_, t_, {a_, b_}, n_] :=
FourierSinCoefficient[expr /. t -> t + a, t, n,
FourierParameters -> {1, Pi/(b - a)}] /. t -> t - a

easyTerm[t_, {a_, b_}, n_] := Sin[Pi/(b - a) n (t - a)]

coe = easyFourierSinCoefficient[tsol, x, {0, 1}, n]

term = easyTerm[x, {0, 1}, n]

solterm = InverseLaplaceTransform[coe term, s, t]


The solution is a summation of solterm from 1 to Infinity:

sol = HoldForm[Sum[#, {n, Infinity}]] &@solterm


$$u(t,x)=\sum _n^{\infty } 2 n \pi \left(\frac{e^{-n^2 \pi ^2 t}}{1+n^4 \pi ^4}+\frac{-\cos (t)+n^2 \pi ^2 \sin (t)}{1+n^4 \pi ^4}\right) \sin (n \pi x)$$

I'm still in v9, you can use Inactivate in or after v10 instead of HoldForm.

With the help of finiteFourierSinTransform, I can now add a solution that doesn't make use of DSolve. We just need to:

1. Make finite Fourier sine transform:

{teq, tic} =
finiteFourierSinTransform[{pde, ic}, {x, 0, 1}, n] /. Rule @@@ bc /.
HoldPattern@finiteFourierSinTransform[a_, __] :> a

2. Make Laplace transform:

tteq = LaplaceTransform[teq, t, s] /. Rule @@@ tic

3. Solve for the transformed solution:

ttsol = Solve[tteq, LaplaceTransform[u[t, x], t, s]][[1, 1, -1]]

4. Transform back:

tsol = InverseLaplaceTransform[ttsol, s, t]
sol = inverseFiniteFourierSinTransform[tsol, n, {x, 0, 1}]


The result is same as above.

As I noted in a comment, one might hope that DSolve[{pde, bc}, u, {x, t}] could handle this problem, but it returns unevaluated. Although this problem can be solved by Laplace transforms, it is more straightforward to solve it numerically.

s = NDSolveValue[{pde, bc}, u, {x, 0, 1}, {t, 0, 2 Pi}];
Plot3D[s[t, x], {t, 0, 2 Pi}, {x, 0, 1}, AxesLabel -> {t, x, u}] • Thanks for the answer. Could you please suggest a way to do Laplace Transform? I just do not want to use NDSolve. Dec 19, 2016 at 2:12
• @Daniel Laplace transform the PDE and x == 0 boundary condition. Then, solve the resulting ODE in x, and perform the inverse Laplace transform. Mathematica has commands for all three steps. If you can wait for about 3 hours, I can take a more thorough look at it. Dec 19, 2016 at 2:18