I have periodic function with certain first impulse function and period value. My task is to decompose it to approximate Fourier series (k_max = 10) using Laplace transform of first impulse function: Method

So, I expect something like this: Result expected

Here is my code:

θ[x_] := Boole[x >= 0]
{R, L, C1, Rl, T, Um} := {500, 0.2, 2*10^(-10), 10^3, 10^(-3), 10};
{t1, t2, t3} := {T/2, T/3, T/6};
u1[t_] := Um/t3*(t*θ[t] - (t - t3) θ[t - t3] - (t - t2) θ[t - t2] + (t - t1) θ[t - t1])

U1[s_] := LaplaceTransform[u1[t], t, s]
Am10[ω_] = Abs[U1[I*ω]];
Ph10[ω_] = Arg[U1[I*ω]];
Am1[ω_] := Piecewise[{{Limit[Am10[ω0], ω0 -> 0], ω == 0}, {Am10[ω], ω != 0}}] (*Removing [0/0]-uncertainity*)
Ph1[ω_] := Piecewise[{{Limit[Ph10[ω0], ω0 -> 0], ω == 0}, {Ph10[ω], ω != 0}}] (*Removing [0/0]-uncertainity*)

{ω1, kmax} := {2*Pi/T, 10};
DAm1 = Table[0, 2*kmax + 1]; (*Container for Am1[k*ω1]; I need it for next part of work*)
DPh1 = Table[0, 2*kmax + 1]; (*Container for Ph1[k*ω1]*)
For[k = -10, k <= 10, k++, DAm1[[k + kmax + 1]] = Am1[k*ω1]];
For[k = -10, k <= 10, k++, DPh1[[k + kmax + 1]] = Ph1[k*ω1]];
elem[k_, t_] = Am1[k*ω1]*E^(I*(Ph1[k*ω1] + k*ω1*t));
fu[t_] = Sum[elem[k, t], {k, -kmax, kmax}];
Plot[fu[t], {t, -5 T, 5 T}, PlotRange -> {{-T, T}, {0, 12}},ImageSize -> 500]

However, I get a lot of errors. Here is the text for copying:

Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating Log[...].

Indeterminate expression ArcTan[0,0] encountered.

Further output of N::meprec will be suppressed during this calculation.

Further output of ArcTan::indet will be suppressed during this calculation.

Real result

How do I fix it? Thanks in advance.

Note: the "*" symbol marks functions related first impulse (0 <= t <= T).

  • $\begingroup$ If I plot "ui[t]", I get a different function. Is something wrong with the definition? Note that there is a built in function "FourierSeries". $\endgroup$ May 16 at 16:45
  • $\begingroup$ @DanielHuber Thanks for answer. First of all, there is no ui[t] function. Did you mean u1[t]? Then, u1[t] has correct plot (like the green one on the post, using Plot[u1[t], {t, 0, T}]). About FourierSeries: I know this one, but I must realise decomposition manually. $\endgroup$ May 16 at 17:52
  • $\begingroup$ If I plot u1[t] from 0 to 100 I first get a positive blob and then a larger negative blob $\endgroup$ May 16 at 19:01
  • $\begingroup$ @DanielHuber Please check it link $\endgroup$ May 16 at 21:13

The trouble/errors are coming in the line where you compute:

    For[k=-10, k <= 10, k++, DPh1[[k+kmax+1]] = Ph1[k*ω1]];

We can shed some light on the error code by doing a stacktrace. Next to the name of the function that generated the error (in red) you can see an ellipsis (three dots). Click it, and choose 'Show Stack Trace', (should be the first option in version 12.2). Internally, there you are evaluating expressions of the form:

    Arctan[(Im[...] + Re[...] + Re[...] + Re[...])/(Re[...] + Im[...] + Im[...] + Im[...])]

Where in the expression above I have used ... to keep the size readable. For certain values that the table is trying to compute, the imaginary and real parts of the expressions that I put as ... are cancelling each other out and simplifying to ArcTan[0,0] (the same as) ArcTan[0/0], which is undefined.

