Prikazali bomo primer uporabe knjižnice Sympy za analitičen izračun diferencialne enačbe prvega reda za nekaj primerov RLC vezij.Pri tem bomo v vsakem primeru pokazali malo drugačen način reševanja ali izrisa rešitve.
Za izvajanje tega zvezka ne potrebujete posebnega znanja programiranja v Pythonu, lahko pa poljubno spreminjate kodo in se sproti učite tudi uporabe programskega jezika.
Rešujemo diferencialno enačbo $R\frac{\text{d}i}{\text{d}t}+\frac{i}{C}=0$.
Vnesemo potrebne knjižnice, pripravimo spremenljivke, tvorimo enačbo in jo rešimo.
import numpy as np
from sympy import *
from IPython.display import *
import matplotlib.pyplot as plt
#matplotlib inline
from IPython.display import Math, HTML
def load_mathjax_in_cell_output(): # Te vrstice so potrebne, da se izpišejo enačbe. Potrebno le za Colab
display(HTML("<script src='https://www.gstatic.com/external_hosted/"
"mathjax/latest/MathJax.js?config=default'></script>"))
get_ipython().events.register('pre_run_cell', load_mathjax_in_cell_output)
import sympy
sympy.init_printing()
init_printing(use_latex=True)
var('R C t C1 a Ug') # priprava spremenljivk
i = Function("i")(t) # funkcija toka
de = Eq(i/C+R*i.diff(t)) # zapis diferencialne enačbe
display(de) # izpis enačbe
des = dsolve(de,i) # rešitev enačbe
display(des) # izpis rešitve
Rešitev diferencialne enačbe je ustrezna, le iz začetnega pogoja je potrebno še določiti konstanto $C_1$. Ta mora izhajati iz pogoja $u_C(t=0^+)=u_C(t=0^-)$. Torej mora biti napetost na kondezatorju ob preklopu enaka nič. Ker mora biti hkrati zadoščeno 2. Kirchoffovemu zakonu, mora biti napetost generatorja ob preklopu enaka napetosti na uporu $u_g(t=0)=U_g=u_R(t=0)=R i(t=0)$ in torej $i(t=0)=U_g/R$.
V tem primeru bomo to upoštevali "ročno", tako, da bomo sami preuredili rešitev enačbe, saj lahko člen $e^{C_1} $ zapišemo z novo konstanto $a$, torej lahko tvorimo enačbo $i(t)=a e^{\frac{-t}{RC}}$. V ta namen bomo uporabili funkcijo $Lambda$.
Ker je $i(t=0)= a$, je $a=U_g/R$. S tem smo prišli do končne oblike rešitve.
des = des.subs(C1,0) # za C1 vstavimo 0 in nato s funkcijo Lambda spremenimo rešitev v bolj praktično obliko
display(des)
f = Lambda((t,R,C,a),(a* des.rhs)) # ki je a*(desna stran enačbe - rhs=right hand side)
#a=Ug/R # a je očitno tok ob času t=0, ki je enak Ug/R
display(f(t,R,C,Ug/R))
# display(Latex('$i(t) = ' + str(latex(f(t,R,C,Ug/R))) + '$')) # izpišemo rešitev z uporabo Latex sintakse, ne dela v Colabu
#Pripravimo za iziris
x = np.linspace(0,15,100) # niz vrednosti za x-os - čas
R=1e5
C=5e-5
Ug=2
plt.grid(True)
plt.xlabel('Čas /s',fontsize=12)
plt.ylabel('Tok / A',fontsize=16)
plt.plot(x,[f(t,R,C,Ug/R) for t in x],color='#008000') # ob izrisu vstavimo vrednosti za R, C in a=Ug/R (1,4,1)
plt.ylim(0,2e-5)
plt.show()
Še mali programerski izziv: Spremeni kodo tako, da bo max vrednost abscisne osi vedno enaka 5RC (5* časovna konstanta).
Rešujemo enačbo ${{U}_{g}}={{u}_{R}}(t)+{{u}_{L}}(t)=iR+L\frac{\text{d}i}{\text{d}t}$. V tem primeru bomo preuredili rešitev tako, da bo možno neposredno uporabiti začetni pogoj, ki je $i(t=0)=0$.
Najprej zapišemo enačbo v obliki $ L \frac{\text{d}i}{\text{d}t} ={{U}_{g}}-i R$ in nato poiščemo njeno rešitev.
i = Function("i") # pripravimo funkcijo in spremenljivke ter zapišemo dif enačbo
L, R, Ug, t, i0, C1 = symbols('L, R, U_g, t, i0 C1')
eq = Eq(L*i(t).diff(t,1), Ug-R*i(t))
display(eq)
dsol = dsolve(eq, i(t)) # rešimo diff enačbo
display(dsol)
Sedaj se želimo "znebiti" konstante $C_1$ in izraziti enačbo z začetnim pogojem, t.j. $i(t=0^+)=0$.
To naredimo tako, da izrazimo enačbo za čas $t=0$ in jo rešimo za $C_1$.
it0=dsol.subs({'t':0})
display(it0)
init_solve=solve(it0, C1) # rešimo enačbo za C1
display(init_solve) # prikažemo rešitev
final = dsol.subs(C1, init_solve[0]) # v rešitev vstavimo C1 in dobimo končno enačbo
final
expand(final) # drugačen izpis
Vstavimo začetni pogoj $i(t=0^+)=i_0=0$ in dobimo končno obliko:
zacetni_pogoji = {i(0): 0}
i2=final.subs(zacetni_pogoji)
simplify(i2)
Poiščemo še rešitev za napetost z odvajanjem toka (desne strani enačbe i2) in množenjem z induktivnostjo L.
u2=symbols('u2')
u2=L*i2.rhs.diff(t) # množenje z L in odvodom desne strani enačbe
expand(u2)
Pripravimo podatke in uporabimo funkcijo "lambdify" za pretvorbo iz simbolnega zapisa v numeričnega.
podatki = {L: 1e-3,C: 1e-5,R: 2, Ug:3}
tau=L/R # časovna konstanta
tau2=tau.subs(podatki) # numeričen izračun časovne konstante z zamenjavo oz. vstavitvijo vrednosti za L in R
print('Tau = ',tau2,' s')
print(podatki)
tok = lambdify(t, i2.args[1].subs(podatki), 'numpy') # tok kot funkcija časa
napetost = lambdify(t,u2.subs(podatki), 'numpy') # napetost kot funkcija časa
print('Tok pri 0s: {:g} A'.format(tok(0)))
print('Napetost pri 0s: {:g} V'.format(napetost(0)))
Tau = 0.000500000000000000 s {L: 0.001, 5e-05: 1e-05, R: 2, U_g: 3} Tok pri 0s: 0 A Napetost pri 0s: 3 V
Sledi priprava enačb za izrise. Naredimo dva niza: cas za izpis gladke krivulje z več točkami in cas2 za izris simbolov v redkih točkah:
cas = np.linspace(0, 1e-3, 100)
cas2 = np.linspace(0, 1e-3, 5)
def slika():
plt.plot(cas, tok(cas), 'C0', label='Tok [A]')
plt.plot(cas, napetost(cas), 'C1', label='Napetost [V]')
plt.plot(cas2, tok(cas2), 'C0o', label='Tok - velik korak')
plt.plot(cas2, napetost(cas2), 'C1o', label='Napetost - velik korak')
plt.xlabel('Čas [s]')
plt.ylabel('Tok [A] / Napetost [V]')
plt.legend(loc=(1.01, 0));
plt.show()
slika()
Še mali programerski izziv za zagnane: Dodaj vrstice kode, s katerimi bi določil časovno konstanto iz analize signala. Npr. tako, da poiščeš čas, ko znaša napetost $u(t=\tau) = U_g e^{-1}$. (Namig: 1. tvori niz vrednosti napetosti (kot je narejeno v plotu). 2. izračunaj vrednost napetosti pri $\tau = RC$. To napetost moraš poiskati v nizu oziroma vrednost, ki je njej najbližje. 3. V ta namen uporabi funkcijo np.where, s katero poiščeš vrednosti večje ali manjše od določene vrednosti. Rezultat je niz vrednosti, ki zadostujejo temu pogoju. 4. Največja vrednost v nizu je tisti indeks, ki ga iščeš. 5.Sedaj le še vstaviš indeks v niz cas in imaš rešitev. Lažje reči kot narediti? Rešitev je na koncu strani.
Rešujemo enačbo ${{u}_{g}}=RC\frac{\text{d}{{u}_{C}}}{\text{d}t}+{{u}_{C}}$, kjer je ${{u}_{g}}(t)={{U}_{g}}\sin {(\omega t})$.
u = Function('u') # pripravimo funkcijo in spremenljivke ter zapišemo dif enačbo
C, R, t, Ug, omega, u0 = symbols('C, R, t, U_g, omega, u0')
eq = Eq(R*C*u(t).diff(t), Ug*sin(omega*t)-u(t))
eq
Najprej poskušamo identificirati tip diferencialne enačbe. V ta namen uporabimo funkcijo classify_ode.
classify_ode(eq, u(t))
('1st_linear', 'Bernoulli', 'almost_linear', '1st_power_series', 'lie_group', 'nth_linear_constant_coeff_undetermined_coefficients', 'nth_linear_constant_coeff_variation_of_parameters', '1st_linear_Integral', 'Bernoulli_Integral', 'almost_linear_Integral', 'nth_linear_constant_coeff_variation_of_parameters_Integral')
Preiskusimo lahko vse zgoraj naštete možne tipe. V našem primeru dobro deluje nth_linear_constant_coeff_undetermined_coefficients
dsol = dsolve(eq, u(t),'nth_linear_constant_coeff_undetermined_coefficients') # rešimo diff enačbo
display(dsol)
Podobno kot v prejšnjem primeru, zamenjamo konstanto $C_1$ za začetni pogoj, $u(t=0)$:
t0=dsol.args[1].subs({'t':0})
eq_init = Eq(u0, t0) # izrazimo enačbo pri t=0
display(eq_init) # prikažemo i0
init_solve=solve(eq_init, C1) # rešimo enačbo za C1
display(init_solve) # prikažemo rešitev
final = dsol.subs(C1, init_solve[0]) # v rešitev vstavimo C1 in dobimo končno enačbo
final
Vstavimo začetni pogoj in dobimo končno rešitev:
zacetni_pogoji = {'u0': 0}
u2=final.subs(zacetni_pogoji)
display((u2))
Izračunamo še tok kot $i=C\frac{\text{d}{{u}_{C}}}{\text{d}t}$
i2=symbols('i2')
i2=C*u2.rhs.diff(t) # množenje s C in odvodom desne strani enačbe
simplify(i2)
Sedaj lahko ocenimo še pravilnost enačbe. Čeprav je zapis zapleten, lahko ocenimo rešitev pri $CR\omega >>1$. Najprej izbrišemo 1 in $cos(ut)$, nato še eksponenta, ki se izničita. Ostane nam
$\frac{C U_g \omega }{C^{2} R^{2} \omega^{2}} \left(\left(C R \omega \sin{\left (\omega t \right )}\right) \right)$,
od koder sledi
$\frac{U_g }{R} \left(\left(\sin{\left (\omega t \right )}\right) \right)$,
kar je tok, ki ga diktira upornost $R$.
Če pa je $CR\omega << 1$, pa se ohrani cosinusni člen, eksponenta se za velike čase izničita, ostane $\frac{C U_g \omega }{C^{2} R^{2} \omega^{2}+1} \left(\left(\cos{\left (\omega t \right )}\right) \right)$ in nato $\frac{C \omega U_g}{1} \left(\left(\cos{\left (\omega t \right )}\right) \right)$. V tem primeru, npr. nizkih frekvencah, tok skozi vezje diktira "upornost" kondenzatorja, ki jo pri izmeničnih signalih imenujemo reaktanca in je enaka $\frac{1}{\omega C}$.
podatki = {C: 5e-5, R: 5e5, Ug: 2, omega:1000}
tau=R*C # časovna konstanta
tau2=tau.subs(podatki) # numeričen izračun konstante z zamenjavo oz. vstavitvijo vrednosti za L in R
print('Tau = ',tau2,' s')
napetost = lambdify(t,u2.args[1].subs(podatki), 'numpy') # napetost kot funkcija časa
tok = lambdify(t, i2.subs(podatki), 'numpy') # tok kot funkcija časa
print('Tok pri 0s: {:g} A'.format(tok(1)))
print('Napetost pri 0s: {:g} V'.format(napetost(1)))
Tau = 25.0000000000000 s Tok pri 0s: 3.30745e-06 A Napetost pri 0s: 3.18755e-05 V
cas = np.linspace(0, 1, 50)
def slika():
plt.plot(cas, tok(cas), 'C0', label='Tok [A]')
plt.plot(cas, napetost(cas), 'C1', label='Napetost [V]')
plt.ylabel('Tok [A] / Napetost [V]')
plt.legend(loc=(1.01, 0));
plt.show()
slika()
Zgornji prikaz ni napačen, problem je le, da sta lahko napetost in tok različnih velikostnih razredov. Zato je bolj primerno izrisati tok in napetost na dveh različnih grafih ali pa na enem grafu ampak z dvema različnima ordinatnima osema. Poleg tega časovni korak ni izbran najbolj ustrezno in ga je potrebno 2zgostiti". Glej spodaj:
cas = np.linspace(0, 2e-2, 5000)
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('Čas (s)')
ax1.set_ylabel('Tok (A)', color=color)
ax1.plot(cas, tok(cas), color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx() # x os za drugi plot naj bo enaka x osi prvega plota
color = 'tab:blue'
ax2.set_ylabel('Napetost (V)', color=color)
ax2.plot(cas, napetost(cas), color=color)
ax2.tick_params(axis='y', labelcolor=color)
plt.show()
Še mali programerski izziv za zagnane: Dodaj vrstice kode, s katerimi bi določil časovno konstanto iz analize signala. Npr. tako, da poiščeš čas, ko znaša napetost $u(t=\tau) = U_g e^{-1}$. (Namig: 1. tvori niz vrednosti napetosti (kot je narejeno v plotu). 2. izračunaj vrednost napetosti pri $\tau = RC$. To napetost moraš poiskati v nizu oziroma vrednost, ki je njej najbližje. 3. V ta namen uporabi funkcijo np.where, s katero poiščeš vrednosti večje ali manjše od določene vrednosti. Rezultat je niz vrednosti, ki zadostujejo temu pogoju. 4. Največja vrednost v nizu je tisti indeks, ki ga iščeš. 5.Sedaj le še vstaviš indeks v niz cas in imaš rešitev. Lažje reči kot narediti? Rešitev je na koncu strani.
Več:
#Rešitev programerskega izziva:
nap=napetost(cas)
display(np.max(nap)*np.exp(-1))
xx=np.where(nap >= np.max(nap)*np.exp(-1))
x2=((np.max(xx)))
display(cas[x2])