how to get the correct plot of function sin(N+1)*p​i*x^2/sin(​pi*x)^2 with different N

3 次查看（过去 30 天）
Daniel Niu
Daniel Niu 2023-9-10
评论： Daniel Niu 2023-9-10
Ran in:
% Define the value of N
format long
N = 1000;
% Define a range of x values
%x = linspace(-0.01*pi,0.01*pi, 1000); % You can adjust the range and number of points as needed
x = -1.5:0.00001:1.5; % You can adjust the range and number of points as needed
% Calculate the function values
y = (sin((N+1) *pi* x).^2) ./ (sin(pi*x).^2);
% Plot the function
plot(x, y)
title('Plot of (sin((N+1) *pi* x).^2) ./ (sin(pi*x).^2)')
grid on
when I change N=10000;
% Define the value of N
format long
N = 10000;
% Define a range of x values
%x = linspace(-0.01*pi,0.01*pi, 1000); % You can adjust the range and number of points as needed
x = -1.5:0.00001:1.5; % You can adjust the range and number of points as needed
% Calculate the function values
y = (sin((N+1) *pi* x).^2) ./ (sin(pi*x).^2);
% Plot the function
plot(x, y)
title('Plot of (sin((N+1) *pi* x).^2) ./ (sin(pi*x).^2)')
grid on
the limit of (sin((N+1) *pi* x).^2) ./ (sin(pi*x).^2) is (N+1)^2, which means for N=1000, the maximum is 1000000, which is consistent with the plot
but for N=10000, the maximum is 1e8, which has something wrong.
  4 个评论
显示 2更早的评论
Torsten
Torsten 2023-9-10
编辑：Torsten 2023-9-10
The precision of usual floating point arithmetics reaches its limits. If you want a more precise plot, use symbolic computation.
And don't evaluate your y exactly at integer points. It will give a division by zero for which the result you get is unpredictable.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

Steven Lord
Steven Lord 2023-9-10
Ran in:
To eliminate the effect of the floating point approximation to π, try using sinpi instead of sin.
% Define the value of N
format long
N = 1000;
% Define a range of x values
%x = linspace(-0.01*pi,0.01*pi, 1000); % You can adjust the range and number of points as needed
x = -1.5:0.00001:1.5; % You can adjust the range and number of points as needed
% Calculate the function values
f = @(x) (sinpi((N+1)*x).^2) ./ (sinpi(x).^2);
y = f(x);
% Plot the function
plot(x, y)
title('Plot of (sin((N+1) *pi* x).^2) ./ (sin(pi*x).^2)')
grid on
Why the spikes at the integer values? Well, those peaks aren't exactly at the integer values but very close to them, since the function is undefined (0/0) at exact integer values of x.
f([-1 0 1])
ans = 1×3
NaN NaN NaN
Sam Chak
Sam Chak 2023-9-10
Ran in:
Hi @Daniel Niu
I checked using the limit() function.
syms x N
y = limit((sin((N + 1)*pi*x)^2)/(sin(pi*x)^2), x, 0)
y = 
N = sym(9999);
y = limit((sin((N + 1)*pi*x)^2)/(sin(pi*x)^2), x, 0)
y = 
100000000
Paul
Paul 2023-9-10
编辑：Paul 2023-9-10
Ran in:
Hi Daniel,
I think you are correct that the limit of y as x approaches any integer is (N+1)^2. However, all the compuations are numeric, and numeric computations don't know about limits, they just compute what's given with the finite precisicion of the data types.
Here's the first case:
format long
N = 1000;
% Define a range of x values
x = -1.5:0.00001:1.5; % You can adjust the range and number of points as needed
% Calculate the function values
y = (sin((N+1) *pi* x).^2) ./ (sin(pi*x).^2);
There is a value of x that is exactly equal to 0. For x == 0, the value of y is
y(x==0)
ans =
NaN
and if you zoom in on the plot around x = 0 you'll see a gap.
Plot what's happening around x = 1, and add a point for the expected value.
plot(x,y,'-o',1,(N+1)^2,'x'),grid,xlim([0.9999 1.0001])
Clearly, we see something fishy at x = 1.
There is a value in x that is exactly equal to 1
i1 = find(x==1)
i1 =
250001
and its corresponding value of y
y(i1)
ans =
5.237999158081999e+03
The problem is that, in double precision, the Matlab function pi cannot exactly represent mathematical pi. Hence, the numerator of y, which should ideally be zero, isn't zero. Same for the denominator
[(sin((N+1) *pi* x(i1)).^2) , (sin(pi*x(i1)).^2)]
ans = 1×2
1.0e-28 * 0.785574047890805 0.000149975978266
Those values get divided and we get what we get, which is basically junk. If you zoom in further, you'll see that the plot at x = 1 is y(i1).
Here's the same thing for the second case, and you'll see that
N = 10000;
% Define a range of x values
%x = linspace(-0.01*pi,0.01*pi, 1000); % You can adjust the range and number of points as needed
x = -1.5:0.00001:1.5; % You can adjust the range and number of points as needed
% Calculate the function values
y = (sin((N+1) *pi* x).^2) ./ (sin(pi*x).^2);
i1 = find(x == 1)
i1 =
250001
y(i1)
ans =
2.599231705387884e+08
figure
plot(x,y,'-o',1,(N+1)^2,'x'),grid,xlim([0.9999 1.0001])
Basically, we're suffering from finite precision and not getting exact values of 0 in the numerator and denominator of y for nonzero, integer values of x.

类别

在 Help Center 和 File Exchange 中查找有关 Line Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by