# How to implement the finiteness condition for a PDE?

I have a PDE

$$u_t=u_{yy}+u_{zz}$$

subject to the following initial and boundary conditions $$u(y=0,z>0,t>0)=1,$$ $$u(y=0,z<0,t>0)=0,$$ $$u(y=10,z,t>0)=0,$$ $$u(y,z=\pm10,t)\,\, \mbox{is finite,}$$ $$u(y>0,z,t=0)=0.$$

I have solved for $u(y,t)$ using NDsolve but having no clue how to do it for $z$.

Here is my try for $u(y,t)$

ClearAll["Global*"];
pdeset = {Derivative[1, 0][U][t, y] == Derivative[0, 2][U][t, y]}
ics = {U[0, y] == 0};
bcs = {U[t, 0] == 1, U[t, 10] == 0};
bcAll = Flatten[{ics, bcs}, 1];
sol = NDSolve[{pdeset, bcAll}, {U}, {t, 0, 10}, {y, 0, 10}];


Can some body help me to solve for $u(y,z,t)$?

• Are you sure the finiteness condition is set on $z=\pm10$? If so, there will be infinite valid solutions. As far as I know, this condition should be set on the whole $z$ direction. – xzczd Sep 24 '16 at 11:52
• @xzczd It was set on $\infty$. In fact, I have reduced the domain from $(0, \infty)$ to $(0, 10)$. – zhk Sep 24 '16 at 11:55
• Yeah, that's it. This condition can't be translated as you did. An approximate way to translate it is to use $\frac{\partial u}{\partial z}=0$ at $z=±10$, I think. I've solved your problem in another way (which is at least more standard), see the answer below. – xzczd Sep 25 '16 at 15:36

# Analytic(?) Solution

Well, it's a pity that DSolve still can't handle this type of problem even after the improvement in v10.3. (Personally I think it should, because as far as I can tell this problem is essentially a mixture of initial problem and Dirichlet problem for the heat equation! ) So we need to solve it in a more manual way, here's my solution based on integral transform and Fourier series expansion.

Here's the equation set:

eqn = D[u[t, y, z], t] == D[u[t, y, z], y, y] + D[u[t, y, z], z, z];
ic = u[0, y, z] == 0;
bc = {u[t, 0, z] == Piecewise[{{1, z > 0}}], u[t, 10, z] == 0};


Do LaplaceTransform first to eliminate the derivative for t:

ltset = LaplaceTransform[{eqn, bc}, t, τ] /.
HoldPattern@LaplaceTransform[a_, __] :> a /. Rule @@ ic;


Then do FourierTransform to eliminate the derivative for z, here we need to utilize the "shell" here to enhance the FourierTransform:

ftltset = ft[ltset, z, ζ] /. HoldPattern@FourierTransform[a_, __] :> a;


Now the PDE becomes an ODE, solve it:

tsol[τ_, y_, ζ_] = u[t, y, z] /. First@DSolve[ftltset, u[t, y, z], y]
(* -((
E^(-y Sqrt[ζ^2 + τ]) (-E^(20 Sqrt[ζ^2 + τ]) + E^(
2 y Sqrt[ζ^2 + τ])) (I + π ζ DiracDelta[ζ]))/((-1 + E^(
20 Sqrt[ζ^2 + τ])) Sqrt[2 π] ζ τ)) *)


Warning: Don't use Simplify or FullSimplify on this step, because DiracDelta@x x // Simplify returns 0 and this will ruin following calculation and already ruined my whole Sunday.

The last step is to do the inverse transform, but sadly InverseFourierTransform and InverseLaplaceTransform can't handle tsol[τ, y, ζ] directly, we need a work around. Inspired by the structure of solution of Dirichlet problem for the heat equation, we expand tsol with Fourier sine series:

term[n_] =Sin[Pi/10 n y];
coe[n_] = FourierSinCoefficient[tsol[τ, y, ζ], y, n, FourierParameters -> {1, Pi/10}];


The following is a illustration for tsol[τ, y, ζ] and its Fourier sine expansion, the approximation is already not bad with only 50 terms:

Plot[{Total[term@# coe@# &@Range@50] /. {τ -> 1, ζ -> 1}, tsol[1, y, 1]} //
Abs // Evaluate, {y, 0, 10}, PlotRange -> All]


term[n] coe[n] can be InverseFourierTransformed:

lttermsol[τ_, y_, z_, n_] =
Piecewise[{{Simplify[#, z > 0], z > 0}}, Simplify[#, z < 0]] &@
InverseFourierTransform[coe[n] term[n], ζ, z]


Piecewise[{{Simplify[#, z > 0], z > 0}}, Simplify[#, z < 0]] &@ isn't necessary, it just make the expression cleaner and speed up the final calculation a little.

Symbolic inverse Laplace transform seems still impossible, but there exists some packages for numeric inverse Laplace transform. Here I use this one:

numberofterm = 50;(* use more terms if you like *)
coreilt[τ_, y_, z_] := Total@Table[lttermsol[τ, y, z, n], {n, numberofterm}]
sol[t_, y_, z_] := FT[coreilt[#, y, z] &, t]


Solution at $t=10$:

cf2 = Compile[{}, Table[#, {y, 0, 10, 2/5}, {z, -10, 10, 4/5}]] &@sol[10, y, z];
ListPlot3D[Transpose@cf2[], DataRange -> {{0, 10}, {-10, 10}}, PlotRange -> All,
AxesLabel -> {"y", "z", "u"}] // AbsoluteTiming


Animation:

Remark: Though (somewhat luckily) it's not that difficult in this case, the numeric inversion procedure can be really troublesome. See here for an example.

# Numeric Solution

Of course, if the more approximate translation for the finiteness condition mentioned in the comment above i.e. setting $\frac{\partial u}{\partial z}=0$ at $z=\pm10$ is acceptable for you, you can use the following:

nsol=NDSolveValue[{eqn, ic, bc}, u, {t, 0, 10}, {y, 0, 10}, {z, -10, 10},
Method -> {"MethodOfLines", "SpatialDiscretization" -> "FiniteElement"}];

Plot3D[nsol[10,y,z],{y,0,10},{z,-10,10}]


Notice I don't explicitly set the Neumann condition on z direction, and NDSolve will automatically use Neumann zero boundary condition in this case when "SpatialDiscretization" -> "FiniteElement" is set.

• Thank you dear for your efforts. If possible I would really appreciate a numerical solution (NDsolve`) to this problem? – zhk Sep 27 '16 at 2:19
• @mmm See my edit. – xzczd Sep 27 '16 at 2:36