NDSolve over derivatives of Heaviside function

Question

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.

Answer 1

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!

Answer 2

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

Addition

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.

Answer 3

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

Answer 4

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.