Wierd wiggle in plot of symbolic expression

Apart from the small wiggles (the biggest near t=8) the plot of the symbolic expression below looks correct. How can I remove the wiggles?

clear; clc; close all;
syms t real;
syms k integer;
syms d integer;

cutoff(t) = piecewise(t < 0, 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;  % dummy variable of integration just for the case d=0
IB0(t) = int(B(tau,0),tau,0,t); 

figure
hold on
fplot(IB0(t),[0,9], 'LineWidth', 2);
xlim([-1,10]);
ylim([-0.2,1.2]);

[Dan]

Although this code uses Matlab's symbolic functionality, it seems (according to a comment of Matt J) that that some floating point computations are being performed under the hood and inaccuracy in these computations are causing the wiggles.

(In that case, I would guess that the function cutoff(t-k)^d might be the culprit.)

It was suggested that I use a different (more stable) algorithm to perform the computation, but I really only interested in the cause of the wiggles. Specifically, I'd like to know that I didn't misuse the symbolic functions.

Thanks Matt.

Answer 1

[Torsten]

Using

IB0(t) = vpaintegral(@(tau)B(tau,0),tau,0,t,'Waypoints',[1]); 

instead of

IB0(t) = int(B(tau,0),tau,0,t);

makes the wiggles disappear.

The reason for the wiggles is that your function B is discontinuous at tau = 1. Using the points of discontinuity as fixed grid points (Waypoints) in the integration gives exact results.

[Dan]

This works like a charm!

The problem indeed stemmed from my incorrect Matlab code.

I have to admit, however, that I don't (yet!) understand how Torsten's solution works. I guess I'll have to learn more about 'Waypoints'.

Thank you Torsten!

[Torsten]

The problem indeed stemmed from my incorrect Matlab code.

Your MATLAB code is not incorrect, but MATLAB does not return an analytical expression for the function

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

Thus when you try to plot the function IB0(t), MATLAB has to use a numerical method to approximate IB0(t). This numerical method is based on the evaluation of B on a grid of tau-values. If the points of discontinuity of B are not part of these grid values for tau where B is evaluated, results can become slightly imprecise. One way to overcome this is to force the grid to contain the tau-values where B is discontinuous using the "Waypoints" option.

[Dan]

I understand. Thank you again for the clear explanation. It's even more useful than your solution!

Answer 2

[Matt J]

Seems like a lot of unnecessary symbolic math gymnastics for such a simple function.

IB0=@(t) min(t,1);

fplot(IB0,[0,9], 'LineWidth', 2);
xlim([-1,10]);
ylim([-0.2,1.2]);

ylim([-0.2,1.2])

●

[Dan]

It does seem that way, but this is only a particular case (d=0) of a more sophisticated function. I believe that all the components of the symbolic expression are needed in the general case.

One might suggest, for the purposes of generating a plot, that your simplified expression suffices. Maybe so, but I'm interested in understanding the reason for the wiggles. I feel that I must have misformulated the symbolic expression (or the plot command).

By the way, the plots (not shown here) for the cases d>0 look just fine.

[Matt J]

It looks like you are trying to reinvent spline interpolation. That is unwise, since there are already stable tools in Matlab for doing spline interpolation, as well as for integrating splines (e.g. fnint).

Regardless, the truncated power spline basis that you are pursuing is known to be unstable numerically. The gold standard is De Boor's algorithm.

[Dan]

Do you think that this "numerical instability" causes the wiggles?

[Matt J]

Not in this case. The wiggles can be removed here by fixing some errors in your symbolic expressions (see my other answer). However, I expect you could run into numerical problems for sufficiently large d, one of several reasons not to do spline analysis this way.

Answer 3

[Matt J]

There appears to be a mistake in your expression for cutoff(). Also, using floor(t) for the upper limit of the symsum makes it harder on the symbolic engine than it needs to be.

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

cutoff(t,d) = piecewise(t < 0, 0, t^d);  %<--- Matt J fixed

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(t,k,d) ,k,0,d+1); %<--- Matt J use d+1 as upper limit of sum

syms tau real;  % dummy variable of integration just for the case d=0
IB0(t) = int(B(tau,0),tau,0,t); 

figure
hold on
fplot(IB0,[0,9], 'LineWidth', 2);
xlim([-1,10]);
ylim([-0.2,1.2]);

●

[Dan]

Matt,

Thank you for pursuing this issue. I'm definitely learning from it. Still, some mysteries remain.

Firstly, I now realize that the definition of B(t,d) in my original post was simply incorrect. I suspect that you caught on to that right away. Not much else to say about it.

The code below generates four graphs. Mathematically, as far as I can tell, they should all be identical (and correct). Obviously, they're not.

