Message Boards Message Boards

Fourier Series of function defined with conditional statement

Kaushik Mallick
Posted 10 years ago

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. enter image description here

Attachments:
POSTED BY: Kaushik Mallick
15 Replies
Sort By:
Kaushik Mallick
Kaushik Mallick
Posted 10 years ago

You are fantastic! Thank you so much. Everything is working fine now.
POSTED BY: Kaushik Mallick

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}]
POSTED BY: Gianluca Gorni
Kaushik Mallick
Kaushik Mallick
Posted 10 years ago

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}]
POSTED BY: Kaushik Mallick

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}]
POSTED BY: Gianluca Gorni
Kaushik Mallick
Kaushik Mallick
Posted 10 years ago

I cut and paste your code in a new notebook and the computation hangs. What can I be doing wrong?

enter image description here
POSTED BY: Kaushik Mallick

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}]
POSTED BY: Gianluca Gorni
Kaushik Mallick
Kaushik Mallick
Posted 10 years ago

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[nmaxInteger, 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.

Fourier series approximation

I tried increasing nmax to 100 and the computation hangs. Any advice?
POSTED BY: Kaushik Mallick

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"}]
POSTED BY: Gianluca Gorni
Kaushik Mallick
Kaushik Mallick
Posted 10 years ago

Yes, I did. Attached is the new file.

Attachments:
POSTED BY: Kaushik Mallick

Did you provide a numeric value for [Chi]?
POSTED BY: Gianluca Gorni
Kaushik Mallick
Kaushik Mallick
Posted 10 years ago

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.
POSTED BY: Kaushik Mallick

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.
POSTED BY: Daniel Lichtblau

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]
POSTED BY: Gianluca Gorni
Kaushik Mallick
Kaushik Mallick
Posted 10 years ago

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?

Attachments:
POSTED BY: Kaushik Mallick

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}]
POSTED BY: Gianluca Gorni
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract