Harmonický ustálený stav

Klasické řešení pomocí NDSolve

Popis obvodu

ClearAll["Global`*"];

ω=314.;
T=2*(Pi/ω);
f0=(1/T);
ampl=10;
fi=0.3;
dos={R->15,c->(10^(-3)),L->0.1,u[t]->ampl*Sin[2*Pi*(t/T)+fi]};
rceSoucastek={uR[t]==iR[t]*R,uL[t]==L*iL'[t],ic[t]==c*uc'[t]};
rceNapeti={0==u[t]-uR[t]-uL[t]-uc[t]};
rceProudy={iR[t]==iL[t],iL[t]==ic[t]};
rcePocatecnichPodminek={uc[0]==0,iL[0]==0};

Řešení

vektorDerivaci={uc'[t],iL'[t]};
nezname=vektorDerivaci/.a_'[t]:>a[t];
prave=vektorDerivaci/.Solve[Union[rceSoucastek,rceProudy,rceNapeti],{uc'[t],iL'[t]},{iR[t],uR[t],uL[t],ic[t]}][[1]];
rceDerivaci=Thread[Equal[vektorDerivaci,prave]]

vse=Union[rceDerivaci,rcePocatecnichPodminek]/.dos;
tmax=18T;
res=NDSolve[vse,nezname,{t,0,tmax}][[1]];
napeti=uc[t]/.res;

Zobrazení výsledků

"... Nejprve se podívejme na výsledné grafy, vidíme, že zpočátku není průběh napětí na kondenzátoru sinusový, po odeznění tzv.přechodného děje se ovšem stane sinusovým, se stejnou frekvencí, jako má zdroj napětí a obecnou amplitudou a fázovým posunem. ..."

Plot[{napeti,0.1u[t]/.dos},{t,0,tmax},AxesLabel->{"t [s]","uc [V],  0.1u[t]"},AxesOrigin->{0,0},PlotRange->All,PlotStyle->{{Hue[0],Thickness[0.004]},{Hue[0.63],Thickness[0.004]}}]

"...Harmonické funkce nás také zajímají proto, že jsou „informačně úsporné “: podíváte - li se v CAO-RLC-NDSolve+HUS.nb na FullForm[res], dostanete odpověď „výstup je příliš rozsáhlý “. ..."

(*FullForm[res]*)

Dvě poznámky pod čarou

1. "... I když je nám to jasné, není špatné si pohrát s Manipulate v CAO-RLC-NDSolve+HUS.nb a osvěžit si vliv amplitudy a fáze. ..."

Manipulate[Plot[{Sin[ω*t],Um*Sin[ω*t+fi]},{t,0,T},AxesLabel->{"t","Sin[ω*t],Um*Sin[ω*t+fi]"},PlotRange->{Automatic,{-2.1,2.1}}],{Um,0.1,2},{fi,-Pi,Pi}]

2. Sečtením dvou sinusovek o stejné frekvenci f (kruhové rychlosti ω) dostaneme sinusovku o téže frekvenci f (kruhové rychlosti ω). 

Funguje to i obráceně: Mějme zadanou sinusovku kruhové rychlosti ω, náhodné amlitudě a náhodné fázi (zadano[t_] = Abs[rnd]*Sin[ω*t + rnd]). Máme-li dvě jiné sinusovky a stejné kruhové rychlosti ω a obecně náhodných amplitudách a fázích (sinusoida1 = Abs[rnd]*Sin[ω*t + rnd]), potom můžeme sinusovku zadano z nich složit: zadano = a*sinusoida1 + b*sinusoida2. Níže  počítáme koeficienty a a b, v příkazu Manipulate si je můžete nastavit jako k1 a k2:

SeedRandom[1];
rnd:=Random[Real,{-1,1}];
zadano[t_]=Abs[rnd]*Sin[ω*t+rnd];
sinusoida1=Abs[rnd]*Sin[ω*t+rnd];
sinusoida2=Abs[rnd]*Sin[ω*t+rnd];
dat=Table[{t,zadano[t]},{t,0,T,(T/20)}];
model=a*sinusoida1+b*sinusoida2;
dosFit=FindFit[dat,model,{a,b},t]
nahr=model/.dosFit
Plot[nahr-zadano[t],{t,0,T},AxesLabel->{"t[s]","chyba nahrady"}]
Manipulate[Plot[{zadano[t],k1*sinusoida1+k2*sinusoida2},{t,0,T},AxesLabel->{"t[s]",""},PerformanceGoal->"Speed",PlotRange->{Automatic,{-2,2}}],Evaluate[{k1,a-1,a+1}/.dosFit],Evaluate[{k2,b-1,b+1}/.dosFit]]

