Maple-ohjelmia ja ``skriptejä'' k. 2001

Viitteitä

Ohjelmien ja koodien käyttöönotto

  cat > oma.txt
  Maalaus hiiren vasemmalla
  Liimaus hiiren keskimmäisellä
  CTR-D   tiedoston sulkeminen

Sitten Maple-istunnossa:

 > read("oma.txt");
 > read("c:\\windows\\mymaple\\oma.txt");  (PC-WIN:ssä)  

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");

Yleiskäyttöisiä apufunktioita

Datan muokkausta, yleisiä apufunktioita:

v2l:=vek->convert(vek,list):    # HAM-assa painovirhe (vec)
l2v:=lis->convert(lis,vector):

vappend:=(v1,v2)->vector([op(convert(v1,list)),op(convert(v2,list))]):

m2l:=M->map(op,convert(M,listlist)):  # Vastaa Matlabin M(:), mutta tässä
                                      # riveittäin 

Diskretointi vrt. Matlab-tuttavamme linspace

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!