# Problem in sharpness of Heaviside Function

The Heaviside function is not sharp, Do not know why?

ss=NDSolve[{D[u[x, t], t] == D[u[x, t], x, x],
u[x, 0] == HeavisideTheta[x - 1] HeavisideTheta[2 - x] ,
u[0, t] == 0, u[Pi, t] == 0}, u, {x, 0, Pi}, {t, 0, 3}]
Plot[u[x, t] /. ss /. t -> 0, {x, 0, Pi}]


• Discretization and interpolation. Jun 17, 2018 at 20:54
• Then how can I invoke sharp Heaviside function in the NDSolve? Jun 18, 2018 at 4:58
• Similar to this post.
– Hugh
Jun 18, 2018 at 17:35
• The more serious problem is, one i.c. Derivative[0,1][u][x,0]=… is missing. Jun 19, 2018 at 6:06
• No, you're not solving diffusion equation. If what you're trying to solve is diffussion equation, then you've made a simple mistake: the left hand side should be D[u[x, t], t] rather than D[u[x, t], t, t]. Jun 22, 2018 at 15:43

In NDSolve you can only approximate a discontinuity in the spatial domain by interpolation over a discrete set. In FEM, this is the only option and the order of interpolation is restricted to 1 or 2. In the method of lines, the time integration may have discontinuities (even if the spatial discretization is done by FEM).

My suggestion is that one can improve the approximation of HeavisideTheta by refining the mesh near x == 1, 2, t == 0. One issue is that HeavisideTheta[0.] is undefined, and therefore x == 1., 2. should be avoided. This is better dealt with by substituting a function for it that is defined everywhere, such as

