Wierd wiggle in plot of symbolic expression

Question

Dan 2025-7-14

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178444-wierd-wiggle-in-plot-of-symbolic-expression

编辑： Dan 2025-7-18

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]);

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Dan 2025-7-14

移动：Matt J 2025-7-14

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 2025-7-14

4
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178444-wierd-wiggle-in-plot-of-symbolic-expression#answer_1567941

编辑：Torsten 2025-7-14

在 MATLAB Online 中打开

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten 2025-7-14

编辑：Torsten 2025-7-14

在 MATLAB Online 中打开

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 2025-7-14

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

请先登录，再进行评论。

Answer 2

Matt J 2025-7-14

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178444-wierd-wiggle-in-plot-of-symbolic-expression#answer_1567906

编辑：Matt J 2025-7-14

在 MATLAB Online 中打开

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

●

4 个评论
显示 2更早的评论隐藏 2更早的评论

Dan 2025-7-14

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

Matt J 2025-7-14

编辑：Matt J 2025-7-15

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 2025-7-14

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178444-wierd-wiggle-in-plot-of-symbolic-expression#answer_1567948

编辑：Matt J 2025-7-14

在 MATLAB Online 中打开

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]);

●

12 个评论
显示 10更早的评论隐藏 10更早的评论

Dan 2025-7-15

编辑：Matt J 2025-7-15

在 MATLAB Online 中打开

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]);

Dan 2025-7-17

在 MATLAB Online 中打开

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 2025-7-18

编辑：Matt J 2025-7-18

在 MATLAB Online 中打开

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 2025-7-18

编辑：Dan 2025-7-18

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.

请先登录，再进行评论。

Wierd wiggle in plot of symbolic expression

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

3 个评论
显示 1更早的评论隐藏 1更早的评论

更多回答（2 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Wierd wiggle in plot of symbolic expression

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

3 个评论 显示 1更早的评论隐藏 1更早的评论

更多回答（2 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

12 个评论 显示 10更早的评论隐藏 10更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论