The graphs in the second and fourth figures look correct.

The graph in the third figure looks correct, except for the wiggle. As explained by Torsten in a comment to this post, the difference has to do with the implementation that Matlab employs for a certain numerical routine that is employed by the fplot function. The graph in figure 4 employs a fix that Torsten recommended.

The graph in the first figure certainly looks very wrong, but the only difference between the code used to specify the first figure and the second figure is that between the cutoff1 and cutoff2 functions. These are mathematically identical, so it must be the case that Matlab is treating them differently under the hood. You probably already know about this because you pointed out that I was mistaken in my specification of that function in your latest post.

In the hopes of avoiding a similar pitfall in the future, I'd like to know:

Why does Matlab treat the cutoff1 and cutoff2 so differently?

clear; clc; clear; close all;

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

cutoff1(t,d) = (piecewise(t < 0, 0, t))^d;
cutoff2(t,d) = piecewise(t < 0, 0, t^d);

parity(k) = 1-2*mod(k,2);

f(t,k,d) = 1/factorial(d)*parity(k)*nchoosek(d+1, k);

f1(t,k,d) = f(t,k,d)*cutoff1(t-k,d);
B1(t,d) = symsum(f1,k,0,d+1);

f2(t,k,d) = f(t,k,d)*cutoff2(t-k,d);
B2(t,d) = symsum(f2,k,0,d+1);

f3(t,k,d) = f(t,k,d)*(t-k)^d;
B3(t,d) = symsum(f3,k,0,floor(t));

IB0_1(t) = int(B1(tau,0),tau,0,t);
IB0_2(t) = int(B2(tau,0),tau,0,t);
IB0_3(t) = int(B3(tau,0),tau,0,t);
IB0_4(t) = vpaintegral(@(tau)B3(tau,0),tau,0,t,'Waypoints',[1]); 

figure
fplot(IB0_1(t),[0,9]);
xlim([-1,10]);
ylim([-0.2,1.2]);

●

figure
fplot(IB0_2(t),[0,9]);
xlim([-1,10]);
ylim([-0.2,1.2]);

figure
fplot(IB0_3(t),[0,9]);
xlim([-1,10]);
ylim([-0.2,1.2]);

figure
fplot(IB0_4(t),[0,9]);
xlim([-1,10]);
ylim([-0.2,1.2]);

[Matt J]

These are mathematically identical, so it must be the case that Matlab is treating them differently under the hood

No, they are mathematically identical only for d>0. For d=0, cutoff1(t,d) equals 1 everywhere, whereas cutoff2(t,d) is the unit step that you really want.

syms t real;
syms k integer;
syms d integer;
syms tau real;
cutoff1(t,d) = (piecewise(t < 0, 0, t))^d;
cutoff2(t,d) = piecewise(t < 0, 0, t^d);


fplot(cutoff1(t,0),[-5,5],'r-'); hold on
fplot(cutoff2(t,0),[-5,5],'bo--'); hold off ;axis padded
legend('cutoff1', 'cutoff2', Location='southeast')

●

[Paul]

cutoff1 and cutoff2 are not the same in the sense that they have different effects on f1 and f2 when shifted by k

syms t real;
syms k integer;
syms d integer;
syms tau real;
cutoff1(t,d) = (piecewise(t < 0, 0, t))^d;
cutoff2(t,d) = piecewise(t < 0, 0, t^d);
parity(k) = 1-2*mod(k,2);
f(t,k,d) = 1/factorial(d)*parity(k)*nchoosek(d+1, k);
f1(t,k,d) = f(t,k,d)*cutoff1(t-k,d)

f1(t, k, d) = 

B1(t,d) = symsum(f1,k,0,d+1)

B1(t, d) = 

f2(t,k,d) = f(t,k,d)*cutoff2(t-k,d);
B2(t,d) = symsum(f2,k,0,d+1)

B2(t, d) = 

The first term in B1(tau,d) has 0^d, which is not 0 when d = 0

sym(0)^sym(0)

ans = 

1

Hence, when evaluating B1(tau,0) we get 0

B1(tau,0)

ans = 

0

Which is not the case for B2(tau,0)

B2(tau,0)

ans = 

Great catch by @Matt J, assuming cutoff2 is really the intended function.

And great insight as well to cap the sum at d+1 instead of floor(t)

[Dan]

No, they are mathematically identical only for d>0. For d=0, cutoff1(t,d) equals 1 everywhere, whereas cutoff2(t,d) is the unit step that you really want.

Yes. I see now.

Thanks again for all the insightful comments.

In summary:

cutoff1 is not the same as cutoff2 and the latter should be used to specify the "correct" function, IB0_2.

