Errors in the symbolic integration of the Dirac Delta function (bug?)

I've been working to numerically integrate some singularity functions and I've encounted an inconsistency in matlab's symbolic integration that occurs when two singularity functions have the same threshold.

To show the issue, I am integrating the following function in three different ways.

In case 1, the function is integrated as written.

In case 2, the two parts of the function are integrated separately and the result is added together.

In case 3, the function is integrated as a whole, but the dirac location is shifted slightly.

Mathematically, all three cases should give exactly the same result.

Even if we accept Matlab's sloppiness with the dirac integration constants (which we shouldn't), the answers should only be off by a constant.

Instead, all three answers are different. Answer #1 is completely wrong as it lacks the dirac function impulse at x=4. Answer #2 and #3 are correct, but different by a constant (which is mathematically valid, but undesirable).

Is there an issue with Matlab or am I missing something in my approach?

%Define variables that will be used symbolically.
syms x c1 c2

%Define a set of q functions
q1=6*heaviside(x-4)+10*dirac(x-4);  %desired function

q2a=6*heaviside(x-4);  %function split into two parts before integration
q2b=10*dirac(x-4);

q3=6*heaviside(x-4)+10*dirac(x-4.000001);  %diraced shift.

%Integrate q
v1=int(q1,x);

v2=int(q2a,x)+int(q2b,x);

v3=int(q3,x);

%Plot for comparison
subplot(3,1,1), fplot(v1,[0,8]),  xlim([0,8]); ylim([-10,25]);
subplot(3,1,2), fplot(v2,[0,8]),  xlim([0,8]); ylim([-10,25]);  hold on, plot([0,8],[0,0],'--r')
subplot(3,1,3), fplot(v3,[0,8]),  xlim([0,8]); ylim([-10,25]);

Answer 1

[Paul]

Hi Geoffrey,

Not sure why there seems to be a problem with the antiderivative v1. You should consider reporting it as bug. I've had trouble in the past with an integral involving dirac and filed a bug report, and the bug was fixed, but haven't done much with antiderivatives.

In the interim, things seems to work as expected for these cases using an integral. Don't know if that will work for other cases you might be looking at.

syms x c1 c2
%Define a set of q functions
q1  = 6*heaviside(x-4) + 10*dirac(x-4);  %desired function
q2a = 6*heaviside(x-4);  %function split into two parts before integration
q2b = 10*dirac(x-4);
q3  = 6*heaviside(x-4) + 10*dirac(x-4.000001);  %diraced shift.

syms t
%Integrate q
%v1 = int(q1,x)
v1a = simplify(int(q1,x,-inf,t),100)

v1a = 

%v2 = int(q2a,x) + int(q2b,x)
v2a = simplify(int(q2a,x,-inf,t) + int(q2b,x,-inf,t))

v2a = 

%v3 = int(q3,x)
v3a = simplify(int(q3,x,-inf,t))

v3a = 

figure
fplot([v1a,v2a,v3a],[0 8])

[Geoffrey Recktenwald]

Paul,

Thank you for the answer and code. It seems that being more explicit in the integral has resolved the problem. This is very helpful.

I've submitted the bug report.

Best,

Geoff

[Sam Chak]

Yup, there seems to be a bug involving the indefinite integral of shifted Dirac delta function.

syms f(x) g(x) h(x) c
assume(c, "real")

%% Case 1: original expression
f(x)= 6*heaviside(x-4) + 10*dirac(x-4);
F   = int(f, x)

F(x) = 

%% Case 2: Dirac delta function is removed
g(x)= 6*heaviside(x-4);
G   = int(g, x)

G(x) = 

%% Case 3: Shift operator is symbolicized as 'c'
h(x)= 6*heaviside(x-c) + 10*dirac(x-c);
H   = int(h, x)

H(x) = 

Hnew= subs(H, c, 4)

Hnew(x) = 

%% Case 4: No shift
f(x)= 6*heaviside(x) + 10*dirac(x);
F   = int(f, x)

F(x) = 

[Paul]

@Geoffrey Recktenwald

