Result of symbolic definite integral is clearly wrong

The plot of the rather complicated function B(t,d) seems correct (see figure 1).

The plot of its definite integral from 0 to t (in the case d=0) is incorrect (see figure 2). It should have constant value 1 on the interval [1,2], but the plot indicates a constant value of 0 on that interval.

What is going on? Did I misuse the symbolic math functions?

PS Both plots look correct for integers d>1.

syms t k d;
cutoff(t) = piecewise(t < 0, 0, t >= 0, t);
parity(k) = 1-2*mod(k,2);
f(t,k,d) = 1/factorial(d)*parity(k)*nchoosek(d+1, k)*cutoff(t-k)^d;
B(t,d) = symsum(f,k,0,floor(t));
IB(t,d) = int(B(t,d),0,t);

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

bad_d=0;

figure(1)
fplot(B(t,bad_d), [-1,2], 'LineWidth', 2);

figure(2)
fplot(IB(t,bad_d), [-1,2], 'LineWidth', 2);

Answer 1

[Paul]

Hi Dan,

Add assumptions on variables, might help with simplifications.

syms t real
syms d integer
syms k integer
cutoff(t) = piecewise(t < 0, 0, t >= 0, t);
parity(k) = 1-2*mod(k,2);
f(t,k,d) = 1/factorial(d)*parity(k)*nchoosek(d+1, k)*cutoff(t-k)^d;
B(t,d) = simplify(symsum(f,k,0,floor(t)))

B(t, d) = 

bad_d=0;
B(t,bad_d)

ans = 

figure
fplot(B(t,bad_d), [-1,2], 'LineWidth', 2);

●

I'm not sure what fplot is doing for t < 0, because it seems like B(t,0) is not well-defined for that situation.

I know that Matlab allows us to use the same variable as the variable of integration and in the limits of integration, but I always think it's more clear to integrate wrt to a different dummy variable.

syms tau real
IB(t,d) = simplify(int(B(tau,d),tau,0,t))

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

IB(t, d) = 

However, tau is still included in IB(t,d), which is a problem (and is troubling as a result).

Instead of treating d as a parameter in the integral, perhaps it will be better to enter it as a value and then integrate

figure
fplot(int(B(tau,bad_d),tau,0,t),[-1,2])

Is it sensible to fplot the integral for t < 0?

[Dan]

Paul,

Thank you for reading my question (my first on the forum) and for the quick reply.

The final plot of your post is indeed the one that I was expecting to generate (the correct one) in my plot of IB(t, d_bad). I don't really understand how this business with the dummy variable tau works, but it seems to fix the problem in the case d=0 so I tried to adapt your solution. The case d=0 works (thanks!!), but the cases d>0 no longer do (see the last 4 commented lines).

By the way, for the record, the function B(t) is not defined for t<0. I think that I addressed your concern in the modified code below, but this is not the issue at the heart of the matter at hand.

Dan

clear; clc; clear all; close all;

syms t real;
syms k integer;
syms d integer;

cutoff(t) = piecewise(t < 0, 0, t >= 0, t);
parity(k) = 1-2*mod(k,2);
f(t,k,d) = 1/factorial(d)*parity(k)*nchoosek(d+1, k)*cutoff(t-k)^d;
B(t,d) = symsum(f,k,0,floor(t));

syms tau real


% The section below causes a compile error
% IB(t,d) = int(B(tau,d),tau,0,t);
% figure
% fplot(IB(t,0),[0,4], 'LineWidth', 2);
% xlim([-1,4]);

% this works! 
IB0(t) = int(B(tau,0),tau,0,t);
figure(1)
fplot(IB0(t),[0,4], 'LineWidth', 2);
xlim([-1,4]);

% The section below causes a compile error
%IB1(t) = int(B(tau,1),tau,0,t);
%figure(2)
%fplot(IB1(t),[0,4], 'LineWidth', 2);
%xlim([-1,4]);

[Paul]

