3
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey Dec 19 '16 at 0:09
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey Dec 19 '16 at 2:51
  • $\begingroup$ "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. $\endgroup$ – xzczd Dec 19 '16 at 3:17
4
$\begingroup$

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

Mathematica graphics

(*use known eigenfunction for homogenouse heat PDE*)
eigFunction = Sin[n Pi x /L0]; 
pde = ReleaseHold[pde /. u[t, x] -> (an[t] eigFunction + uEquib)]

Mathematica graphics

(*replace the new source by eigenfunction expansion*)
pde = pde /. (D[uEquib, t] -> qn eigFunction)

Mathematica graphics

 (*find qn by orthogonality*)
 qn = 2/L0 Assuming[Element[n, Integers], Integrate[ D[uEquib, t] eigFunction, {x, 0, L0}]]

Mathematica graphics

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

Mathematica graphics

L0 = 1; (*length*)
sol[t_, x_, k_] := Evaluate[Sum[ c eigFunction + uEquib , {n, 1, k}]]

Mathematica graphics

Plot the solution for 15 terms

Plot3D[sol[t, x, 15], {t, 0, 2 Pi}, {x, 0, L0}, AxesLabel -> {t, x, u}]

Mathematica graphics

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

enter image description here

$\endgroup$
3
$\begingroup$

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

Mathematica graphics

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.

$\endgroup$
2
$\begingroup$

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

enter image description here

$\endgroup$
  • $\begingroup$ Thanks for the answer. Could you please suggest a way to do Laplace Transform? I just do not want to use NDSolve. $\endgroup$ – Daniel Dec 19 '16 at 2:12
  • $\begingroup$ @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. $\endgroup$ – bbgodfrey Dec 19 '16 at 2:18

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.