I'm curious about this statement: "Even if we accept Matlab's sloppiness with the dirac integration constants (which we shouldn't), ...."

Assuming there was no bug and all three, two-argument calls to int returned a correct antiderivative as documented, why would such be considered sloppy?

[Geoffrey Recktenwald]

The sloppiness is in the integration of the dirac delta function.

It is generally accepted that

Matlab prefers to use the signum notation, so

The difference of course is mathematically permissible because the 1/2 can be carried by the arbitarary constant of integration. However, it doesn't follow standard convention and thus is 'sloppy'.

To add insult to injury, the 1/2 shift doesn't always occur. See the difference between case 2 and case 3, the shift is not consistently applied.

As a side note, including limits does produce the expected +1/2.

syms x t, simplify(int(dirac(x-4),x,-inf,t),100)

ans = 

Three quick follow up points.

1) I did submit a bug ticket, I've been told it is a "known limitation".

2) The use case for my question is using singularity functions to solve beam deflection problems. I'm kinda shocked that Matlab breaks for such a standard engineering problem.

3) Finally, the answer you provided above worked great for the first integral, but failed once I added a constant of integration (as one would expect). However, going back to the regular int(x) worked fine.

syms t x c1 c2
assume(c1 ~= 0)
%Define a set of q functions
q1  = 6*heaviside(x-4) + 10*dirac(x-4);  %desired function

%Integrate q
v1 = simplify(int(q1,x,-inf,t),100);
v  = subs(v1,t,x)+c1;

%Integrate v
m_fails = simplify(int(v,x,-inf,t),100)

m_fails = 

m_works = int(v)

m_works = 

[Walter Roberson]

@Geoffrey Recktenwald

It is interesting that you got floating point output there. Normally you would get rational output.

I have observed that sometimes here on MATLAB Answers that sympref() defaults to weird states instead of the defaults.

syms t x c1 c2
assume(c1 ~= 0)
%Define a set of q functions
q1  = 6*heaviside(x-4) + 10*dirac(x-4);  %desired function

%Integrate q
v1 = simplify(int(q1,x,-inf,t),100);
v  = subs(v1,t,x)+c1;

%Integrate v
m_fails = simplify(int(v,x,-inf,t),100)

m_fails = 

m_works = int(v)

m_works = 

●

[Paul]

1) I too have gotten the "known limitation" response. Don't know why they don't just call it a bug.

3) I can definitely see why Matlab returns that result.

2) All I know about using singularity functions to solve beam deflection problems is what I read in the last five minutes after a quick Google search. Bearing that in mind ...

If Matlab always returns the correct result for the definite integral, maybe an approach to get around the known limitation would be to separate the problem into integrating q1(x) and dealing with the constants of integration manually. For this example

sympref('default');
syms t x c1 c2 real
assumeAlso(c1 ~= 0)
%Define a set of q functions
q1(x)  = 6*heaviside(x-4) + 10*dirac(x-4);  %desired function
%Integrate q
v1(x) = simplify(int(q1(t),t,-inf,x),100)  % could also make lower bound = 0

v1(x) = 

v1(x) = rewrite(v1(x),'heaviside')

v1(x) = 

v1c(x) = v1(x) + c1;

Verify that v1c(x) is the antiderivative of q1(x)

simplify(diff(v1c(x),x))

ans = 

m1(x) = simplify(int(v1(t),t,-inf,x),100)  % could also make lower bound = 0

m1(x) = 

m1(x) = rewrite(m1(x),'heaviside')

m1(x) = 

m1c(x) = m1(x) + c1*x + c2

m1c(x) = 

Verify that m1c(x) is the antiderivative of v1c(x)

expand(v1c(x))

ans = 

simplify(diff(m1c(x),x))

ans = 

Are v1c(x) and m1c(x) correct answers for this problem?

Alternatively, we can start the integration somewhere to the left of the left edge of the beam (I didn't want to start at 0 in case it's possible to have a Dirac at x = 0) and add the constants of integration as we go, which we'd have to do anyway even if int didn't have the known limitation.

v1(x) = rewrite(int(q1(t),t,-1,x),'heaviside') + c1

v1(x) = 

m1(x) = rewrite(int(v1(t),t,-1,x) + c2,'heaviside')

m1(x) = 

Yields the same result as above, albeit with different values of c2 and c1.