To fix it, you can numericize your input when you are defining the functions:

    Am1[ω_] := Piecewise[{{Limit[Am10[ω0], ω0 -> N[0], ω == 0}, {Am10[N[ω]], ω != 0}}];
    Ph1[ω_] := Piecewise[{{Limit[Ph10[ω0], ω0 -> N[0], ω == 0}, {Ph10[N[ω]], ω != 0}}];

This squelches the error messages on my rig, but you still won't be able to plot the functions elem[k,t] or fu[t] as they have been defined since they will have both real and imaginary components (as given) unless you seperate them out first. But if you plot it on the unit square of the complex plane with the function ComplexPlot you can see the ripples:

Fourier Ripples.

I can't promise that solves all issues, but hopefully that gets you going in the right direction...

  • $\begingroup$ Thanks for your answer! Excuse me, but elem[k_, t_] = Am1[k*ω1]*E^(I*(Ph1[k*ω1] + k*ω1*t)) function is mathematically proven to have only real values since its imaginary component is always zero... Does that problem with real plot occure due to the way machine perceives expressions like a + 0*I? $\endgroup$ May 17 at 9:09
  • $\begingroup$ N[] usage solved my problem. Furthermore, I have turned to be right about real plot: I've got it without any additional manipulations cause all Im-components are reduced to zero! $\endgroup$ May 17 at 13:46

With this code I got some output in the plot:

\[Theta][x_] := Boole[x >= 0]
{R, L, C1, Rl, T, Um} := {500, 0.2, 2*10^(-10), 10^3, 10^(-3), 10};
{t1, t2, t3} := {T/2, T/3, T/6};
u1[t_] := 
Um/t3*(t*\[Theta][t] - (t - t3) \[Theta][t - t3] - (t - t2) \[Theta][
  t - t2] + (t - t1) \[Theta][t - t1])

U1[s_] := LaplaceTransform[u1[t], t, s]
Am10[\[Omega]_] = Abs[U1[I*\[Omega]]];
Ph10[\[Omega]_] = Arg[U1[I*\[Omega]]];
Am1[\[Omega]_] := 
 Piecewise[{{Limit[Am10[\[Omega]0], \[Omega]0 -> N[0]], \[Omega] == 
 0}, {Am10[\[Omega]], \[Omega] != 
 0}}] (*Removing[0/0]-uncertainity*)

Ph1[\[Omega]_] := 
Piecewise[{{Limit[Ph10[\[Omega]0], \[Omega]0 -> N[0]], \[Omega] == 
0}, {Ph10[\[Omega]], \[Omega] != 
0}}] (*Removing[0/0]-uncertainity*)

{\[Omega]1, kmax} := {2*Pi/T, 10};
DAm1 = Table[0, 
2*kmax + 
1];(*Container for Am1[k*\[Omega]1]I need it for next part of \
DPh1 = 
Table[0, 2*kmax + 1];(*Container for Ph1[k*\[Omega]1]*)
For[k = -10, 
k <= 10, k++, DAm1[[k + kmax + 1]] = Am1[k*\[Omega]1]];
For[k = -10, k <= 10, k++, DPh1[[k + kmax + 1]] = Ph1[k*\[Omega]1]];
elem[k_, t_] := 
Am1[k*\[Omega]1]*E^(I*(Ph1[k*\[Omega]1] + k*\[Omega]1*t));
fu[t_] = Sum[
elem[k, t], {k, {-10, -8, -7, -5, -4, -2, -1, 0, 1, 2, 4, 5, 7, 8, 
Show[ListPlot[Table[{t, fu[t]}, {t, 0, T, 10^-6}]], 
Plot[u1[t]/1000, {t, 0, T}, PlotStyle -> {Orange}]]

I have no idea whether this is correct, probably needs much tweaking. I left out multiples of 3 in the sum over k (result indeterminate, limit->0enter image description here)

  • $\begingroup$ I have executed your code, but I still get heaps of errors like those I described in my question. Furthermore, when I attempt to print DAm1 and DPh1 I get list of some strange symbolic expressions instead of list of decimal numbers. $\endgroup$ May 17 at 10:29
  • $\begingroup$ @CppNosavvier Yes, but I won't do further debugging..., trace the errors one by one $\endgroup$
    – Andreas
    May 17 at 10:42
  • $\begingroup$ I started to study WM recently so I don't understand stack traces. I don't think that to post answer with uncorrected code, appealing to uninformative stack trace is effective problem solution way. $\endgroup$ May 17 at 12:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.