Ukázka nalezení amplitudy a fáze hledaného napětí a proudu "hrubou silou", pomocí hledání minima pro zvídavější čtenáře.

"...V notebooku CAO-RLC-NDSolve+HUS.nb je ukázka nalezení řešení „hrubou silou “ (a bez znalosti dalšího teoretického aparátu), kdy do obvodových rovnic dosadíme za hledaný proud a napětí harmonické funkce s neznámou amplitudou a fází. Každou rovnici nahradíme druhou mocninou rozdílu pravé a levé strany..."

ClearAll[rceHUS,ciselneHUS];
dosHus={iL->(ILm*Sin[ω*#+fiIL]&),uc->(Ucm*Sin[ω*#+fiUc]&),u[t]->Uz*Sin[ω*t+fiUz]};
rceHUS=rceDerivaci/.dosHus;
ciselneHUS[t_]=rceHUS//.dos/.{fiUz->0,Uz->10};
proMinimum=Plus@@(Flatten[Table[ciselneHUS[t],{t,0,T,(T/(20+Random[]))}]]/.a_==b_:>(((a-b))^2));
resHus=NMinimize[{proMinimum,ILm≥0&&Ucm≥0&&fiUc>-Pi&&fiUc≤Pi},{{ILm,0.1,20},{Ucm,0.1,20},{fiUc,0,Pi},{fiIL,0,Pi}},MaxIterations->(10^4),AccuracyGoal->20,PrecisionGoal->18,WorkingPrecision->10]
ucZHus=(uc[t]/.dosHus)/.resHus[[2]]
Plot[{napeti,ucZHus},{t,0,tmax/2},AxesLabel->{"t [s]","uc [V] ,  ucHUS[t]"},AxesOrigin->{0,0},PlotRange->All,PlotStyle->{{Hue[0],Thickness[0.004]},{Hue[0.63],Thickness[0.004]}}]

Funkce Arg

"... Některé vlastnosti funkce Arg jsou ukázány v CAO-RLC-NDSolve+HUS.nb. ..."

Print["Arg[-∞-I]= ",Arg[-∞-I],"     Arg[-∞+I]= ",Arg[-∞+I]]
Print["Arg[-1.03*(10^6)-I]= ",Arg[-1.03*(10^6)-I],"     Arg[1.03*(10^6)+I]= ",Arg[1.03*(10^6)+I]]
Print["Arg[0]= ",Arg[0],"     Arg[163]= ",Arg[163]]
Print["Arg[0.*I]= ",Arg[0.I],"     Arg[161*I]= ",Arg[161*I]]
Print["Arg[-158.*I]= ",Arg[-158.*I],"     ((Arg[0.  + 161. *I])/π)= ",((Arg[0.+161.*I])/π)]

(*Plot3D[Arg[x+I y],{x,-3,3},{y,-3,3},AxesLabel->{"Re[z]","Im[z]","Arg[z]"},PerformanceGoal->"Quality"]
Plot3D[Abs[x+I y],{x,-3,3},{y,-3,3},AxesLabel->{"Re[z]","Im[z]","Abs[z]"},PerformanceGoal->"Quality"]*)

Převod popisu do HUS

rceHUS=(((Union[rceSoucastek,rceProudy,rceNapeti]/.a_'[t]:>I*ω*a)/.a_[t]:>a)//.dos)/.u->ampl*Exp[I*fi]
pocetRovnic=Length[rceHUS];
Print["rovnic mame:  ",pocetRovnic];
nezname=Cases[rceHUS,_Symbol,{0,∞}]//Union;
Print["nezname jsou:  "Sequence@@nezname]
Print["neznamych je:  ",Length[nezname]];
resHUS=Solve[rceHUS]

ucHUS[t_]:=Im[(uc/.resHUS[[1]])*(E^(I*ω*t))];
Plot[{napeti,ucHUS[t]},{t,0,tmax/2},AxesLabel->{"t [s]","uc [V] ,  ucHUS[t]"},AxesOrigin->{0,0},PlotRange->All,PlotStyle->{{Hue[0],Thickness[0.004]},{Hue[0.63],Thickness[0.004]}}]