Numerical integration of function to determine Fourier coefficients

Hello,

Looking for an efficient and succinct way to numerically calculate N Fourier coefficients with integrate().

For example, I can easily calculate a_k for k=1 as shown below, however I would like to succinctly calculate an array of a_k for k=1:N, without using a for loop.

k =1;
fun = @(x) cos(k*x)./sin(x);
a_k = integral(fun,x1,x2);

Any suggestions?

Cheers!

Edit: Sorry, I meant integral() function in the first sentence.

Answer 1

[John D'Errico]

Of course, the limits will matter. I will guess if they include 0, for example, you may have a problem.

syms k x
x1 = pi/4;
x2 = pi/2;
F = cos(k*x)/sin(x);

Fcoef = arrayfun(@(K) int(subs(F,k,K),x1,x2),1:10)
Fcoef =
[ log(2)/2, log(2^(1/2) + 1) - 2^(1/2), log(2)/2 - 1, log(2^(1/2) + 1) - (2*2^(1/2))/3, log(2)/2, log(2^(1/2) + 1) - (7*2^(1/2))/15, log(2)/2 - 1/3, log(2^(1/2) + 1) - (64*2^(1/2))/105, log(2)/2 - 1/3, log(2^(1/2) + 1) - (227*2^(1/2))/315]

Note the odd terms show one pattern, the even terms another.

Fcoef(1:2:end)
ans =
[ log(2)/2, log(2)/2 - 1, log(2)/2, log(2)/2 - 1/3, log(2)/2 - 1/3]
Fcoef(2:2:end)
ans =
[ log(2^(1/2) + 1) - 2^(1/2), log(2^(1/2) + 1) - (2*2^(1/2))/3, log(2^(1/2) + 1) - (7*2^(1/2))/15, log(2^(1/2) + 1) - (64*2^(1/2))/105, log(2^(1/2) + 1) - (227*2^(1/2))/315]

I don't want to forget the zero'th order term.

int(subs(F,k,0),x1,x2)
ans =
log(2^(1/2) + 1)

If you prefer to do this numerically using integral, you can do so easily enough. But the integrate function would not apply in any case.

[Derry Lappin]

Thanks, this was helpful. Computation time is important as I need the option to use many coefficients. Numerically integrating seems to execute much faster than the symbolic solution.

I'm glad you showed me an example using syms, now I've got both in my toolbox.

Thanks again!

[John D'Errico]

The version for numerical integration is essentially the same, just using integral.

x1 = pi/4;
x2 = pi/2;
F = @(x,k) cos(k*x)./sin(x);
Fcoef = arrayfun(@(K) integral(@(x) F(x,K),x1,x2),0:10)

Fcoef =
  Columns 1 through 4
         0.881373587019543         0.346573590279973        -0.532839975353552        -0.653426409720027
  Columns 5 through 8
       -0.0614354545625204         0.346573590279973         0.221407257912099        0.0132402569466394
  Columns 9 through 11
        0.0193767490016565        0.0132402569466393        -0.137758091262021
        
log(2^(1/2) + 1)
ans =
         0.881373587019543
         
log(2)/2
ans =
         0.346573590279973

I've made the first term here the zero'th order term. As you can see they are the same in numerical form as the analytically computed values.