Převod do Fourierovy řady + HUS

V tomto notebooku demonstrujeme dva přístupy řešení jednoduchého obvodu:
1) Klasické řešení pomocí diferenciálních rovnic (NDSolve)
2) Řešení pomocí rozkladu na součet harmonických funkcí s následnou aplikací HUS.
Poznamenejme, že pro to, abychom mohli použít rozklad na součet harmonických funkcí, musí mít vstupní napětí periodický průběh.

ClearAll["Global`*"];
$HistoryLength=2;
Off[NMaximize::"cvmit",NIntegrate::"slwcon",NIntegrate::"ncvb"]

Zadání obvodu

Uvažujme RC článek, což je známý jednoduchý obvod, který jsme už řešili.

Rezistor R nechť má hodnotu 30 kΩ, kondezátor C kapacitu 10 μF. Obvod nechť je buzen obdélníkovým napětím u(t) s periodou 2s (frekvence 0,5 Hz).

soucastky={R->30*(10^3),c->10*(10^(-6))};
perioda=2.;
frekv=(1/perioda);
u[t_]:=Tanh[100Sin[2Pi*t/perioda]]+1;
pl1=Plot[u[t],{t,-perioda,perioda},PlotStyle->{Thick,Hue[0.63]}];
nb=500;
pl2=ListPlot[{#,u[#]}&/@((perioda/nb)*(Range[2nb]-nb)),PlotStyle->{Hue[0.63],PointSize[0.007]}];

Vstupní napětí u(t) má tento průběh :

Show[pl1,AxesLabel->{"t","u[t]"},PlotRange->{Automatic,{-0.4,2.4}}]

1. Klasické řešení - NDSolve

Protože obvodem teče jeden proud, můžeme použít zjednodušený popis, kdy vynecháme rovnice pro uzly. Proud obvodem označíme i(t). Rovněž využijeme faktu, že pro napětí uzlu 1 platí: u1(t)==u(t).

Rovnice součástek :

rezistorR=u[t]-u2[t]==R*i[t];
kondC=i[t]==c*u2'[t];

Počáteční podmínky :

poc1=u2[-perioda]==0;

Řešení :

rovnice={rezistorR,kondC,poc1};
rceDosad=rovnice/.soucastky;
nezname={i[t],u2[t]};
reseni=NDSolve[rceDosad,nezname,{t,-perioda,perioda}];
nDSolvePlot=Plot[u2[t]/.reseni,{t,-perioda,perioda},PlotStyle->{Thick,Red}];

Zobrazme výsledky. Vstupní napětí je modře, výstup je červeně.

Show[pl1,nDSolvePlot, AxesLabel->{"t", "u(t), u2(t)"}, PlotRange->{Automatic,{-0.4,2.4}}]

2. Fourierova řada + HUS

Vstupní napětí u(t) rozložíme na součet harmonických funkcí. Potom na obvod aplikujeme každou harmonickou zvlášť a spočítáme pro ni průběh na výstupu. Jednotlivé výstupní průběhy poté poskládáme (posčítáme), čímž získáme výstup pro u(t).

Rozklad u(t) do Fourierovy řady

Budeme počítat koeficienty pro prvních 50 harmonických funkcí.

nHarm=50;
rn=Range[nHarm];

Zajímá-li vás, co obsahují seznamy (angl. Lists) siny a cosiny, stačí odstranit středník na konci řádku.

sinus[n_]:=Sin[((2Pi)/perioda)*t*n];
siny=Map[sinus,rn];
cosinus[n_]:=Cos[((2Pi)/perioda)*t*n];
cosiny=Map[cosinus,rn];

Spočítáme koeficienty an a bn a uložíme je do seznamů aenka a beenka.

a0=(1/perioda)*NIntegrate[u[t],{t,0,perioda}];
a[n_]:=(2/perioda)*NIntegrate[u[t]*cosinus[n],{t,0,perioda}]; (* Definice funkce pro výpočet a[n]. Zatím se výpočet neprovádí. *)
aenka=Map[a,rn];(* Teprva teď se provádí výpočet. 50x se zavolá a[n] a výsledky se uloží do seznamu aenka. *)
b[n_]:=(2/perioda)*NIntegrate[u[t]*sinus[n],{t,0,perioda}];
beenka=Map[b,rn];

Složíme harmonické funkce, abychom vytvořili aproximaci u(t). 

Použití vzorce
nahradaVstup[n_]:=a0+(∑_(i=1)^n)(a[i]*cosinus[i]+b[i]*sinus[i]);
by ale bylo výpočetně náročné - při každém zavolání funkce nahradaVstup by se znova a znova počítaly všechny koeficienty a[i] a b[i], tedy znova a znova bychom integrovali. Proto jsme všechny koeficienty (jednou) předpočítali a uložili do seznamů aenka a beenka a teď je použijeme:

nahradaVstup[n_]:=a0+Apply[Plus,Take[aenka*cosiny+beenka*siny,n]];

Tím je hotový rozklad u(t) na součet harmonických funkcí.

Zobrazení výsledku převodu u(t) do Fourierovy řady

Tato část není nutná. Uvádíme ji pouze pro ty, koho zajímá, jak vypadá aproximace u(t) pomocí součtu harmonických funkcí. Můžete přeskočit rovnou na HUS.

Příprava na zobrazování.

harmonicke=aenka*cosiny+beenka*siny;
pln[n_]:=Plot[nahradaVstup[n],{t,-perioda,perioda},PlotStyle->{Thick,GrayLevel[0]},AxesLabel->{"t","u[t]"}];
pln1[n_]:=Plot[nahradaVstup[n-1],{t,-perioda,perioda},PlotStyle->{Thick,GrayLevel[0.5]}];
plHarm[n_]:=Plot[harmonicke[[n]],{t,-perioda,perioda},PlotStyle->{Thick,Green}];
vys[n_]:=Show[pl1,pln[n],  AxesLabel->{"t", "u(t)"}, PlotRange->{Automatic,{-0.4,2.4}}];
vys2[n_]:=Show[pl1,plHarm[n],pln1[n],pln[n],  AxesLabel->{"t", "u(t)"}, PlotRange->{Automatic,{-1.4,2.4}}];

Modře je napětí u(t), černě je jeho aproximace pomocí součtu n harmonických funkcí. Nultá harmonická je stejnosměrná složka, tedy koeficient a0.

Manipulate[vys[n],{n,0,nHarm,1}]

V následujícím grafu je pro názornost navíc šedivě zobrazena aproximace pomocí n-1 harmonických funkcí a zeleně je zobrazna n-tá harmonická funkce. 
   
Přičtením zelené čáry k šedivé čáře obdržíme černou čáru. Povšimněme si, že obdélníkový průběh se skládá ze součtu pouze lichých harmonických funkcí (sudé harmonické jsou nulové).

Manipulate[vys2[n],{n,0,nHarm,1}]

HUS

Fázory harmonických funkcí, ze kterých je skládáno vstupní napětí u(t), uložíme do seznamu (listu) fazoryVstup.

fazoryVstup=beenka+I*aenka;

Kruhová frekvence základní harmonické je ω. 

Kruhová frekvence ωn n-té harmonické funkce je n-násobkem ω.

omega1=2*Pi*frekv;
omega[n_]:=n*omega1;

Obvod budeme řešit jako impedanční dělič pro jednotlivé frekvence.

Koeficient a0 představuje stejnosměrnou složku ve vstupním signálu. Musíme proto stanovit přenos pro obvod ve stejnosměrném stavu. Tedy, rozpojíme kondenzátory a zkratujeme cívky a obdržíme:

Přes rezistor neteče žádný proud, není tedy na něm žádný úbytek napětí, a proto u20==a0. Stejnosměrný přenos je proto jednotkový.

prenos0=1;

Impedance kondenzátoru je závislá na frekvenci. Pro ωn je rovna

Zc[n_]:=(1/(I*omega[n]*c));

Přenos obvodu pro frekvenci ωn je tedy roven:

prenos[n_]:=((Zc[n])/(R+Zc[n]))/.soucastky;

Přenosy pro jednotlivé frekvence spočítáme a uložíme do seznamu prenosy:

prenosy=Map[prenos,rn];

Pro každou frekvenci spočítáme fázor výstupního napětí jako součin fázoru vstupního napětí a přenosu na dané frekvenci. 

Součiny můžeme provést hromadně vynásobením seznamů fazoryVstup a prenosy; dostaneme seznam součinů (vždy i-tá položka prvního seznamu se vynásobí s i-tou položkou druhého seznamu):

fazoryVystup=fazoryVstup*prenosy;

Harmonickou funkci o kruhové frekvenci ωn obdržíme tak, že její fázor vynásobíme výrazem (E^(I*omega[n]*t)) a vybereme imaginární část. Jednotlivé harmonické funkce, z nichž se skládá výstupní napětí u2(t), uložíme do seznamu harmVystupy.

e[n_]:=(E^(I*omega[n]*t));
eenka=Map[e,rn];
harmVystupy=Im[fazoryVystup*eenka];

Výstupní napětí u2(t) složíme z n harmonických funkcí.

nahradaVystup[n_]:=a0*prenos0+Apply[Plus,Take[harmVystupy,n]];

Zobrazení u2(t) složené z harmonických funkcí

Zbývá zobrazit výsledek. Proveďme přípravu k zobrazení.

fourierPlot[n_]:=Plot[nahradaVystup[n],{t,-perioda,perioda},PlotStyle->{Thick,GrayLevel[0]},AxesLabel->{"t","u[t]"}];
fourierPlotn1[n_]:=Plot[nahradaVystup[n-1],{t,-perioda,perioda},PlotStyle->{Thick,GrayLevel[0.5]},AxesLabel->{"t","u[t]"}];plHarmVystup[n_]:=Plot[harmVystupy[[n]],{t,-perioda,perioda},PlotStyle->{Thick,Green}];
vys3[n_]:=Show[nDSolvePlot,fourierPlot[n],PlotRange->{Automatic,{-1,2}}];
vys4[n_]:=Show[nDSolvePlot,fourierPlotn1[n],plHarmVystup[n],fourierPlot[n],PlotRange->{Automatic,{-1,2}}];

Červeně je napětí u2(t), jak ho spočítala NDSolve. Černě je aproximace u2(t) pomocí součtu n harmonických funkcí. Nultá harmonická je stejnosměrná složka, tedy součin a0*prenos0.

Manipulate[vys3[n],{n,0,nHarm,1}]

V následujícím grafu je pro názornost navíc šedivě zobrazena aproximace pomocí n-1 harmonických funkcí a zeleně je zobrazna n-tá harmonická funkce. 
   
Přičtením zelené čáry k šedivé čáře obdržíme černou čáru.

Manipulate[vys4[n],{n,0,nHarm,1}]