The expression in IB0_4 is also valid mathematically, but requires a special Waypoint-adjusted specification of the definite integral in order to avoid wiggle artifacts in the Matlab plot.

[Matt J]

The expression in IB0_4 is also valid mathematically, but requires a special Waypoint-adjusted specification of the definite integral in order to avoid wiggle artifacts in the Matlab plot.

Valid in its current form for d=0. For general d, I believe that you would need,

f3(t,k,d) = f(t,k,d)*cutoff2(t-k,d);

B3(t,d) = symsum(f3,k,0,floor(t));

and I believe you would need waypoints at all k=1,2,...,d+1.

[Paul]

Hi Matt,

Why would floor(t) be needed for general d? I think the expression for IB(t,d) in this comment is same using either d+1 or floor(t).

[Matt J]

@Paul

floor(t) is not needed, and note that my answer does not propose it. It is being considered here only because @Dan was seeking to compare my answer to several other proposals.

Avoiding floor(t) is better, because then the symbolic engine is then able to do the integral analytically, rather than numerically. However, this then prevents us from exploring how to accomplish the thing with vpaintegral.

[Dan]

If floor(t) is used, then use of the cutoff function can be avoided. There are advantages to this. For the Matlab implementation, it seems to me that Matt J's solution is the cleanest and quickest for Matlab. His solution avoids the floor(t) function but uses the cutoff function.

For instructional purposes, the following code generates two different figures in three different ways. Perhaps one solution appeals to you more than the other two.

clear; clc; close all;
syms t real;
syms k integer;
syms d integer;
syms tau real;  % dummy variable of integration

parity(k) = (-1)^k;%1-2*mod(k,2);%
f(t,d,k) = 1/factorial(d)*parity(k)*nchoosek(d+1,k);

f1(t,d,k) = f(t,d,k)*(t-k)^d;
B1(t,d) = symsum(f1,k,0,floor(t));
W = [1]; %W = [1,2];    % Waypoints (How to choose these or to know when to use Waypoints in general?)
IB1_0(t) = vpaintegral(@(tau)B1(tau,0),tau,0,t,'Waypoints',W);
IB1_1(t) = vpaintegral(@(tau)B1(tau,1),tau,0,t,'Waypoints',W);
IB1_2(t) = vpaintegral(@(tau)B1(tau,2),tau,0,t,'Waypoints',W);
IB1_4(t) = vpaintegral(@(tau)B1(tau,4),tau,0,t,'Waypoints',W);
IB1_8(t) = vpaintegral(@(tau)B1(tau,8),tau,0,t,'Waypoints',W);

cutoff(t) = piecewise(t<0,0,t);
f2(t,d,k) = f(t,d,k)*cutoff(t-k)^d; 
B2(t,d) = piecewise(d==0, B1(t,0), d>0, symsum(f2,k,0,d+1));
%B2(t,d) = symsum((piecewise(d==0, f1(t,0,k), f2(t,d,k))),k,0,d+1); %
%misses the case d=0
IB2(t,d) = int(B2(t,d), 0,t);
IB2_0(t) = vpaintegral(@(tau)B2(tau,0),tau,0,t,'Waypoints',W);

% Matt J's solution:
cutoff(t,d) = piecewise(t<0,0,t^d);
f3(t,d,k) = f(t,d,k)*cutoff(t-k,d); 
B3(t,d) = symsum(f3,k,0,d+1);
IB3(t,d) = int(B3(t,d), 0,t);

first_default_colour = [0    0.4470    0.7410];

figure(1)
hold on
fplot(B1(t,0),[0,9], 'LineWidth', 2); % first_default_colour
fplot(B1(t,1),[0,9], 'LineWidth', 2);  
fplot(B1(t,2),[0,9], 'LineWidth', 2);  
fplot(B1(t,4),[0,9], 'LineWidth', 2);  
fplot(B1(t,8),[0,9], 'LineWidth', 2);  
plot([0 0], [0 1], 'color', first_default_colour,'LineWidth', 2); % explicitly plot left vertical leg of square pulse
xlim([-1,10]);
ylim([-0.2,1.2]);
legend('B1(t,0)','B1(t,1)','B1(t,2)','B1(t,4)','B1(t,8)');

figure(2)
hold on
fplot(B2(t,0),[0,9], 'LineWidth', 2); % first_default_colour
fplot(B2(t,1),[0,9], 'LineWidth', 2);  
fplot(B2(t,2),[0,9], 'LineWidth', 2);  
fplot(B2(t,4),[0,9], 'LineWidth', 2);  
fplot(B2(t,8),[0,9], 'LineWidth', 2);  
plot([0 0], [0 1], 'color', first_default_colour,'LineWidth', 2); % explicitly plot left vertical leg of square pulse
xlim([-1,10]);
ylim([-0.2,1.2]);
legend('B2(t,0)','B2(t,1)','B2(t,2)','B2(t,4)','B2(t,8)');

