Maple-ohjeita Y3 syksy -96

Harj. 12

Teht. 1

#Olk. n = x-osavälien lkm, p=n-1=sisäsolmujen lkm (x-suunnassa)
#     kmax= aikaindeksin suurin arvo
n:=5:p:=n-1:kmax:=5:u:='u':h:=0.2:
# Alkuehdot:
for j from 0 to n do u[j,0]:=evalf(sin(j*Pi*h)) od;
# Reunaehdot:
for k from 0 to kmax do u[0,k]:=0: u[n,k]:=0 od:
#
# Muodostetaan kutakin k-arvoa kohti C-N-yhtälö (r=1)ja ratkaistaan samantien.
# Saadut arvot ovat tunnettuja seuraavassa vaiheessa. Tässä on hyvin
# kätevää käyttää assign-komentoa, joka muuttaa solve-komennolla saatavat
# yhtälöt sijoituslauseiksi (eli voidaan ajatella, että "=" muuttuu ":=" -
# operaatioksi)
#
for k from 0 to kmax do
     for j from 1 to p do
        yht[j,k]:=4*u[j,k+1]-u[j+1,k+1]-u[j-1,k+1]=u[j+1,k]+u[j-1,k]
                       od:
         assign(solve({seq(yht[j,k],j=1..p)},{seq(u[j,k+1],j=1..p)}))
                     od:
# Järjestetään tulokset matriisiksi:
with(linalg):
umatr:=matrix(n+1,kmax+1,[seq(seq(u[j,k],j=0..n),k=0..kmax)]):
op(umatr);

Harj. 11

Teht. 1

n:=5:h:=2/n:x5:=[seq(-1+k*h,k=0..n)];  # väli [-1,1], pituus = 2
# Sitten voit muutella n:ää.
f:=x->1/(1+25*x^2);
y5:=map(f,x5);                         # Tämä on kätevä ja yleispätevä tapa
p5:=interp(x5,y5,x);
plot({f(x),p5},x=-1..1);

Huom! Erilaisiin funktioiden approksimointeihin on Maplessa kiintoisan näköinen pakkaus numapprox (kts. ?numapprox) Esim. pade on usein käytännössä hyödyllinen, se muodostaa rationaalifunktio- approksimaation.

Teht. 2

 
# alkuviikon teht. 5 vastaus:
> u:=100*sin(Pi*x/80)*exp(-(1.158*Pi/80)^2*t)
# Pintapiirros:
>  plot3d(u, x=0..80, t=0..500,axes=boxed, title=`tehtava 2`);

> # lämpötilaprofiileja: annetut ajanhetket t=0.1 jne käsitetään
> minuuteissa annetuiksi, muuten käyrät ovat turhan lähekkäin.  
> Eli 0.1 <-> 6, 0.5 <-> 30 jne.
> plot({subs(t=0,u),subs(t=6,u),subs(t=30,u)},x=0..80):

Kuvaa voi kierrellä tarttumalla hiiren vasemmalla ja
sitten oikella valita "redraw"

Jos mieluummin teet display-tyylillä, niin mikäs siinä:

> p1:=plot(subs(t=0,u),x=0..80):
> p2:=plot(subs(t=6,u),x=0..80):
> ...
> with(plots):
> display({p1,p2,p3,p4,p5});

Teht. 3

  u:='u':n:='n';  # Varmuuden vuoksi
> osdy:=a^2*diff(u(x,t),x$2)=diff(u(x,t),t);
>  eval(subs(u(x,t)=X(x)*T(t), osdy));
 # Jaetaan puolittain (a^2*X(x)*T(t)):llä:
> separoitu:=simplify("/(a^2*X(x)*T(t)));
> xyht:=lhs(separoitu)=-p^2;   # Vasen puoli = vakio (negat)
> tyht:=rhs(separoitu)=-p^2;   # Oikea puoli = sama vakio
> xyht:=xyht*X(x);  tyht:=a^2*tyht*T(t); # Kerrotaan nimittäjät pois
> xx:=rhs(dsolve(xyht,X(x)));
> dsolve(tyht,T(t)); tt:=rhs(");                                    
> tt:=subs(_C1=1,tt);  # (Vakio kertaa vakio) yhdistetään yhdeksi vakioksi.
                       # Katso, onko _C1 vai _C2, dsolve voi vaihdella
                       # järjestystä.
> # RE1 => ... RE2 => ..., sijoita nämä suoraan (Tässä ei ole mieltä "solvata")
> # sin:n 0-kohdat ovat Pi:n monikertoja, kuten tiedetään.
> xx:=subs({_C1=0,_C2=Bn,p=...},xx);
> tt:=subs(p=...,tt);
> uu:=xx*tt

# Sitten vielä alkuehto:
# Sinisarjan kertoimet:

> Bn:=simplify((2/l)*int(F*sin(n*Pi*x/l),x=0..l)); # Tässä F on alkuarvofunktio
>                                             # Sijoita se F:n paikalle.
Huom! Maple ei havaitse, että n:n arvolla 1 tulee 0. No sen voi
ottaa huomioon alkamalla alempana olevan summauksen 2:sta

> l:=...;a:=...; # Sijoita annetut numeeriset arvot
> UU:=sum(uu,n=1..10);  # tai sum:n sijasta  add  (vm. on kai nykyisin
                        #                                           suositus)
                        # Otetaan siis 10 termiä (voit kokeilla muutakin)
> plot3d(UU, x=0..10, t=0..5,axes=boxed);

plot({subs(t=0,UU),xsubs(t=0.1,UU),subs(t=0.5,UU),subs(t=2,UU),subs(t=5,UU},
x=0..10):

 Teht. 4 

Painovirhe: p.o.: a2=... (eli potenssi 2 on pudonnut
tehtäväpaperista)
-----------------------------------------------------

Maple-alustuskomentoja

with(inttrans)  # V.4  # readlib(laplace)  #V.3
alias(H=Heaviside,L=laplace,IL=invlaplace):
############
### Muistathan vielä, jos matriiseja tarvitset:
with(linalg):

Harj. 10

Periaate:

Kehitetään oikea puoli r(x) Fourier-sarjaksi r(x)=r1(x)+r2(x)+...
ja etsitään erikseen diff. yhtälöiden

vp=r1(x), vp=r2(x) ,...

erityisratkaisut käyttämällä yritteitä yr1,yr2,...

Tällöin yr1+yr2+...  on yhtälön  vp=r1(x)+r2(x)+...  ratkaisu.

Harj. 9 ja 10

Kannattaa käsitellä vasen puoli (vp) ja oikea puoli (op) erikseen.

with(inttrans):alias(L=laplace,IL=invlaplace,H=Heaviside):
vp:=R*diff(i(t),t)+(1/C)*i(t);
Lvp:=L(vp,t,s);
Lvp:=subs(i(0)=0,L(vp,t,s));

Lop:=int((e0/epsilon)*exp(-s*t),t=0..epsilon); # Lasketaan suoraan L-muunnoksen
                                               # määritelmästä.
Ldy:=Lvp=Lop     # Laplace-muunnettu yhtälö
II:=solve(Ldy,L(i(t),t,s));

assume(a>0):L(H(t-a),t,s);
assume(epsilon>0);#  Tätä tarvitaan, jotta vastaavasti käänteismuunnokset
                  # jatkossa onnistuisivat.
a:='a':     # Tällä pääsee eroon "assume-ominaisuudesta".

Fourier-sarjoista Luentoesimerkki

Otetaan vielä pelkät inputit (valikoiden/modifioiden)

Eulerin Fourier-kerroinkaavat: 

## Jätetään Maple-prompti pois, silloin voi leikata/liimata kokonaisuuden.

 a0:=(1/(2*l))*int(f(x),x=-l..l);  ## Tässä ei ole mitään olennaista
 an:=(1/l)*int(f(x)*cos(n*Pi*x/l),x=-l..l); # etua a[n]-tyylistä, se
 bn:=(1/l)*int(f(x)*sin(n*Pi*x/l),x=-l..l); # voi pikemminkin olla
                                             # harhaanjohtava.

## Nyt voidaan laskea mitä esimerkkejä halutaan. Vaikkapa harj.9 AV teht. 6:

 f:=x->x;
 assume(n,integer):
 a0;an;
 bn; # V.4 osaa käyttää integer-tietoa, V.3 ei osaa - no tehdään sitten itse:
 bn:=subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n,l=Pi},bn);

 s:=seq(sum(bn*sin(n*x),n=1..N),N=1..5);
 plot({x,s},x=-Pi..Pi);

## Jos halutaan laskea jokin erillinen osasumma vaikka suurella n:llä
## Gibbsin ilmiön havainnollistamiseksi tai muuten, niin luonnollisesti:

 s50:=sum(bn*sin(n*x),n=1..50):plot({x,s50},x=-Pi..Pi);
 plot(x-s50,x=-Pi+.1..Pi-.1); #Erotus,mutta ei ihan "hyppypisteeseen" saakka

Harj. 8

Teht. 1 ja 2

Teht. 3-6

• exp -funktio lasketaan exp(t) eikä e^t. Maple ei tunne symbolia e (ellei tehdä e:=exp(1), mitä en suosittele (Aivan sama pätee Matlabiin nähden.)

> with(inttrans); #Versiossa V.4: 
> readlib(laplace);#Versiossa V.3:
> alias(H=Heaviside,L=laplace,IL=invlaplace); # Tämä on aika mukava

Luentoesimerkkejä L-muunnoksen käytöstä

Harj. 5

Teht. 1

dyht:=diff(y(t),t)=sin(t)-y(t)/2;

ratk:=dsolve(dyht,y(t));
Y:=rhs(ratk);  # dsolve palauttaa yhtälön, otetaan sen oikea puoli (rhs).
Y:=subs(_C1=c,Y); # helpompi symboli vakiolle
parvi:=seq(Y,c=-10..10): # Voit kokeilla eri c:n arvoja ja
plot({parvi},t=0..10);   #   aikaskaalaa (huom: tämä vaikuttaa hyvin
                         #   olennaisesti kuvaan)

> 
> # Teht. 4
> 
> with(linalg):
> A:=matrix(2,2,[3,-1,-1,3]);
> esys:=eigenvects(A);
> v1:=op(esys[1][3]);v2:=op(esys[2][3]);
> lam:=eigenvals(A);
> y:=evalm(c1*exp(lam[2]*t)*v1+c2*exp(lam[1]*t)*v2);
> # Katsotaan pääakselikoordinaatistossa:
> z1:=exp(lam[2]*t);z2:=exp(lam[1]*t);
> # z2=K*z1^2
> # Piirretään ensin z1-z2-koordinaatistoon (sanotaan (u,v))
> 
> v:=seq(K*u^2,K=-10..10);
> plot({v},u=-2..2);
> op(v1);op(v2);
> y1aks:=plot([[0,0],convert(v1,list)]):y2aks:=plot([[0,0],convert(v2,list)]):
> with(plots):

> parvi:=seq(seq([y[1],y[2]],c1=-5..5),c2=-5..5):

> varit:=[aquamarine, black , blue ,     navy,   coral, cyan , \
> brown ,     gold ,  green ,    gray ,  grey , khaki ,\
> magenta   , maroon, orange   , pink ,  plum  ,red,\
> sienna ,    tan  ,  turquoise, violet, wheat ,white ,yellow];
>  
> ## pp:=seq(plot([op(parvi[i]),t=-1.. .4],color=varit[i]),i=1..25):
> pp:=seq(plot([op(parvi[i]),t=-1.. .4],y1=-5..5,y2=-5..5,color=varit[i]),i=1..25):
> display({pp,y1aks,y2aks});

Alkuviikko, teht. 3

Harj. 4

Teht. 2

A:=randmatrix(20,20);
A:=convert(op(A),float);     # jotta eigenvals laskisi numeerisesti
eig:=eigenvals(A);
c:=seq([Re(eig[i]),Im(eig[i])],i=1..20);  # [x1,y1],[x2,y2],...
plot([c],style=point);       # Piirtää listassa olevat pisteet
r:=max(op(map(abs,[eig])));  # spektraalisäde

-------------------
plot([r*cos(t),r*sin(t),t=0..2*Pi]):  # parametrikäyrien syntaksi
plot([c],x=-100..100,y=-100.100,style=point): # x-ja y-välit kannattaa antaa
                                              # ainakin Re- tai Im-akselilla
                         # sijaitsevia pistejoukkoja piirrettäessä
Kuvat saadaan samaan näin:

with(plots): # Ladataan lisägrafiikkapakkaus, jossa mm. display-funktio
p1:=plot([r*cos(t),r*sin(t),t=0..2*Pi]):  # muista kaksoispiste, eikä 
p2:=plot([c],x=-100..100,y=-100.100,style=point): # puolipistettä (tai kärsi)
display({p1,p2}):

Harj. 2

Teht. 1,3,4,5

> with(linalg):   # Aina, kun käsitellään matriiseja
> A:=matrix([[rivi1],[rivi2],...[rivim]]);
> ##   tai
> A:=matrix(3,4,[a11,a12,a13,a14,a21,a22,...,a34]); # koko ja sitten data
>                            # pötkönä vaakarivijärjestyksessä.
> b:=vector([b1,b2,b3]);
> Ab:=augment(A,b);  # Vektori liitetään matriisin perään.
> # Voidaan liittää myös matriisi, kuten teht. 5:
> Id:=diag(1,1,1); AI:=augment(A,Id);

Miten poimitaan matriisin osia, miten liitettäisiin alle

> osa:=submatrix(A,i1..i2,j1..j2)
> osa2:=submatrix(A,[i1,i2,i3],[j1,j2,j3,j4]); # Sieltä täältä poiminta
> stack(A,b)  # alle liittäminen (näissä tehtävissä ei tarvita)

Teht. 6 , GramSchmidt, normeeraus, map

Miten sitten operoidaan listaan vektoreita? normalize ei suoraan sovella itseään alkioittain. Tehdään se yleispätevällä komennolla map

Esim:

> lista:=[a,b,c];

                               lista := [a, b, c]
------------------------------------------------------------------------------
> map(f,lista);

                               [f(a), f(b), f(c)]
> 

Tässä tapauksessa

> gslista:=GramSchmidt([v1,v2,v3]);  # ensin tietysti v1:=...,v2:=...,v3:=...;
> ortoNORMAALIkanta:=map(normalize,gslista);  # kaunista!

Siis, kun haluat itse päättää järjestyksen, niin käytä listaa , eli hakasulkuja [...] .