Harj10LV
22.11.02
Alustukset
> | restart:with(inttrans):alias(L=laplace,IL=invlaplace,u=Heaviside): |
1.
> | restart: with(inttrans): alias(LL=laplace,IL=invlaplace,u=Heaviside): # Joudumme käyttämään LL-nimea, koska L on tässä tehtävässä varattu kelan induktanssille. |
> | RC:=int(i1(tau),tau=0..t)/C+R*i2(t)=e(t); |
> | RL:=L*diff(i1(t),t)-L*diff(i2(t),t)=R*i2(t); |
> | L:=2;C:=50*10^(-6);R:=100;e:=t->50*sin(100*t); |
> | LRC:=LL(RC,t,s); |
> | LRL:=LL(RL,t,s); |
> | AE:=i1(0)=0,i2(0)=0; |
> | LRC:=subs(AE,LRC);LRL:=subs(AE,LRL); |
> | solve({LRC,LRL},{LL(i1(t),t,s),LL(i2(t),t,s)}); |
> | I12:=subs(%,[LL(i1(t),t,s),LL(i2(t),t,s)]); |
> | I1:=I12[1]: I2:=I12[2]: |
> | I1:=convert(I1,parfrac,s);I2:=convert(I2,parfrac,s); |
> | i1:=IL(I1,s,t); |
> | i2:=IL(I2,s,t); |
> | plot([e(t)/R,i1,i2],t=0..3/10,color=[gray,red,blue],legend=["Syöttöjännite normeerattuna R:llä","Virta i1","Virta i2"]); |
2.
> | restart:with(inttrans):alias(L=laplace,IL=invlaplace,u=Heaviside): |
> | dy1:=m1*diff(y1(t),t,t)=-k1*y1(t)+k2*(y2(t)-y1(t)); dy2:=m2*diff(y2(t),t,t)=-k2*(y2(t)-y1(t))-k3*y2(t); |
> | AE1:=y1(0)=0,D(y1)(0)=1; AE2:=y2(0)=0,D(y2)(0)=-1; |
> | Ldy1:=L(dy1,t,s): Ldy1:=subs(AE1,Ldy1); |
> | Ldy2:=L(dy2,t,s): Ldy2:=subs(AE2,Ldy2); |
> | YY:=solve({Ldy1,Ldy2},{L(y1(t),t,s),L(y2(t),t,s)}); |
> | Y1:=subs(YY,L(y1(t),t,s)); Y2:=subs(YY,L(y2(t),t,s)); |
Maplella kun lasketaan, niin voidaan kivuttomasti viedä tähän pisteeseen yleisin symbolein. Tässä nyt viimeistään täytyy sijoittaa arvot.
> | m1:=10: m2:=10: k1:=20: k2:=20: k3:=40: |
> | Y1; Y1:=normal(Y1); # Supistaminen onnistuus sekä komennolla normal ... |
> | Y2; Y2:=simplify(Y2); # ... että simplify. |
Nimittäjä pitää jakaa tekijöihin. (Tarvitaan sitä kuuluisaa 2. asteen yhtälön ratkaisukaavaa!)
> | nimittaja:=denom(Y1); |
> | juuret:=solve(nimittaja=0,s); |
> | nim1:=expand((s-juuret[1])*(s-juuret[2])); |
> | nim2:=expand((s-juuret[3])*(s-juuret[4])); |
Tästä päästään kiinni osamurtohajoitelmaan. Annetaan tässä nyt Maplen laskea.
Maple ei selviä ilman tietoa siitä, että kertoimien joukossa esiintyy . Tämän ilmoitamme näin:
> | Y1:=convert(Y1,parfrac,s,sqrt(5)); |
> | Y2:=convert(Y2,parfrac,s,sqrt(5)); |
Tästä onnistuu käänteismuuntaminen neliöksi täydentäen, s-siirtolausetta ja cos/sin-muunnoksia soveltaen.
Annetaan Maplen hoitaa yksityiskohdat:
> | y1:=IL(Y1,s,t); |
> | y2:=IL(Y2,s,t); |
> | y1:=convert(y1,exp):evalc(%):simplify(%); |
> | y2:=convert(y2,exp):evalc(%):simplify(%); |
> | y2:=collect(%,{sin((5+5^(1/2))^(1/2)*t),sin((5-5^(1/2))^(1/2)*t)}); |
> | y2:=map(simplify,y2); |
> | plot([y1,y2],t=0..10); |
3.
> | restart:with(inttrans):alias(L=laplace,IL=invlaplace,u=Heaviside): |
> | dy:=diff(y(t),t,t)+2*diff(y(t),t)+y(t)=Dirac(t)-u(t-2*Pi); |
> | AE:=y(0)=0,D(y)(0)=0: |
> | Ldy:=L(dy,t,s); |
> | Ldy:=subs(AE,Ldy); |
> | Y:=solve(Ldy,L(y(t),t,s)); |
> | nim:=factor(denom(Y)); |
> | Y1:=s/nim; |
> | y1:=IL(Y1,s,t); # |
:n muunnos ja s-siirto (s - (-1)) .
> | Y2perexp:=-1/nim; |
> | Y2perexp:=convert(%,parfrac,s); |
Keskimmäinen jo muunnettiin, muut suoraan:
> | f2:=IL(Y2perexp,s,t); |
aiheuttaa t-siirron
> | y2:=subs(t=t-2*Pi,f2)*u(t-2*Pi); |
> | y:=y1+y2; |
> | L(y,t,s): simplify(%); Y; # Tarkistus, oikein meni. |
> | plot([-u(t-2*Pi),y],t=0..25,numpoints=100,axes=boxed); |
> | plots[display](%,view=[5..8,-0.4..0.1]); |
> |
4.
> | restart:with(inttrans):alias(L=laplace,IL=invlaplace,u=Heaviside): |
> | konv:=(f,g,t)->int(f(tau)*g(t-tau),tau=0..t); |
a)
> | F:=1/(s-a); |
> | f:=t->exp(a*t); |
> | tulos:=konv(f,f,t); |
> | L(%,t,s); #Tarkistus. |
b)
> | F1:=1/s^2; F2:=1/(s^2+omega^2); |
> | f1:=IL(F1,s,t); f2:=IL(F2,s,t); |
> | f1:=unapply(f1,t): f2:=unapply(f2,t): |
> | f:=konv(f1,f2,t); # Tämä on siis kysytty käänteismuunnos: |
> | L(f,t,s);simplify(%); # Tarkistus, ok. |
c)
> | F1:=s/(s^2+omega^2) ; F2:=1/(s^2+omega^2); |
> | f1:=unapply(IL(F1,s,t),t); |
> | f2:=unapply(IL(F2,s,t),t); |
> | f:=konv(f1,f2,t); # Tässä vastaus. |
> | L(f,t,s); #Tarkistus, oikein meni. |
5.
> | restart:with(inttrans):alias(L=laplace,IL=invlaplace,u=Heaviside): |
> | dy:=diff(y(t),t,t)+omega^2*y(t)=r(t); |
> | AE:=y(0)=y0, D(y)(0)=v0; |
> | Ldy:=L(dy,t,s); |
> | Ldy:=subs(AE,Ldy); |
> |
> | Ldy:=subs(L(r(t),t,s)=R(s),Ldy); |
> | Y:=solve(Ldy,L(y(t),t,s)); |
> | Ytr:=(s*y0+v0)/(s^2+omega^2); Yss:=(R(s))/(s^2+omega^2); |
> | ytr:=IL(Ytr,s,t); |
> | Q:=1/(s^2+omega^2); |
> | q:=unapply(IL(Q,s,t),t); |
> | yss:=konv(q,r,t); |
> | y:=ytr+yss; |
Kaunis ratkaisukaava.
6.
> | restart: with(inttrans): alias(L=laplace,IL=invlaplace,u=Heaviside): plots[setoptions](axes=boxed): |
> | interface(showassumed=2): |
a) KRE Exa 4 ss. 271-272
> | dyht:=diff(y(t),t,t)+3*diff(y(t),t)+2*y(t)=r(t); |
> | r:=t->u(t-1)-u(t-2); |
> | plot(r,0..3,title="Herätteenä neliöaalto, välin [1,2] karakt. fkt",scaling=constrained); |
> | AE:=y(0)=0,D(y)(0)=0; |
> | Ldyht:=L(dyht,t,s); |
> | Ldyht:=subs(AE,Ldyht); |
> | Y:=solve(Ldyht,L(y(t),t,s)); |
> | F:=1/denom(Y); es:=numer(Y); |
> | F:=convert(F,parfrac,s); |
Tästäpä on helppo käänteismuuntaa t-siirtolauseen avulla. Antaa nyt Maplen hoitaa:
> | yy:=IL(es*F,s,t); |
> | plot([r(t),yy],t=0..6,legend=["Heräte","Vaste"]); |
> |
b)
Lasketaan pisteeseen t=1 keskitetyn -jonon vaste. Ts. otetaan jono suorakulmiopulsseja, joiden ja .
eli (jos kyseessä on voima). Olemme käyttäneet luennolla ja www-viitteissämme merkintää . Tässä siis .
Määritellään mieluummin symmetrisesti samalla tavoin kuin luennolla (eikä KRE-tavalla).
> | d1:=(1/epsilon)*(u(t-1+epsilon/2)-u(t-1-epsilon/2)); |
> | plot(subs(epsilon=0.2,d1),t=0..2,scaling=constrained); |
> | Ld1:=1/(2*epsilon*s)*(exp(-s)*exp(epsilon*s/2)-exp(-epsilon*s/2)); |
> | dyht:=diff(y(t),t,t)+3*diff(y(t),t)+2*y(t)=d1; |
> | AE:=y(0)=0,D(y)(0)=0; |
> | assume(epsilon >0 , epsilon<1): # Muuten Maple ei L-muunna oikeaa puolta (eikä pidäkään). |
> |
> | Ldyht:=L(dyht,t,s); |
> | Ldyht:=subs(AE,Ldyht); |
> | Y:=solve(Ldyht,L(y(t),t,s)); |
> | y:=IL(Y,s,t); |
Tästä voitaisiin määrätä raja-arvo, kun -> samaan tapaan kuin tehtiin yleisesti (otetaan exp-funktion sarjasta sopiva määrä termejä).
Demonstroimme numeerisesti ja kuvilla tilannetta.
> | y1:=subs(epsilon=1,y):y05:=subs(epsilon=0.5,y):y02:=subs(epsilon=0.2,y):y01:=subs(epsilon=0.1,y):y005:=subs(epsilon=0.05,y): |
> | plot([y1,y05,y02,y01,y005],t=0..4,color=[red,blue,black,green,cyan,gold],title="Systeemin vaste kapeneviin yksikköpulsseihin",legend=["eps=1","eps=0.5","eps=0.2","eps=0.1","eps=0.05"]); |
> |
> |
> |
> | F:=1/(s*(s+1)*(s+2)); |
c)
Lasketaan nyt sama lasku, kun herätteenä on yksikköimpulssi hetkellä t=1, eli käytetään suoraan Diracin deltaa herätteenä.
Merkitään nyt ratkaisua x(t):llä ja vastaavasti L.muunnosta X(s):llä.
> | dy:=diff(x(t),t,t)+3*diff(x(t),t)+2*x(t)=Dirac(t-1); |
> | AE:=x(0)=0,D(x)(0)=0;Ldy:=L(dy,t,s); |
> | Ldy:=subs(AE,Ldy); |
> | X:=solve(Ldy,L(x(t),t,s)); |
> | x:=IL(X,s,t); |
Voitaisiin tietysti käyttää tätä muotoa, mutta saadaan sievempään (ja tutumpaan) muotoon näin:
> | convert(x,exp);expand(%);simplify(%,symbolic);x:=%; |
> | plot(x,t=0..20); |
> | plot([y1,x],t=0..4,color=[blue,gold]); |
> | plot([y02,x],t=0..4,color=[black,gold]); |
> | plot(y005-x,t=0..3); |
> |
b)-kohdan jono näkyy lähestyvän delta-herätettä vastaavaa tulosta. Tämä ei ole yllättävää. deltan Laplace-muunnos määrättiin juuri
tällaisen rajaprosessin avulla (L-muunnosten raja-arvo, muunnettaessa -jonoa).