figure(3)
hold on
fplot(B3(t,0),[0,9], 'LineWidth', 2); % first_default_colour
fplot(B3(t,1),[0,9], 'LineWidth', 2);  
fplot(B3(t,2),[0,9], 'LineWidth', 2);  
fplot(B3(t,4),[0,9], 'LineWidth', 2);  
fplot(B3(t,8),[0,9], 'LineWidth', 2);  
plot([0 0], [0 1], 'color', first_default_colour,'LineWidth', 2); % explicitly plot left vertical leg of square pulse
xlim([-1,10]);
ylim([-0.2,1.2]);
legend('B3(t,0)','B3(t,1)','B3(t,2)','B3(t,4)','B3(t,8)');

figure(4)
hold on
fplot(IB1_0,[0,9], 'LineWidth', 2);  
fplot(IB1_1,[0,9], 'LineWidth', 2);  
fplot(IB1_2,[0,9], 'LineWidth', 2);  
fplot(IB1_4,[0,9], 'LineWidth', 2);  
fplot(IB1_8,[0,9], 'LineWidth', 2);  
xlim([-1,10]);
ylim([-0.2,1.2]);
legend('IB1_0','IB1_1','IB1_2','IB1_4','IB1_8', 'Location','southeast');

figure(5)
hold on
fplot(IB2_0,[0,9], 'LineWidth', 2);  
fplot(IB2(t,1),[0,9], 'LineWidth', 2);  
fplot(IB2(t,2),[0,9], 'LineWidth', 2);  
fplot(IB2(t,4),[0,9], 'LineWidth', 2);  
fplot(IB2(t,8),[0,9], 'LineWidth', 2);  
xlim([-1,10]);
ylim([-0.2,1.2]);
legend('IB2_0','IB2_1','IB2_2','IB2_4','IB2_8', 'Location','southeast');

figure(6)
hold on
fplot(IB3(t,0),[0,9], 'LineWidth', 2);  
fplot(IB3(t,1),[0,9], 'LineWidth', 2);  
fplot(IB3(t,2),[0,9], 'LineWidth', 2);  
fplot(IB3(t,4),[0,9], 'LineWidth', 2);  
fplot(IB3(t,8),[0,9], 'LineWidth', 2);  
xlim([-1,10]);
ylim([-0.2,1.2]);
legend('IB3_0','IB3_1','IB3_2','IB3_4','IB3_8', 'Location','southeast');

[Matt J]

If floor(t) is used, then use of the cutoff function can be avoided.

I don't see how, for d>0.

[Dan]

See B1 and IB1 above.

[Matt J]

OK, I see. However, since you're willing to accept a numerical solution from vpaintegral, I don't know why you don't just do the whole thing numerically, which is a fair bit faster:

timeit(@()methodSym)

ans = 0.0758

timeit(@()methodNum)

ans = 0.0013

fplot(methodSym,[0,9]); hold on
[T,IB]=methodNum;
plot(T,IB,'o')
legend Sym Num; hold off

function IB=methodSym

    syms t real;
    syms k integer;
    syms d integer;
    syms tau real;  % dummy variable of integration
    parity(k) = (-1)^k;
    f(t,d,k) = 1/factorial(d)*parity(k)*nchoosek(d+1,k);
    F(t,d,k) = f(t,d,k)*(t-k)^d;
    B1(t,d) = symsum(F,k,0,floor(t));
    W = [1]; %W = [1,2];   
    
    IB(t) = vpaintegral(@(tau)B1(tau,2),tau,0,t,'Waypoints',W);

end

function [T,IB]=methodNum

    d=2;

    parity = @(k) (-1).^k; 
    f = @(t,k) parity(k).*(d+1)./factorial(d+1-k)./factorial(k);
    F= @(t,k) f(t,k).*(t-k).^d;
    B = @(t) arrayfun( @(tau) sum( F(tau,0:min(floor(tau),d+1 ) ) )    ,t);
    
    T=linspace(0,9,1e2);
    
    IB=cumtrapz(T,B(T));

end

[Dan]

why you don't just do the whole thing numerically

My ignorance.

I care more about clarity and simplicity than speed (although some of the mixed symbolic/numeric solutions I tried took much too long).

I don't know how to determine (a priori) Matlab's capability to handle mathematical expressions symbolically. I also don't know how to determine appropriate 'Waypoints'. Unless and until I become more proficient at those, I'm hesitant to trust Matlabs symbolic capabilities.

Now that I am aware of the purely numeric solution (works for both the formula with the floor and the one with the cutoff) then I'll certainly stick with that. It seems much more straightforward.