sympref('default'); % Sometimes needed on Answers
syms t real;
syms k integer;
syms d integer;
cutoff(t) = piecewise(t < 0, 0, t >= 0, t);
parity(k) = 1-2*mod(k,2);
f(t,k,d) = simplify(1/factorial(d)*parity(k)*nchoosek(d+1, k)*cutoff(t-k)^d);
B(t,d) = simplify(symsum(f,k,0,floor(t)))

B(t, d) = 

syms tau real
IB0(t) = int(B(tau,0),tau,0,t)

IB0(t) = 

IB0(t) is an unevaluated integral, though it's not clear to me why int couldn't find a closed form solution, insofar as a closed form solution is found for the general case with d as a parameter.

I think what happens here with fplot is that it might resort to a numerical evaluation using vpaintegral, but I'm not really sure.

%figure(1)
%fplot(IB0(t),[0,4], 'LineWidth', 2);
%xlim([0,4]);
% The section below causes a compile error

With IB1 it's a different story

IB1(t) = int(B(tau,1),tau,0,t)

IB1(t) = 

We have problem because IB1(t) is still a function of tau (tau shows up in the upper limit of the summation), which can't be the case because tau was the variable of integration. Hence, fplot throws an error because it doesn't know what to do with the tau. It migth be tempting to substitute t for tau in IB1(t) and go from there. Hold that thought.

%figure(2)
%fplot(IB1(t),[0,4], 'LineWidth', 2);
%xlim([-1,4]);

Let's try to compute the integral by first computing the anti-derivative of B(t;d)

intB(tau,d) = int(B(tau,d),tau)

intB(tau, d) = 

And then evaluate at the limits of integration

IB(t,d) = intB(t,d) - intB(0,d)

IB(t, d) = 

IB(t,d) = simplify(IB(t,d))

IB(t, d) = 

Note that the upper bound on the sum is floor(t), not floor(tau)

Let's check at d = 0

IB0(t) = IB(t,0)

IB0(t) = 

Our previous expression for IB1(t)

IB1(t)

ans = 

If we assume that int just "forgot" to sub t for tau in the summation, we can do that ourselves (and rename for later use)

originalIB1(t) = simplify(subs(IB1(t),tau,t))

originalIB1(t) = 

●

Compare to our current expression evaluated at d = 1

IB1(t) = IB(t,1)

IB1(t) = 

We see that originalIB1 and IB1 are the same for t >0, except that the former is expressed in terms of nchoosek whereas the latter is not, presumably due to simplifications. We'll come back to this later.

Plot B(t) for a few values of d

figure
fplot(B(t,[0,1,2]),[0,4])

figure
fplot(subplot(3,1,1),IB(t,0),[0,5]);
fplot(subplot(3,1,2),IB(t,1),[0,5]);
fplot(subplot(3,1,3),IB(t,2),[0,5]);
set(get(gcf,'children'),'YLim',[0,1.1])

The plots all cut off at t = d+2. That's because we have a term in the denominator of IB(t) that is factorial(d+1-k), so things are undefined if when k > d+1. k is summed up to floor t, so all is well for d + 1 <= t < d + 2, in which case floor(t) = d + 1, but things are undefined for t >= d+2. Is this expected behavior for this problem?

Returning now to originalIB1 and comparing it to IB1

figure
fplot([IB1(t),originalIB1(t)],[0,5])
ylim([0,1.1])

We see that orignalIB1 doesn't cut off at t = 4. The reason for this is that originalIB1 is still defined in terms of nchoosek. A feature of that function for symbolic inputs is that it returns for zero for cases that ordinarily would be undefined

nchoosek(sym(2),sym(-1))

ans = 

0

nchoosek(sym(2),sym(3))

ans = 

0

This behavior is documented, but I don't know if it's what you want. Maybe allowing simplify to chew on originalIB1 with more steps would get rid of it. In any case, originalIB1 as constituted is happy to be evaluated for cases where IB1 is not.

Now we return to the original approach, which was to not distinguish between the variable of integration and the limits of ingtegration. BTW, when the ingtegrand depends on more than one variable it's always clearer and safer to explicilty specify the variable of integration as done here

