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

Question

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?

Answer 1

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[0] == 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[0] == 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}]

Answer 2

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.

Answer 3

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