NDSolve over derivatives of Heaviside function

I am trying to numerically solve for the strain of a Maxwell material in response to a step stress. The governing equations are $$\dot{\sigma} + \sigma = \dot{\varepsilon}$$ and I want to find $$\varepsilon$$ with $$\sigma = H(t-1)$$. I have tried

maxwellEqns = {D[sigma[t], t] + sigma[t] == x'[t], x[0] == 0};
maxwellSol = NDSolve[maxwellEqns, x, {t, 0, tf}, Method -> {"PDEDiscretization" -> "FiniteElement"}]


which tells me the PDE is convection dominated. This unfortunately does not reproduce an accurate result as it does not capture the jump in strain at $$t=1$$. I am new to solving PDEs in Mathematica and am looking for advice on how to improve the methods of NDSolve to accurately capture the behavior I am looking for.

• There's some information in this section and the following one in the documentation. FEMDocumentation/tutorial/FiniteElementBestPractice#107053695
– user87932
Commented Oct 17, 2023 at 20:11
• This deals with stabilizing the convection-dominated behavior but still does not lead to the correct physics. Commented Oct 17, 2023 at 21:01
• Well, your equation may contain approximations to the correct physics. What the message may be telling you is that you're neglecting an effect which is actually present in the actual physical situation. Adding artificial diffusion seems to be a more or less standard technique in fields like computational fluid dynamics, but your case may be different.
– user87932
Commented Oct 17, 2023 at 22:16
• The FEM is not well suited for these type of problems, have tried without the option? Commented Oct 17, 2023 at 23:20
• In fact, you deal with HeavisideTheta which is a distribution, not a usual function. As the documetation to DiracDelta says "Numerical routines will typically miss the contributions from measures at single points:". Commented Oct 18, 2023 at 12:14

With MeshOptions->"MaxCellMeasure" -> .1  you get a smooth result and a sharp continuation at t=1.

Try

tf = 3;
sigma[t_] = HeavisideTheta[t - 1]
maxwellEqns = {D[sigma[t], t] + sigma[t] == x'[t], x[0] == 0};
maxwellSol =
NDSolveValue[maxwellEqns, x, {t, 0, tf},
Method -> { "FiniteElement" ,"MeshOptions" -> "MaxCellMeasure" -> .1 } ]

Plot[maxwellSol[t], {t, 0, 3}]


Hope it helps!

• Thanks! Unfortunately the solution should have a discontinuous jump at t=1 Commented Oct 19, 2023 at 22:21
• @JamesVR No, only x'[t] shows this discontiinuos jump at t==1. You asked for solution x[t] Commented Oct 20, 2023 at 6:35
Clear[x]

sigma[t_] = HeavisideTheta[t - 1]

maxwellEqns = D[sigma[t], t] + sigma[t] == x'[t]
(*DiracDelta[t - 1] + HeavisideTheta[t - 1] == x'[t]*)


Integrate the equation to get x[t]

x[t_] = Integrate[DiracDelta[t - 1] + HeavisideTheta[t - 1], t] + c
(*c + (t - 1) HeavisideTheta[t - 1] + HeavisideTheta[t - 1]*)


Solve for the constant c

x[0] == 0
(*c==0*)


So we have

x[t_] = (t - 1) HeavisideTheta[t - 1] + HeavisideTheta[t - 1]

Plot[x[t], {t, 0, 3}]


Numerical solutions rarely, if ever work with a true DiracDelta. It takes an infinitely small step size to capture it. We can use a numerial approximation for the DiracDelta that gets us very close in numerical solutions. One I like to use is

dd[t_] = 100 Sqrt[10/Pi] E^(-100000 t^2)


Its a tall, thin spike whose integral is 1 over the entire t axis as is a true DiracDelta.

we then have for the differential equation:

maxwellEqns = {dd[t - 1] + HeavisideTheta[t - 1] == x'[t], x[0] == 0}

NDSolve[maxwellEqns, x[t], {t, 0, 3}] // Flatten

x[t_] = x[t] /. %

Plot[x[t], {t, 0, 3}]


And the numerical solution looks very close to the analytical one.

• +1. DSolve[{HeavisideTheta[t - 1] + HeavisideTheta'[t - 1] == x'[t], x[0] == 0}, x[t], t] is simpler. Commented Oct 28, 2023 at 9:10
• @user64494 Yes, but a lot less fun. Commented Oct 28, 2023 at 18:38
• (t - 1) HeavisideTheta[t - 1] + HeavisideTheta[t - 1] is a distribution, not a usual function, so its value at t==0 makes no sense. This is an invention of Mathematica developers. Hope, I am clear. Commented Oct 29, 2023 at 7:41
• But it has a limits above an below. @user64494 Commented Oct 29, 2023 at 18:22
• Can you kindly ground your claim? Deep regard. Commented Oct 30, 2023 at 20:06

What is the problem, really?

 sigma[t_]=HeavisideTheta[t-1]
maxwellEqns = {D[sigma[t], t] + sigma[t] == x'[t], x[0] == 0};
maxwellSol =
y = NDSolveValue[maxwellEqns, x, {t, 0, 3},
Method -> {"PDEDiscretization" -> "FiniteElement"}]



The problem is the first continuation over 1 from the right straight line to (0.9,1). The algorithm continues this small second order spline to the left in order to have a zero at x=0.

Look for handling singularity crossing in NDSolve, Methods

• The solution should have a discontinuous jump at t=1 Commented Oct 19, 2023 at 22:21
• It seldomly makes sense to solve a time dependent problem with a spatial solution technique. Commented Oct 28, 2023 at 10:58

Numerically integrating over derivatives of the Heaviside function can be challenging due to the discontinuous nature of the function and its derivatives. The Heaviside function, typically denoted as H(x), is defined as:

H(x) = { 0, for x < 0 1, for x ≥ 0 }

When you calculate the derivative of the Heaviside function, it results in the Dirac delta function (also known as the delta distribution). The Dirac delta function is not a traditional function but rather a distribution that is zero for all x except at x = 0, where it becomes "infinite" in such a way that its integral equals 1.

To numerically integrate over derivatives of the Heaviside function, you would typically need to use a numerical integration method that can handle functions with discontinuities. One common approach is to use a numerical solver that can handle piecewise-defined functions or distributions, like the Dirac delta.

Here's a basic outline of the process:

Discretization: Discretize your integration domain into small intervals.

Evaluate Derivatives: Calculate the derivatives of the Heaviside function over these intervals. For example, if you are working with the first derivative (Dirac delta function), assign appropriate values to these derivatives over the intervals. For the Dirac delta function, it's usually treated as infinite at x = 0 and zero elsewhere, but this needs to be approximated for numerical purposes.

Numerical Integration: Use a numerical integration method suitable for piecewise-defined functions, such as the trapezoidal rule or Simpson's rule. You'll be integrating over these intervals, and the behavior of the derivatives at the points of discontinuity needs to be considered in the integration method.

The specific implementation may depend on the software or programming language you are using for your numerical calculations.

Keep in mind that handling the Dirac delta function or its derivatives numerically can be complex and might require specialized numerical methods or custom algorithms, depending on the specific problem you are trying to solve.

• Your definition of H[x]isn't correct: HeavisideTheta[t] isn't defined for t==0! Commented Oct 18, 2023 at 15:51
• It seems like an answer generated by AI. Commented Oct 28, 2023 at 9:41