Harjoitustehtäviin 2 liittyviä ohjeita ja ehdotuksia

Yleistä

Alkuviikko

Kaikissa tehtävissä tarvitaan tämäntyylistä hommaa. (Toisinaan 10 lukua kertoo enemmän kuin yksi kuva.)

g:=x->...;
x0:=0.0: X:=seq((evalf(g@@k)(x0)),k=0..10);

porras:=x->plot([[x,x],[x,g(x)],[g(x),g(x)]],color=black):
newtonkulma:=x->plot([[x,0],[x,f(x)],[N(x),0]],color=black);
N:=x->x-f(x)/D(f)(x);

Kiintopisteiteraatiokuva

gjax:=plot([x,g(x)],x=3..6):
display([gjax,seq(porras(X[k]),k=1..11)]);

Newtoniteraatiokuva

f:=x->...; # Ajatellaan, että N on määritelty yllä.
x0:=2.0: X:=seq(evalf((N@@k)(x0)),k=0..10);
display([plot(f(x),x=-2..2),seq(newtonkulma(X[k]),k=1..11)]);

Teht. 1

Kohdassa (c) etsitään pienintä positiivista, kuten (b):ssä.

Jääköön tässä vapaaehtoisuuden varaan kiintopistelauseen tarkka soveltaminen, mutta kuva on syytä piirtää ja perustella ainakin sen avulla suppeneminen.

Teht. 2

Teht. 3

Teht. 4

Tietenkin päättelyyn kuuluu myös kiintopistelauseen käyttö, muutenhan johtopäätös on vain, että jos suppenee, niin suppenee aivan hirvittävän lujaa.

Teht. 5

Kummassakin tapauksessa kannattaa hakea iterointifunktion minimi. Siitä oikealle toimii KP-lause. Vasen taas rävähtää ensi-iteraatiossa (hyvälle) oikealle (samaan tapaan, mutta vielä "kiltimmin") kuin tehtävässä 3.

Teht. 6

Loppuviikko

Maple-koodien kokoaminen ja dokumentointi

 read("v201.mpl"):
 Myös suoraan:
 read("/p/edu/mat-1.414/maple/v201.mpl"):

   cat > oma.mpl
   Liimas hiiren kesimmäisellä
   CTR-D   tiedoston loppu.
 

   read("oma.mpl");
   print(newt3);    # Katsotaan, miltä koodi näyttää.

Ohjelmankehitys

plot

 plot(x^2,x=-1..1);         # 1. argumentti on lauseke  (ei funktio)
 plot([x,x^2,x^3],x=-1..1); # useita kuvaajia samaan kuvaan.

  plot(f(x),x=-1..1);

  plot([x->x^2,x->x^3],-1..1);

Datapisteiden piirto

pisteet:=[[x1,y1],[x2,y2],[x3,y3],[x4,y4]];
plot(pisteet);  # jos style=point-määrettä ei anneta, yhdistetään janoilla.

Iteraatiot

N:=x->x-f(x)/D(f)(x);
f:=x-> ... ;

  fN:=evalf@N   # "floatingpointNewton"

 fN(2.1);fN(%);fN(%);fN(%);fN(%);  # Askel kerrallaan edelliseen tulokseen sovellettuna.

n:=5:
X1:=seq(fN@@k)(x0)),k=0..n);       # n iteraatiota kerrallaan
X2:=seq(fN@@k)(X1[n+1])),k=0..n);  # Jatketaan siitä, mihin edellä päästiin.
X3:=seq(fN@@k)(X2[n+1])),k=0..n);  #  -- " --
X:=X1,X2,X3;                       # Koko jono.

Iteraatio for-silmukalla

Siis näin:

  x[0]:=2.2: n:=5:
  for k from 1 to n do
     X[k]:=evalf(N(X[k-1])) 
  od;

                                X[1] := N(X[0])
                              X[2] := N(N(X[0]))
                             X[3] := N(N(N(X[0])))
                           X[4] := N(N(N(N(X[0]))))
                          X[5] := N(N(N(N(N(X[0]))))) 

Teht. 1

Teht. 2

Teht. 3

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:

# j1:=1.:j2:=evalf(exp(I*2*Pi/3)):j3:=evalf(exp(-I*2*Pi/3)):
z2xyw:=(z,w)->[Re(z),Im(z),w]: # "Kertakäyttöfkt."
# juuret3d:=spacecurve([z2xyw(j1,0),z2xyw(j2,1),z2xyw(j3,2)],style=point,symbol=circle,color=black,symbolsize=20):

#with(plots):
# Piirretään "alkuarvokuva", jossa kutakin juurta kohti suppeneva alkupiste
# nostetaan sille tasolle, joka on juuren numero. Väritetään korkeuden mukaan.
# Värikoodattu versio saadaan katsomalla suoraan ylhäältä
#aakuva:=plot3d('newt3'(x+I*y),x=-2..2,y=-2..2,grid=[40,40],
view=[-2..2,-2..2,0..3],axes=FRAME,style=PATCHNOGRID,orientation=[-130,50]):
# Valitse Color-valikosta Z HUE

Otetaan jossain vaiheessa Matlab-tyyli, silloin voidaan hienontaa.

Toinen mahdollisuus

> g := z-> z - (z^3-2)/(3*z^2);
> complexplot3d(g,-3-3*I..3+3*I,view=-4..4,grid=[20,20],style=patch);
# 2. iteraatiota:
> complexplot3d(g@@2,-3-3*I..3+3*I,view=-4..4,grid=[20,20],style=patch,a
> xes=BOX);

Teht. 4

Tuo kahden pituinen sykli voidaan hakea näin: Vaaditaan, että N(x)=-x, johtuu N:n parittomuudesta, mieti!

   solve(N(x)=-x,x);

Jos aikaa ja kokeilunhalua on, voisi kehitellä jotain tämäntyylistä.

newt2:=proc(x)
 local xx,j0,j1,j2,m,g;
 g := x->2*x^3/(3*x^2-1) ;
 j0:=0.; j1:=1/sqrt(5.); j2:=-1/sqrt(5.); 
 xx:=evalf(g(x));
 for m from 0 to 50 while abs(xx-g(xx))>=10^(-8)   do
     xx:=g(xx);
  od;
if m > 45 then RETURN(10) fi;
  if abs(xx-j0) < 0.1 then RETURN(0)
  elif abs(xx-j1) < 0.1 then RETURN(1)
  elif abs(xx-j2) < 0.1 then RETURN(2)
  else RETURN(-1)
  fi;
end:
seq(newt2(0.1*k),k=-10..10);
  -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -1
newt2(1/sqrt(5));
                                  10
plot('newt2',1/sqrt(5)..1/sqrt(3));

map(newt2,linspace(0,1/2,100));  # Hyvä tekniikka
                                 # linspace on oma "matlab-tyylinen" funktio.
map(newt2,linspace(1/sqrt(5),1/sqrt(3),100));  # 0:lla jako

Sinänsä 1-ulotteinen tilanne on hyvä "fraktaali-ideoiden" testaamiseen, koska laskenta on paljon kevyempää kuin tasossa.

Teht. 5

fsolve-komentoon kannattaa samalla tutustua (ellet jo ennen ...).