# Initial-boundary value problem of 2nd order linear PDE with variable coefficient

Today I was studying partial differential equations and I tried to check if my solution was correct, I thought Mathematica code would be easy, but it wasn't. Here is my problem:

Solve $$u_{xx}-2tu-u_t=0$$ for $$0 and $$t>0$$ with the boundary conditions $$u_x(0,t)=u(1/2,t)=0$$ and initial conditions $$u(x,0)=1-2x$$.

By hand I obtained $$u(x,t)=\frac{8}{\pi^2}\sum_{n=0}^{\infty}\frac{1}{(2n+1)^2}e^{-t^2-(2n+1)^2\pi^2t}\cos((2n+1)\pi x)$$

Then I tried with the following code in Mathematica:

pde = D[u[x, t], {x, 2}] - 2 t*u[x, t] == D[u[x, t], {t, 1}];
ic = {u[x, 0] == 1 - 2 x};
bc = {u[1/2, t] == 0, Derivative[1, 0][u][0, t] == 0};
sol = First@(u[x, t] /. DSolve[{pde, ic, bc}, u[x, t], {x, t}]);


I don't know why I am getting error. Thanks in advance for your help.

• The boundary condition doesn't look correct to me. It looks to me like it's saying that the derivative in the x-direction for x = 0 and any time is zero. This means that there must be no change along the x-direction, but your initial condition u[x, 0] = 1- 2x clearly shows that there is a linear slope at x = 0, t = 0, so there must be some change in the x-direction. In fact, if I switch your code to NDSolve, I get an error "Warning: boundary and initial conditions are inconsistent." It looks to be like it deals with it by simply ignoring one of your boundary conditions. Jan 8, 2021 at 6:39
• @MassDefect The problem isn't the inconsistent boundary condition (which is a common feature of many Fourier problems). It's that Mathematica doesn't seem able to solve this particular PDE. If we put pde = D[u[x, t], t] == D[u[x, t], {x, 2}] (heat equation), the code works. Jan 8, 2021 at 7:43
• @yawnoc Ah, okay, I didn’t realize that. Jan 8, 2021 at 8:05
• I tried to check the hand solution by substituting it into the original DE. Mathematica struggled with the infinite sum so I looked at it for a little bit and thought I might be able to justify checking for arbitrary n. u[x_,t_]:=8/Pi^2 1/(2n+1)^2 E^(-t^2-(2n+1)^2 Pi^2 t) Cos[(2n+1)Pi x]; Simplify[D[u[x, t], {x, 2}] - 2 t*u[x, t] == D[u[x, t], {t, 1}]] returns True If I haven't made any mistakes then this might offer some support that the hand solution is correct.
– Bill
Jan 8, 2021 at 8:20
• Maple confirms the solution you found
– rmw
Jan 8, 2021 at 9:55

DSolve can solve the problem symbolically if you help it a bit. Just introduce the transform $$u=e^{-t^2}v$$:

newsys = {pde, ic, bc} /. u -> ({x, t} |-> Exp[-t^2] v[x, t]) // Simplify[#, t >= 0] &

(tsol = v[x, t] /. DSolve[newsys, v[x, t], {x, t}][[1]]) // AbsoluteTiming
(* {31.8537, ……} *)
sol = Exp[-t^2] tsol // Simplify[#, K[1] ∈ Integers] &

sol /. K[1] -> n


$$u(x,t)=4 e^{-t^2} \underset{n=0}{\overset{\infty }{\sum }}\frac{2 e^{-(2 \pi n+\pi )^2 t} \cos ((2 n+1) \pi x)}{(2 \pi n+\pi )^2}$$

Tested in v12.2.

To check your solution perhaps a numerical solution of the PDE helps. Without warning Mathematica is able to solve it if you transform your equation to a first order system:

pde = {D[u[x, t], {x, 1}] == v[x, t],D[v[x, t], {x, 1}] - 2 t*u[x, t] == D[u[x, t], {t, 1}]};
ic = {u[x, 0] == 1 - 2 x };
bc = {u[1/2, t] == 0, v[0, t] == 0};
U = NDSolveValue[Join[{pde, ic, bc}], u , {x, 0, 1/2}, {t, 0, 1}]
Plot3D[U[x, t], {x, 0, 1/2}, {t, 0, 1 }, PlotRange -> All,AxesLabel -> { x, t, None}]


check the solution:

xt = U["Grid"];
xtu = Map[{x, t,Sum[8/Pi^2 1/(2 n +1)^2 Exp[-t^2 - (2 n + 1)^2 Pi^2 t] Cos[(2 n + 1) Pi x], {n, 0, 20}]} /. {x -> #[[1]], t -> #[[2]]} &,
xt ]; // Quiet
Show[{Plot3D[U[x, t], {x, 0, 1/2}, {t, 0, 1 }, PlotRange -> All,AxesLabel -> { x, t, None}], ListPointPlot3D[xtu]}]


• First code dosen't work for me on MMA 12.2.0 ? Works on MMA 12.1.1.0 Jan 8, 2021 at 10:29
• @MariuszIwaniu On Mathematica 11.1 it runs without problems! Jan 8, 2021 at 10:31
• @MariuszIwaniu Sorry, my version is 12.0.0.0 (not 11.1) Jan 8, 2021 at 11:04

An alternative method to solve this problem numerically without the NDSolveValue::ibcinc warning (which appears to discard one of the boundary conditions) is to modify the boundary condition at x = 0 slightly, so that it is consistent with the initial condition there while not changing the computed solution in any noticeable way.

NDSolveValue[{D[u[x, t], {x, 2}] - 2 t*u[x, t] == D[u[x, t], {t, 1}],
u[x, 0] == (1 - 2 x), u[1/2, t] == 0,
Derivative[1, 0][u][0, t] == -2 Exp[-1000 t]},
u[x, t], {x, 0, 1/2}, {t, 0, 1}];
Plot3D[%, {x, 0, 1/2}, {t, 0, 1}, PlotRange -> All, AxesLabel -> {x, t, u},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]


I too find it perplexing that Ulrich Neumann's code does not yield an answer for Mathematica "12.2.0 for Microsoft Windows (64-bit) (December 12, 2020)".