HeavisideTheta -> UnitStep              (* asymmetric *)
HeavisideTheta -> (1 - UnitStep[-#] &)  (* asymmetric *)
HeavisideTheta -> (1/2 + Sign[#]/2 &)   (* symmetric *)


The question is what value to assign HeavisideTheta[0]: 1, 0, 1/2, etc. The first two are asymmetrical and would result in the sorts of parabolas seen in the OP (if using order 2 interpolation). The value 1/2 would yield a symmetrical result if the grid points were symmetrically placed about the discontinuity. It would result in a straight, steep, but not vertical, line.

Here is one way to construct the symmetrical solution:

Needs@"NDSolveFEM";
With[{dx = 10^-3},
bmesh = ToBoundaryMesh[
"Coordinates" -> {{0, 0}, {1 - dx, 0}, {1 + dx, 0},
{2 - dx, 0}, {2 + dx, 0}, {Pi, 0}, {Pi, 3}, {0, 3}},
"BoundaryElements" -> {LineElement[Partition[Range@8, 2, 1, 1]]},
"MeshOrder" -> 2]
]

emesh = ToElementMesh[bmesh, MaxCellMeasure -> 0.01]
emesh["Wireframe"]


Block[{HeavisideTheta = 1/2 + Sign[#]/2 &},
ss = NDSolve[{D[u[x, t], t, t] == D[u[x, t], x, x],
u[x, 0] == HeavisideTheta[x - 1] HeavisideTheta[2 - x],
u[0, t] == 0, u[Pi, t] == 0}, u, {x, t} ∈ emesh]
]

Plot[u[x, t] /. ss /. t -> 0., {x, 0, Pi}, MaxRecursion -> 15]


A close-up look reveals the lines are not vertical:

Plot[u[x, t] /. ss /. t -> 0., {x, 0.99, 1.01}, MaxRecursion -> 15]


The quality of the solution to the whole PDE seems pretty poor, but that is a separate question.

• Thank you so much, but after adding this emesh, this step function should spread at the subsequent time for Diffusion equation, but /t->1 is not giving the correct behavior Jun 19, 2018 at 8:14
• @SachinKaushik That's what my last sentence meant to hint at, that there's a problem with the mathematics/numerics. See @ andre's answer and @ xzczd's comment. You can combine my trick with method of lines: Use With[{dx = 10.^-3}, With[{c = List /@ Union[Subdivide[0., Pi, 100], {1 - dx, 1 + dx, 2 - dx, 2 + dx}]}, emesh = ToElementMesh["Coordinates" -> c, "MeshElements" -> {LineElement[Partition[Range@Length@c, 2, 1]]}, "MeshOrder" -> 2] ]] and specify {x} \[Element] emesh instead of {x, 0, Pi} in Andre's code, but it still looks wrong. Jun 19, 2018 at 15:00
• the solution you provided is working fine. while editing you make it a wave equation whereas initial conditions and boundary conditions are suitable for diffusion equation (heat equation) and I am finding quality of solution is also good. Only problem is if I add u[x,t] (x-Integrate[x u[x,t],{x,0,Xmax}]) to the diffusion equation in your code with emesh then it shows error, this error is not present if we simply use NDSolve without emesh, however then there is issue of the sharpness of Heaviside function. If you can suggest something how can I invoke this in your code? Jun 20, 2018 at 11:40
• your code is good as long as I have linear equation, for any non linearity it does not work. Jun 20, 2018 at 11:54
• @SachinKaushik FEM is implemented for linear equations only. For nonlinear, NDSolve can offer only the method of lines, which requires a regular grid (or a Chebyshev grid for "Pseudospectral"). That will give you an asymmetric situation leading to a slight overshoot like for your original code, unless you move the unit bump to points that happen to land symmetrically between the regular grid points. With a fine grid spacing, it should still give an decent approximation to the IC, imo. Jun 20, 2018 at 17:44

This is not a answer, only some investigations that are worth to be shared.

Here is the code of the OP, exactly :

ss=NDSolve[{D[u[x, t], t, t] == D[u[x, t], x, x],
u[x, 0] == HeavisideTheta[x - 1] HeavisideTheta[2 - x] ,
u[0, t] == 0, u[Pi, t] == 0}, u, {x, 0, Pi}, {t, 0, 3}]
Plot[u[x, t] /. ss /. t -> 0, {x, 0, Pi}]


The aim of the following approach is to discover :

• the method used by NDSolve

So far I know there are 3 possibilites :

1) Method of Lines + Tensor Product Grid

2) Method of Lines + Finite Element Method

3) Finite Element Method for all the independent variables, including what is generally considered as time (here t) and is treated with the Method of Lines.

• What is the real mesh used by FEM. Is it responsible of the problem of the OP (answer : Yes)

extraction of the interpolatingFunction :

ssFunction=ss[[1,1,2]]


First question : Was FEM used ?

This subject is discuted here

ssFunction["ElementMesh"]


ElementMesh[{{0., 3.14159}, {0., 3.}}, {QuadElement["<" 420 ">"]}]

So FEM was used.

Second Question : Was it a purely finite element method or was the temporal variable treated with the Method of Lines ?

The following code permits to have the answer. It comes from the documentation "finite element Programming":

{state}=NDSolveProcessEquations[{D[u[x, t], t, t] == D[u[x, t], x, x],
u[x, 0] == HeavisideTheta[x - 1] HeavisideTheta[2 - x] ,
u[0, t] == 0, u[Pi, t] == 0}, u, {x, 0, Pi}, {t, 0, 3}]


"SteadyState" indicates that there no temporal variable. If t were considered as a temporal variable, then the code above would have returned :

{NDSolveStateData["<" 0. ">"]}

What is the Mesh used by FEM ?

Since the problem is treated without special consideration for the variable t, the mesh is of dimension 2.
In fact we are mainly interested in the spacial grid at time t=0.

To analyse this, I propose to change the problem from kind 3) to kind 2). (that has the problem of the OP, I mean the bumps, too).
Here is the new code :

ssValue=NDSolve[{D[u[x, t], t, t] == D[u[x, t], x, x],
u[x, 0] == HeavisideTheta[x - 1] HeavisideTheta[2 - x] ,
(D[u[x, t],t] /. t -> 0) == 0,
u[0, t] == 0, u[Pi, t] == 0}, u, {x, 0, Pi}, {t, 0, 3},
Method -> {"MethodOfLines",
"TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement","MeshOptions"->
{"MaxCellMeasure"-> 3. 10^-1}}}]

ssValueFunction=ssValue[[1,1,2]]


Note that I have added (D[u[x, t],t] /. t -> 0) == 0 in your code. This is necessary for the wave equation : see here.
I'm surprised that you can obtain an answer with your code.

Now, one can retrieve the mesh and the coordinates :

Point00=Point[ssValueFunction["ElementMesh"] ["Coordinates"] //
(First /@ # &) //
({#,ssValueFunction[#,0]}& /@ # &)]


Finally, a plot explains what happens in the OP's question :

Plot[ ssValueFunction[x,0.], {x, 0, Pi},Mesh-> All,ImageSize->900,AspectRatio->0.2,
Epilog-> {AbsolutePointSize[8],Red,Point00}]


The red points are the points of the mesh used by the FEM.
It's the interpolation in the function returned by NDSolve that is the origin of the bumps.

• "I'm surprised that you can obtain an answer with your code. " Not that surprising, because "FiniteElement" method has used zero NeumannValue for t==3` and solved the BVP of wave equation, which is a typical ill-posed problem. Jun 19, 2018 at 6:10