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

% 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.

[Torsten]

Have you learned on how to use L'Hospital to compute the limit of (sin((N+1) *pi* x).^2) ./ (sin(pi*x).^2) as x approaches an integer (...,-2,-1,0,1,2,...) ? I get (N+1)^2 ~ (1e4)^2 = 1e8. For N = 1000, the maximum should be approximately 1e6.

[Daniel Niu]

yes,Torsen. but why the maximum is different at 0 and 1 when N=10000 in the plot thank you

[Torsten]

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.

[Daniel Niu]

how do you perform that and get the right plot. thank you

Answer 1

[Steven Lord]

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

[Bruno Luong]

For older version one can use sind

f = @(x) (sind(180*(N+1)*x).^2) ./ (sind(180*x).^2)

Answer 2

[Sam Chak]

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

[Daniel Niu]

thank you,Sam,you are right I can not understand why the maximum is different at 0 and 1 when N=10000

Answer 3

[Paul]

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.