![[Up]](/gadgets/trup.gif)
cat > oma.txt Maalaus hiiren vasemmalla Liimaus hiiren keskimmäisellä CTR-D tiedoston sulkeminen(Tai emacs tms.)
Sitten Maple-istunnossa:
> read("oma.txt");
> read("c:\\windows\\mymaple\\oma.txt"); (PC-WIN:ssä)
Huom: Joskus näyttää leikkaus/liimaus toimivan emacsista, vaikka
www-sivulta ei toimi.
Kaikki kurssilla talletetut ohjelmat voi lukea näin: (yritetään pitää ajan tasalla)
>read("/p/edu/mat-1.414/maple/v201.mpl");
>read("/home/apiola/opetus/peruskurssi/v2-3/maple/v201.mpl");
Jälimmäinen on HA:n privaatti.
v2l:=vek->convert(vek,list): # HAM-assa painovirhe (vec) l2v:=lis->convert(lis,vector):Liitetään kaksi vektoria peräkkäin:
vappend:=(v1,v2)->vector([op(convert(v1,list)),op(convert(v2,list))]):Matriisin jonouttaminen pitkäksi vektoriksi
m2l:=M->map(op,convert(M,listlist)): # Vastaa Matlabin M(:), mutta tässä
# riveittäin
linspace:=proc() local i,n,a,b;
a:=args[1];b:=args[2];
if nargs=2 then n:=10 else n:=args[3] fi;
[seq(a+i*(b-a)/(n-1),i=0..n-1) ] end:
Esim:
x10:=linspace(0,1); x20:=linspace(0,1,20);
Kts. 5.3 Diskretointi [HAM] s. 117. Vaihdoimme oletusarvon 10:een,
Matlab:n 100:n sijasta.
Kompleksiluvut ja -funktiot
Konformikuvauksia
Tässä on
ratkaisutyöarkin harj1.mws lopusta
poimittuja kuvausvälineitä ja esimerkkiskriptejä hieman täydennetyin
selityksi, mutta muuten tiiviimmässä muodossa.
>
# Vaihtoehtoisia ja paranneltuja Maple-ratkaisutapoja
# Aikalailla kätevää on tehdä "cobweb"-tyylillä. Saamme varsin
# yleispäteviä työkaluja. Tarkoittaa sitä, että määrittelemme muutaman
# grafiikka-arvoisen funktion.
# Sovitaan ne g-alkuisiksi (ainakin nyt tässä).
#
# gympyra on funktio, jolle annetaan keskipiste ja säde , tuloksena on
# vastaava ympyrägrafiikka.
# gympyrakuva on vastaava kuva sovellettaessa funktiota f, joka myös
# annetaan argumenttina.
> gympyra:=(z0,r)->complexplot(z0+r*exp(I*Theta),Theta=0..2*Pi):
> gympyrakuva:=(z0,r,f)->complexplot(f(z0+r*exp(I*Theta)),Theta=0..2*Pi)
> :
# Testataanpa. Katsotaan, mitä exp-funktio tekee 1-ympyralle
#
> display([gympyra(0,1),gympyrakuva(0,1,exp)]);
# Katsotaan nyt se Joukowski-kuvaus ja nähdään toden totta, että
# reaaliakselin väli [-2,2] on kuvana.
> display([gympyra(0,1),gympyrakuva(0,1,z->z+1/z)]);
# Voidaan tietysti kokeilla kaikenlaisia funktioita. Tässä on tapa
# sijoittaa kuvat vierekkäin, tehdään yleinen taulukko array .
# (Kuten matriisi, mutta sallii kaikenlaisia olioita alkioikseen.)
>
> f:=z->sin(z):display(array([[gympyra(0,1),gympyrakuva(0,1,z->z+sin(1/z
> ))]]));
>
# Ympyräparvet syntyvät käden käänteessä, pannaan samalla liikettä
# mukaan:
> display([seq(gympyra(I+0.1*k,1),k=0..20)],insequence=true);
> display([seq(gympyra(0,0.1*k),k=0..10)]);
> gsade:=(Theta,r1,r2)->complexplot(r*exp(I*Theta),r=r1..r2):
> gsadekuva:=(Theta,r1,r2,f)->complexplot(f(r*exp(I*Theta)),r=r1..r2):
> display([seq(gsade(k*2*Pi/10,0.5,3),k=0..9)]);
Mielivaltaisen ympyrärenkaiden ja sektorisäteiden verkon kuvaaminen
Tässä on pari valmista skriptiä. Tarvitsee vain muuttaa muuttujia:
nymp - ymyröiden lkm
nsateet - moneenko osaan täyskulma 2*Pi jaetaan
rmax - suurimman ympyrän säde (aloitetaan origosta
(rmin - pienimmän ympyrän säde (kts. alempana))
f:=z-> ... - funktio, jolla kuvataan
Tässä esimerkkiajoja
> nymp:=7: nsateet:=15: rmax:=4: hr:=rmax/nymp: hTheta:=2*Pi/nsateet:
> display([seq(gsade(k*hTheta,0,rmax),k=0..nsateet),seq(gympyra(0,k*hr),
> k=0..nymp)]);
> f:=z->z+sqrt(z):
> display([seq(gsadekuva(k*hTheta,0,rmax,f),k=0..nsateet),seq(gympyrakuv
> a(0,k*hr,f),k=0..nymp)]);
Joukowski airfoil
Otetaan joukko 1+I-keskisiä ympyröitä, joiden säde vaihtelee, kuvataan
vain ympyräparvi, eikä säteitä. Animaatio on verraton, kun ympyräparvea
ei animoida, vaan pelkkästään kuvaparvea, silloin on helppo nähdä vastaavuudet.
> nymp:=20: rmin:=0.5:rmax:=2: hr:=(rmax-rmin)/nymp:
> display([seq(gympyra(1+I,rmin+k*hr),k=0..nymp)]);
> display([seq(gympyrakuva(1+I,rmin+k*hr,z->z+1/z),k=1..nymp)],insequenc
> e=true);
Tietysti voitaisiin vierittää parvea laittamalla keskipiste muuttumaan
jne. jne.
Iteraatiot
# Työarkista Lv7iter.mws
# Kiintopisteiteraatio
# Cobweb: (Robert Israelin tapaan)
porras:=x->plot([[x,x],[x,g(x)],[g(x),g(x)]]):
newtonkulma:=x->plot([[x,0],[x,f(x)],[N(x),0]],color=black);
N:=x->x-f(x)/D(f)(x);
Cobweb-piirroksia funktion "porras" avulla
Esim. 1
g:=x->(1/3)*(x^2+1);
X[0]:=1.0:n:=10:for i to n do X[i]:=g(X[i-1]) od:
gjaxkuva:=plot([x,g(x)],x=0 .. 1,color=[blue,black]):
alkup:=plot([[X[0],X[0]]],style=point,symbol=circle,symbolsize=20,color=blue):
display([alkup,gjaxkuva,seq(porras(X[i]),i=0..n)]);
KRE EXA 1 s. 926-927 (mikä painos?)
g1:=x->(1/3)*(x^2+1);
plot([x,g1(x)],x=0..5);
X:='X':f:='f':
f:=g1:
X[0]:=1.0:
n:=10:for i to n do X[i]:=g1(X[i-1]) od:
fjaxkuva:=plot([x,f(x)],x=0 .. 1,color=[blue,black]):
display({fjaxkuva,seq(porras(X[i]),i=0..n)});
Tässä vaiheessa voidaan muuttaa alkupistettä ja viedä kursori takaisin
(tai leikata/liimata) ja
laskea ja piirtää uudestaan. Muutetaan samalla x-skaala sopivaksi.
X[0]:=2.0;
n:=10:for i to n do X[i]:=f(X[i-1]) od:
fjaxkuva:=plot([x,f(x)],x=0 .. 4,color=[blue,black]):
display({fjaxkuva,seq(porras(X[i]),i=0..n)});
X[0]:=3.0:
n:=3:for i to n do X[i]:=f(X[i-1]) od:
fjaxkuva:=plot([x,f(x)],x=0 .. 4,color=[blue,black]):
display({fjaxkuva,seq(porras(X[i]),i=0..n)});
Suurempi n-arvo räjäyttää kuvan jo liian isoksi.
Vaihdetaan iterointifunktiota (edelleen saman funktion nollakohdista
on kyse.)
g2:=x->3-1/x;f:=g2:
X[0]:=1.0:
n:=10:for i to n do X[i]:=f(X[i-1]) od:
fjaxkuva:=plot([x,f(x)],x=0 .. 5,color=[blue,black]):
display({fjaxkuva,seq(porras(X[i]),i=0..n)});
X[0]:=2.0;n:=10:for i to n do X[i]:=f(X[i-1]) od:
fjaxkuva:=plot([x,f(x)],x=0 .. 4,color=[blue,black]):
display({fjaxkuva,seq(porras(X[i]),i=0..n)});
X[0]:=3.0:n:=10:for i to n do X[i]:=f(X[i-1]) od:
fjaxkuva:=plot([x,f(x)],x=2.5 .. 4,color=[blue,black]):
display({fjaxkuva,seq(porras(X[i]),i=0..n)});
Newtonin iteraatio ja fraktaalit Maplella
kts: V2-mappi: Projektit
# Kts. myös Kofler complexplot3d
newt:=proc(z)
local zz,m;
zz:=evalf(z);
for m from 0 to 50 while abs(zz^3-1)>=0.001 do
zz:=zz-(zz^3-1)/(3*zz^2);
od;
m
end:
newt3:=proc(z)
local zz,m,j1,j2,j3;
j1:=1.;
j2:=evalf(exp(I*2*Pi/3));
j3:=evalf(exp(-I*2*Pi/3));
zz:=evalf(z);
for m from 0 to 50 while abs(zz^3-1)>=0.001 do
zz:=zz-(zz^3-1)/(3*zz^2);
od;
if m > 45 then RETURN(10) fi;
if abs(zz-j1) < 0.1 then RETURN(0)
elif abs(zz-j2) < 0.1 then RETURN(1)
elif abs(zz-j3) < 0.1 then RETURN(2)
else RETURN(-1)
fi;
end:
# Tää on ihan nätti ja havainnollinen:
#plot3d('newt3'(x+I*y),x=-2..2,y=-2..2,grid=[20,20],
# view=[-2..2,-2..2,0..3], axes=FRAME, style=PATCH,orientation=[-130,50]);
# Valitse Color-valikosta Z HUE
# Tämä ei ole hyvä, värien hallinta?
#densityplot('newt3'(x+I*y), x=-2..2, y=-2..2,colorstyle=HUE,
# grid=[40,40], style=PATCHNOGRID, scaling=constrained, axes=box);
# Newtonin iteraatio ja complexplot3d:
#f := z-> z - (z^3-2)/(3*z^2);
#complexplot3d(f@@4,-3-3*I..3+3*I,view=-4..4,grid=[50,50],style=patch);
# Mandelbrot, Julia, ym, iterointifkt.
# Ande s. 386 ->
#
Mandelbrot:=proc(c)
local z,n,r,nmax;
z:=0.0;n:=0;r:=10000; nmax:=25;
while abs(z) < r and n < nmax
do
z:=z^2+c;
n:=n+1;
od;
return(n);
end:
# plot3d('Mandelbrot'(x+I*y),x=-2..2,y=-2..2,grid=[20,20],
# view=[-2..2,-2..2,0..25], axes=FRAME, style=PATCH,orientation=[-130,50]);
#densityplot('Mandelbrot'(x+I*y), x=-1.5..1.5, y=-1.5..1.5,colorstyle=HUE,
#grid=[40,40], style=PATCHNOGRID, scaling=constrained, axes=box);
Interpolaatio
Kirjoitetaan Lagrangen polynomin koodi.
Periaate:
Kirjoitetaan lauseke (x-x0)*(x-x1)...(x-xn) ja jaetaan se
(x-xj):llä.
Tämä antaa osoittajan.
Nimittäjä saadaan tästä korvaamalla x arvolla
xj. Näin johdutaan koodiin:
L:=proc(j,xd,x)
local oso,nimi,i,j1;
j1:=j+1;
oso:=product((x-xd[i]),i=1..nops(xd))/(x-xd[j1]);
nimi:=subs(x=xd[j1],oso);
oso/nimi;
end:
Parametrit
j -- antaa L:n indeksin, eli kyseessä on Lj
j=0..n, missä n on xd-listan pituus - 1
xd -- xdata, [x0,x1,...,xn], lista (n+1 pistettä)
x -- polynomin muuttuja
Lagrangen polynomijono muodostetaan nyt näin, kun data on listassa
xd
seq(L(j,xd,x),j=0..nops(xd)-1);
Esim
> L(0,[a,b,c],x);
(x - b) (x - c)
---------------
(a - b) (a - c)
# Yleisemmin:
> xd:=[seq(X[j],j=0..3)];
xd := [X[0], X[1], X[2], X[3]]
> seq(L(j,xd,x),j=0..nops(xd)-1);
(x - X[1]) (x - X[2]) (x - X[3])
-----------------------------------------,
(X[0] - X[1]) (X[0] - X[2]) (X[0] - X[3])
(x - X[0]) (x - X[2]) (x - X[3])
-----------------------------------------,
(X[1] - X[0]) (X[1] - X[2]) (X[1] - X[3])
(x - X[0]) (x - X[1]) (x - X[3])
-----------------------------------------,
(X[2] - X[0]) (X[2] - X[1]) (X[2] - X[3])
(x - X[0]) (x - X[1]) (x - X[2])
-----------------------------------------
(X[3] - X[0]) (X[3] - X[1]) (X[3] - X[2])
Interpolaatiopolynomi
> add(y[j]*L(j,xd,x),j=0..nops(xd)-1);
y[0] (x - X[1]) (x - X[2]) (x - X[3])
-----------------------------------------
(X[0] - X[1]) (X[0] - X[2]) (X[0] - X[3])
y[1] (x - X[0]) (x - X[2]) (x - X[3])
+ -----------------------------------------
(X[1] - X[0]) (X[1] - X[2]) (X[1] - X[3])
y[2] (x - X[0]) (x - X[1]) (x - X[3])
+ -----------------------------------------
(X[2] - X[0]) (X[2] - X[1]) (X[2] - X[3])
y[3] (x - X[0]) (x - X[1]) (x - X[2])
+ -----------------------------------------
(X[3] - X[0]) (X[3] - X[1]) (X[3] - X[2])
Määrittelemme funktioksi
lagint:=(xd,yd,x)->add(yd[j+1]*L(j,xd,x),j=0..nops(xd)-1):
Kokeillaan:
> lagint([a,b,c],[ya,yb,yc],x);
ya (x - b) (x - c) yb (x - a) (x - c) yc (x - a) (x - b)
------------------ + ------------------ + ------------------
(a - b) (a - c) (b - a) (b - c) (c - a) (c - b)
> p:=unapply(%,x);
p := x ->
ya (x - b) (x - c) yb (x - a) (x - c) yc (x - a) (x - b)
------------------ + ------------------ + ------------------
(a - b) (a - c) (b - a) (b - c) (c - a) (c - b)
> p(a);
ya
> p(b);
yb
> p(c);
yc
Hyvin toimii!
Interpolaatiovirhe:
interR:=(f,xd,x)->(D@@nops(xd))(f)(xi)/(nops(xd))!*product((x-xd[i]),
i=1..nops(xd)):
Kokeillaan:
interR(f,[a,b,c],x);
(3)
1/6 (D )(f)(xi) (x - a) (x - b) (x - c)
Kaunista!