How to get correct values of a(n) from this equation. Attached is my code but it gives wrong values. where n=1...5,

clc; close all; clear d=0.5; N=10; for n=1:N/2 sum=0; fderiv = @(f,x) (f(x+eps)-f(x))/eps; % Derivative of a function xderiv = @(x) ((x+eps)-x)/eps; % Derivative of an expression for k = 1:(N/2) x(k)=sin(2*(k+n-1)*pi/N)/(k+n-1) + sin(2*(k-n)*pi/N)/(k-n); if isnan( x(k) ) % Test for ‘NaN’ x(k) = fderiv(@sin, 2*(k-n)*pi/N)/(xderiv(k-n)); % L'Hospital's Rule end end sum=sum+x; an = (2/(N*pi))*(pi-sum); end display(an);% this is not giving the correct values????? an_manually=[-0.5736 -0.4197 -0.1352 0.2383 0.6486] % these are correct values

[Jan]

Please use the "{} Code" button to format your code. You want it to be readable, don't you?

Answer 1

[alice]

Here is what I have understood of your problem:

When you calculate a(n), you have a problem in the sum because when k=n the expression sin(2*(k-n)*pi/N)/(k-n) is undefined.

So you try to replace it by its limit when k tends toward n. To evaluate this limit you use L'Hospital's rule.

What I suggest you:

• Maybe you should not test for NaN values, but test if k==n, it would be clearer.
• As Jan told you, change the name of your sum, something like mySum for example.In the loop on k, you add x to mySum, so define x and not x(k):

...
for n=1:Nele/2
    mySum = 0;
    for k = 1:(Nele/2)
        x = sin(2*(k+n-1)*pi/Nele)/(k+n-1) + sin(2*(k-n)*pi/Nele)/(k-n);
        ...
        mySum = mySum + x;
    end
    an = ...
end
...

• Then, when k=n, you forgot the first part of your expression in your code and it is unclear what is your variable. We consider the two functions corresponding to the numerator and the denominator of the expression you want to find the limit of:

f_num = @(k) sin(2*(k-n)*pi/Nele);
f_den = @(k) k-n;

You will use the derivative when k=n, so I think you don't need to evaluate it numerically, because the expression of the derivative taken in n is really simple: it is 2*pi/Nele for the numerator and 1 for the denominator. So your x is:

if k==n
    x = sin(2*(k+n-1)*pi/Nele)/(k+n-1) + 2*pi/Nele;
end

I hope it helps

[Jan]

+1. An exhaustive explanation.

[Tariq]

still there is some thing wrong. Please note that I want to find the values of an according to equation (see posted picture) and n is ranging from 1, 2, ...5. It's need to find a1, a2, a3, a4, then a5

[alice]

What is the problem exactly?

Answer 2

[Jan]

Do not redefine "sum" as a variable, because this shadowing of builtin functions causes troubles frequently.

This is a bad choice for calculating a derivative:

fderiv = @(f,x) (f(x+eps)-f(x))/eps;

Not the eps is the smallest number with 1+eps ~= eps. But 2+eps is exatcly 2 already, such that your formula will fail if abs(x) > 1.0. The same for xderiv. Use a function instead:

function df = fderiv(f, x)
e  = sqrt(eps(x));
df = (f(x+e)-f(x)) / e;
end

Withusing the squareroot of eps(x), you reduce the cancellation error.

What do you expect as output of: ((x+eps)-x)/eps? Is this either 1 or 0?

There is no need to redefine the fdreiv and xderiv in each iteration.

We cannot guess what you think the correct values are. The posted picture does not contain derivatives, so the intention of your code is not clear.

[Tariq]

Thanks Jan for your comments. Actually I want to find the values of a(n) according to the equation (see posted picture) for n=1, 2, ...5.