IBB(t,d) = simplify(int(B(t,d),t,0,t))

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

IBB(t, d) = 

IBB(t,0)

ans = 

IBB(t,1)

ans = 

figure
fplot(IBB(t,0:2),[0,5])

The plots for d = 1 and d = 2 look the same as above (I didn't overlay). I don't know why the case with d = 0 is not generating a plot.

Frankly, this whole problem is making me a bit nervous. Perhaps there is an alternative approach ....

[Paul]

As an alternative, instead of taking the integral of the sum, try taking the sum of the integrals.

sympref('default'); % Sometimes needed on Answers
syms t real;
syms k integer;
syms d integer;
cutoff(t) = piecewise(t < 0, 0, t >= 0, t);
parity(k) = 1-2*mod(k,2);
f(t,k,d) = simplify(1/factorial(d)*parity(k)*nchoosek(d+1, k)*cutoff(t-k)^d);
syms tau real
IF(t,k,d) = simplify(int(f(tau,k,d),tau,0,t))

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

IF(t, k, d) = 

Not sure what to think of that warning.

Continuing

%B(t,d) = simplify(symsum(f,k,0,floor(t)))
IB(t,d) = symsum(IF(t,k,d),k,0,floor(t))

IB(t, d) = 

IB(t,d) = simplify(IB(t,d))

IB(t, d) = 

figure
fplot(IB(t,[0,1,2]),[0,5])
ylim([0,1.1])

Don't know why IB(t,0) doesn't fplot.

It seems to be well defined

IB(t,0)

ans = 

And evaluates and plot with numerical values

tval = 0:.01:2;
figure
plot(tval,double(IB(tval,0)))

[Paul]

Using

the second approach

redefining cutoff as suggested by @Matt J in this answer

using d+1 for the upper limit of the symsum as suggeste by Matt J in that same answer

and assuming k >= 0 and d >= 0

and NOT simplifying so that the nchoosek term isn't simplified to factorials

sympref('default'); % Sometimes needed on Answers
syms t real;
syms k integer;
syms d integer;
cutoff(t) = piecewise(t < 0, 0, t >= 0, t^d);
parity(k) = 1-2*mod(k,2);
f(t,k,d) = (1/factorial(d)*parity(k)*nchoosek(d+1, k)*cutoff(t-k))

f(t, k, d) = 

syms tau real
assumeAlso([k,d] >= 0)
IF(t,k,d) = (int(f(tau,k,d),tau,0,t))

IF(t, k, d) = 

IB(t,d) = symsum(IF(t,k,d),k,0,d+1)

IB(t, d) = 

figure
fplot(IB(t,[0:4]),[0,10])

Answer 2

[Dan]

The following code does what I wanted. Unfortunately,

1) It's a bit of a hack (why the special treatment for the case d=0?)

2) The plot for the curve d=0 (in blue) has two small wiggles (near t=5 and t=8) that I can't explain.

clear; clc; clear all; close all;

syms t real;
syms k integer;
syms d integer;
syms tau real;  % used as a dummy variable of integration

cutoff(t) = piecewise(t < 0, 0, t);
parity(k) = 1-2*mod(k,2);

f0(t,k,d) = 1/factorial(d)*parity(k)*nchoosek(d+1, k);
f(t,k,d) = 1/factorial(d)*parity(k)*nchoosek(d+1, k)*cutoff(t-k)^d;
B(t,d) = symsum(f,k,0,floor(t));

IB(t,d) = int(B(t,d),t,0,t); % plotting this causes an error when d=0
IB0(t) = int(B(tau,0),tau,0,t);  % plot this function for the case d=0

figure
hold on
fplot(IB0(t),[0,9], 'LineWidth', 2);
fplot(IB(t,1),[0,9]);
fplot(IB(t,2),[0,9]);
fplot(IB(t,4),[0,9]);
fplot(IB(t,8),[0,9]);
xlim([-1,10]);

●