3
$\begingroup$

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)$?

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – xzczd
    Sep 24, 2016 at 11:52
  • $\begingroup$ @xzczd It was set on $\infty$. In fact, I have reduced the domain from $(0, \infty)$ to $(0, 10)$. $\endgroup$
    – zhk
    Sep 24, 2016 at 11:55
  • $\begingroup$ 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. $\endgroup$
    – xzczd
    Sep 25, 2016 at 15:36

1 Answer 1

Reset to default
3
$\begingroup$

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]

Mathematica graphics

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

Mathematica graphics

Animation:

enter image description here

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

enter image description here

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you dear for your efforts. If possible I would really appreciate a numerical solution (NDsolve) to this problem? $\endgroup$
    – zhk
    Sep 27, 2016 at 2:19
  • 1
    $\begingroup$ @mmm See my edit. $\endgroup$
    – xzczd
    Sep 27, 2016 at 2:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.