Plotting Partial Sums of Fourier Series
In the definition of s
you're summing from k==0
. Since the summand has a term 1/k
this gives a divide-by-zero error when calculating the partial sums. The sum should in fact start from k==1
(the zeroth coefficient is taken care of by the constant term in front of the sum). The first few approximations then look like
s[n_, x_] :=
8/4 + 3/(9 π) Sum[(6 (-1)^k)/(k π) Cos[(k π x)/
2] + (16 (-1)^k + 13)/(π k) Sin[(k π x)/2], {k, 1, n}]
partialsums = Table[s[n, x], {n, 1, 5}];
Plot[Evaluate[partialsums], {x, -4, 4}]
To compare this with f
we plot the s[5,x]
and f
in the same plot:
Plot[{s[5, x], f[x]}, {x, -2, 2}]
which doesn't seem right to me, so I suspect you made a mistake somewhere in calculating the coefficients.
Calculating coefficients by hand
We could use FourierSeries
to calculate the partial sums, but this is very slow, and doesn't give produce the general equation for the coefficients. Therefore it's better to calculate the coefficients by hand. The $n$-th coefficient can be calculated according to
coeff[0] = 1/4 Integrate[f[x], {x, -2, 2}];
coeff[n_] = 1/4 Integrate[f[x] Exp[I Pi n x/2], {x, -2, 2}]
(1/(2 n^4 π^4))E^(-I n π) (-48 + 48 E^(I n π) - 48 I n π + 28 n^2 π^2 - 6 E^(I n π) n^2 π^2 + 2 E^(2 I n π) n^2 π^2 + 12 I n^3 π^3 - I E^(I n π) n^3 π^3 - I E^(2 I n π) n^3 π^3)
Then the partial sums are given by
series[m_, x_] := Sum[Exp[-I Pi n x/2] coeff[n], {n, -m, m}]
Plotting the first few approximations:
Plot[Evaluate[Table[series[j, x], {j, 0, 5}]], {x, -6, 6}]
To see how this compares with the original function f
:
Plot[Evaluate[{series[5, x], f[Mod[x, 4, -2]]}], {x, -4, 4}]
which looks a lot better than the before.
Edit: Real coefficients
Here, coeff[n]
are the coefficients for the Fourier series in exponential form, but these can be easily converted to the coefficients for the $\cos$ and $\sin$ series, a_n
and b_n
, by doing something like
a[0] = coeff[0];
a[n_] = Simplify[ComplexExpand[coeff[n] + coeff[-n]]];
b[n_] = Simplify[ComplexExpand[I (coeff[n] - coeff[-n])]];
You can use FourierCoefficient
to pre-calculate the Fourier coefficients to arbitrary degree and then use the result very effectively.
There may be some issues with zero-th degree, therefore I excluded this using Piecewise
. Here's the main code block:
f[x_] := Piecewise[{
{-x^3 - 2 x, -2 < x < 0},
{-1 + x, 0 <= x <= 2}},
0
];
Module[{x, fp},
(* Set parameters so that the integration runs
from -2 to 2 *)
fp = {0, -Pi/2};
fc = Piecewise[{
{FourierCoefficient[f[x], x, 0, FourierParameters -> fp], # == 0},
{ComplexExpand@FourierCoefficient[f[x], x, #, FourierParameters -> fp], True}
}] &;
fc = Evaluate /@ fc;
];
The output (fc
) of this is something very unpleasant to look at; however, there's nothing Fourier-related left there, all the hard math has been done already, and only a bunch of elementary functions remain. fc
is now a function of one argument that gives you the $n$-th Fourier coefficient of f
in no time.
(* Calculate the first 2001 Fourier coefficients *)
AbsoluteTiming[Table[fc[n], {n, -1000, 1000}]] // First
0.777059 seconds
To convert this back to the function, you have to do the partial sum with your hands, for example here are the second, fourth and eightieth partial sums:
myPartialSums = Table[
(* The Pi/2 compensates for the custom FourierParameters, see
documentation of FourierSeries/FourierParameters
under "more info" *)
Re[Sum[fc[k] Exp[-Pi/2 I k t], {k, -n, n}]],
{n, {2, 4, 80}}
];
A very large output has been generated, but we're luckily not interested in it anyway but would rather plot it
Plot[
{f[t]} ~Join~ myPartialSums,
{t, -10, 10},
PlotRange -> All,
Evaluated -> True,
PlotStyle -> {Thick, Automatic, Automatic, Automatic}
]
If the summation index goes from 1 to n, as you suggested (?) above, and you make the mod to the Plot function I mentioned
s[n_, x_] := 8/4 + 3/(9 π) Sum[(6 (-1)^k)/(k π) Cos[(k π x)/
2] + (16 (-1)^k + 13)/(π k) Sin[(k π x)/2], {k, 1, n}]
partialsums = Table[s[n, x], {n, 1, 5}];
f[x_] = Piecewise[{{-x^3 - 2 x, -2 < x < 0}, {-1 + x, 0 <= x <= 2}}]
Plot[Evaluate@Tooltip[Append[partialsums, f[x]]], {x, -4 Pi, 4 Pi}]
Then you get the following figure: