I am trying to determine the Fourier Series coefficients of flowrate pulsations from a triplex reciprocating pump. I have the flowrate defined in terms of a function that has conditional statements. The evaluation of the FourierSeries of the function hangs up. I would really appreciate some help in figuring out how to evaluate the Fourier Series.
I would use radiants instead of degrees. Also, the If constructs do not play well with Integrate. I would write
A = 1; f = 5; r = 1; L = 5 r;
\[Theta] =.; \[Theta]2 = \[Theta] +
120 Degree; \[Theta]3 = \[Theta]2 + 120 Degree;
\[Delta] = -30 Degree;
P1 = -r Sin[\[Theta]] - (r^2 Sin[\[Theta]] Cos[\[Theta]])/Sqrt[
L^2 - r^2 (Sin[\[Theta]])^2];
P2 = -r Sin[\[Theta]2] - (r^2 Sin[\[Theta]2] Cos[\[Theta]2])/Sqrt[
L^2 - r^2 (Sin[\[Theta]2])^2];
P3 = -r Sin[\[Theta]3] - (r^2 Sin[\[Theta]3] Cos[\[Theta]3])/Sqrt[
L^2 - r^2 (Sin[\[Theta]3])^2];
tot[\[Theta]_] = Max[P1, 0] + Max[P2, 0] + Max[P3, 0]
Even with this, FourierSeries seems too slow. You can calculate the Fourier coefficients directly fast enough with NIntegrate:
NIntegrate[
tot[\[Theta]] E^(-3 I \[Theta])/(2 Pi), {\[Theta], 0, 2 Pi/3,
4 Pi/3, 2 Pi}]
Thank you for your help. Much appreciated!
I have better understanding of the problem. Per your suggestion I have also moved to using radians and have defined the function in terms of Max[f,0]. I am now trying to calculate the individual Fourier coefficients. But I am still having problem in evaluating the integrals for the Fourier coefficients.
I have attached my new file. Any advice?
The command NIntegrate only works for fully numerical functions. Symbolic parameters are not accepted. You can try like this:
a[0, \[Chi]_] :=
a[0, \[Chi]] = (1/\[Chi]) NIntegrate [fsum[t], {t, 0, \[Chi]} ];
a[n_Integer, \[Chi]_] :=
a[n] = (1/\[Chi]) NIntegrate [
fsum[t] Cos[( \[Pi] n t)/\[Chi]], {t, 0, \[Chi]}];
b[n_Integer, \[Chi]_] :=
b[n] = (1/\[Chi]) NIntegrate [
fsum[t] Sin[( \[Pi] n t)/\[Chi]], {t, 0, \[Chi]}];
a[0, 1]
b[2, 3]
For purposes of integration or FourierTransform and related functions, it is almost always better to use Piecewise than If, Which, Conditional, or other programmatic and pattern constructs. Also HeavisideTheta is sometimes useful in this context.
Thanks again for the suggestions and help.
I modified my file to explicitly compute the Fourier coefficients like this:
a[n_Integer] : = (1/[Chi]) NIntegrate [
fsum[t] Cos[( [Pi] n t)/[Chi]], {t, 0, [Chi]}]
b[n_Integer] : = (1/[Chi]) NIntegrate [
fsum[t] Sin[( [Pi] n t)/[Chi]], {t, 0, [Chi]}]
Print[{a0, a[1], b[1], a[2], b[2]}]
and the computation still hangs and provides errors on convergence.
Did you provide a numeric value for [Chi]?
Yes, I did. Attached is the new file.
You can suppress the convergence messages by setting and AccuracyGoal and giving the corner points explicitly. The following computes without errors, but fsum and fourier do not match at all
a0 = (1/[Chi]) NIntegrate [
fsum[t], {t, 0, 1/6, 1/3, 1/2, 2/3, 5/6, [Chi]} ,
AccuracyGoal → 8];
a[n_Integer] :=
a[n] = (1/\[Chi]) NIntegrate [
fsum[t] Cos[( \[Pi] n t)/\[Chi]], {t, 0, 1/6, 1/3, 1/2, 2/3,
5/6, \[Chi]}, AccuracyGoal -> 8]
b[n_Integer] :=
a[n] = (1/\[Chi]) NIntegrate [
fsum[t] Sin[( \[Pi] n t)/\[Chi]], {t, 0, 1/6, 1/3, 1/2, 2/3,
5/6, \[Chi]}, AccuracyGoal -> 8]
Print[{a0, a[1], b[1], a[2], b[2]}]
fourier[nmax_] :=
a0 / 2 + Sum[
a[n] Cos[ (\[Pi] n t)/\[Chi]] + b[n] Sin[(\[Pi] n t)/\[Chi]], {n,
1, nmax}]
Plot [Evaluate[{fsum[t], fourier[2]}], {t, 0, 1}, Frame -> True,
Axes -> True, PlotRange -> { {0, 1}, {-1.25, 1.25}},
FrameLabel -> {"Time", "Flow rate"}, GridLines -> Automatic,
AspectRatio -> .65,
PlotStyle -> {{Red, Thickness[0.005]}, {Green,
Thickness[0.0051]}, {Blue, Thickness[0.005]}, {Gray,
Thickness[0.01]}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 15,
FontFamily -> "Arial"}]
I modified the final section and got it to converge for 10 Fourier series terms (nmax).
a[n_Integer] := (1/[Chi]) NIntegrate [
fsum[t] Cos[( [Pi] n t)/[Chi]], {t, 0, 1/6, 1/3, 1/2, 2/3,
5/6, [Chi]}, AccuracyGoal → 8]
b[n_Integer] := (1/[Chi]) NIntegrate [
fsum[t] Sin[( [Pi] n t)/[Chi]], {t, 0, 1/6, 1/3, 1/2, 2/3,
5/6, [Chi]}, AccuracyGoal → 8]
Clear[fourier, nmax];
fourier[nmax_Integer, t_] :=
a0 / 2 + Sum[
a[n] Cos[ ([Pi] n t)/[Chi]] + b[n] Sin[([Pi] n t)/[Chi]], {n,
1, nmax}]
nmax = 10;
This is the result I get and the match between the original function and Fourier series approximation is nowehere near close.
I tried increasing nmax to 100 and the computation hangs. Any advice?
I am familiar with Fourier series only in radians. Using radians I get a good match with this code:
r = 1;
L = 5 r;
Clear[t];
P1[t_] = -r Sin[\[Theta]] - (r^2 Sin[t] Cos[t])/Sqrt[
L^2 - r^2 Sin[t]^2];
P2[t_] = P1[t + 2 Pi/3];
P3[t_] = P1[t + 4 Pi/3];
fsum[t_] = Max[P1[t], 0] + Max[P2[t], 0] + Max[P3[t], 0];
coeff[n_Integer] :=
coeff[n] =
1/(2 Pi) NIntegrate[
fsum[s] Exp[-I n s], {s, 0, Pi/3, 2 Pi/3, Pi/2, 4 Pi/3, 5 Pi/3,
2 Pi}];
fourier[nmax_Integer] :=
fourier[nmax] =
Function[t,
Evaluate[Re@Sum[coeff[n] Exp[I n t], {n, -nmax, nmax}]]];
Animate[Plot[Evaluate[{fsum[t], fourier[n][t]}], {t, 0, 2 Pi},
PlotRange -> {0, 1.3}], {n, 1, 18, 1}]
Sorry, there was a stray theta left over when I replaced theta with t.
This should work:
Clear[t, fsum, coeff, fourier, n]
r = 1;
L = 5 r;
P1[t_] = -r Sin[t] - (r^2 Sin[t] Cos[t])/Sqrt[L^2 - r^2 Sin[t]^2];
P2[t_] = P1[t + 2 Pi/3];
P3[t_] = P1[t + 4 Pi/3];
fsum[t_] = Max[P1[t], 0] + Max[P2[t], 0] + Max[P3[t], 0];
coeff[n_Integer] :=
coeff[n] =
1/(2 Pi) NIntegrate[
fsum[s] Exp[-I n s], {s, 0, Pi/3, 2 Pi/3, Pi/2, 4 Pi/3, 5 Pi/3,
2 Pi}];
fourier[nmax_Integer] :=
fourier[nmax] =
Function[t,
Evaluate[Re@Sum[coeff[n] Exp[I n t], {n, -nmax, nmax}]]];
Animate[Plot[Evaluate[{fsum[t], fourier[n][t]}], {t, 0, 2 Pi},
PlotRange -> {0, 1.3}], {n, 1, 18, 1}]
Thanks. That code works. I really appreciate you helping out.
Now that we know that the series evaluation works in the radian domain, it should be easy to extend it to the time domain. I tried modifying the code to this, but no luck. I think I am close, but missing something simple.
[Chi] = 1; n =.; t =.; f = 1/[Chi]; [Omega] = 2 [Pi] f; r = 1.; L = 5 r; [Theta] = [Omega] t;
Clear[fourier, n, nmax];
[Theta]2 = [Theta] + 2 [Pi]/3.; [Theta]3 = [Theta]2 + 2 [Pi]/3;
Clear[P1, P2, P3, coeff]; P1[[Theta]_] := -r Sin[[Theta]] - (r^2 Sin[[Theta]] Cos[[Theta]])/ Sqrt[L^2 - r^2 (Sin[[Theta]])^2]; P2[[Theta]_] := P1[[Theta] + 2 Pi/3];
P3[[Theta]_] := P1[[Theta] + 4 Pi/3];
fsum[t] := Max[P1[[Theta]], 0] + Max[P2[[Theta]], 0] + Max[P3[[Theta]], 0];
coeff[n_Integer] := coeff[n] = (1/[Chi]) NIntegrate[ fsum[s] Exp[-I ([Pi] n s)/[Chi]], {s, 0, [Chi]/6, [Chi]/3, [Chi]/2, 2 [Chi]/3, 5 [Chi]/6, [Chi]}];
fourier[nmaxInteger , t_] := Evaluate[Re@Sum[coeff[n] Exp[I ([Pi] n t)/[Chi]], {n, 0, nmax}]];
Manipulate[ Plot[Evaluate[{fsum, fourier[n, t]}], {t, 0, [Chi]}, PlotRange -> {0, 1.3}], {n, 10, 20, 1}]
You were missing some omegas:
Clear[fsum, fourier, n, nmax, coeff];
\[Chi] = 1; n =.; t =.; f = 1/\[Chi]; \[Omega] =
2 \[Pi] f; r = 1.; L = 5 r; \[Theta] = \[Omega] t;
\[Theta]2 = \[Theta] + 2 \[Pi]/3.;
\[Theta]3 = \[Theta]2 + 2 \[Pi]/3;
Clear[P1, P2, P3, coeff];
P1[\[Theta]_] := -r Sin[\[Theta]] - (r^2 Sin[\[Theta]] Cos[\[Theta]])/
Sqrt[L^2 - r^2 (Sin[\[Theta]])^2];
P2[\[Theta]_] := P1[\[Theta] + 2 Pi/3];
P3[\[Theta]_] := P1[\[Theta] + 4 Pi/3];
fsum[\[Theta]_] :=
Max[P1[\[Theta]], 0] + Max[P2[\[Theta]], 0] + Max[P3[\[Theta]], 0];
coeff[n_Integer] :=
coeff[n] = (1/\[Chi]) Quiet@
NIntegrate[
fsum[\[Omega] s] Exp[-I \[Omega] n s], {s,
0, \[Chi]/6, \[Chi]/3, \[Chi]/2, 2 \[Chi]/3,
5 \[Chi]/6, \[Chi]}, AccuracyGoal -> 8];
fourier[nmax_Integer, t_] :=
Evaluate[Re@Sum[coeff[n] Exp[I \[Omega] n t], {n, -nmax, nmax}]];
Manipulate[
Plot[Evaluate[{fsum[\[Omega] t], fourier[n, t]}], {t, 0, \[Chi]},
PlotRange -> {0, 1.3}], {n, 10, 20, 1}]
You are fantastic! Thank you so much. Everything is working fine now.