Function decomposition to Fourier series using first impulse function

Question

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:

So, I expect something like this:

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.

How do I fix it? Thanks in advance.

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

Answer 1

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:

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

Answer 2

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 \
work*)
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, 
 10}}];
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->0)