Pienimm\303\244n neli\303\266summan sovitus V2 19.3. 2002

with(LinearAlgebra):alias(Tr=Transpose):with(plots): xd:=[-1.3,-0.1,0.2,1.3]; yd:=[0.103,1.099,0.808,1.897]; n:=nops(xd); plot([seq([xd[i],yd[i]],i=1..n)],style=point,symbol=circle,symbolsize=15); datakuva:=%: C:=<<1,1,1,1>|<op(xd)>>; M:=Tr(C).C; B:=Tr(C).Vector(yd); ab:=LinearSolve(M,B); suora:=ab[1]+ab[2]*x; display(plot(suora,x=-2..2,color=blue),datakuva);

Muistellaan aluksi interpolaatioteht\303\244v\303\244\303\244. Esim: xd:=[$1..10];yd:=[1.3,3.5,4.2,5.0,7.0,8.8,10.1,12.5,13.0,15.6]; plot([seq([xd[i],yd[i]],i=1..10)],style=point,symbol=cross,symbolsize=18); p:=interp(xd,yd,x); display(%%,plot(p,x=0..10,color=blue)); Interpolaatio pakottaa polynomin usein voimakkaaseen heilahteluun, jolloin "mets\303\244 katoaa n\303\244kyvist\303\244 puilta". **** Kts esim. Runge tai Census ***** Luovutaan vaatimuksesta, ett\303\244 polynomin tulisi kulkea datapisteiden kautta. Valitaan alempiasteinen polynomi, joka sopivassa mieless\303\244 approksimoi parhaiten annettua dataa. Kriteeriksi valitaan usein funktiomallin (t\303\244ss\303\244 polynomi) ja datapisteiss\303\244 laskettujen erotuksten neli\303\266summan minimointi. T\303\244t\303\244 kautta menetelm\303\244 liittyy funktioiden \303\244\303\244riarvoteemaan. Polynomi-PNS-teht\303\244v\303\244ss\303\244 on annettu datapisteit\303\244 ( LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUYvNiVRImlGJ0YyRjUvRjZRJ25vcm1hbEYnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRictSSNtb0dGJDYtUSIsRidGPS8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R0ZILyUqc3ltbWV0cmljR0ZILyUobGFyZ2VvcEdGSC8lLm1vdmFibGVsaW1pdHNHRkgvJSdhY2NlbnRHRkgvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GLDYlLUYvNiVRInlGJ0YyRjVGOEY/Rj0=) (paljon) ja teht\303\244v\303\244n\303\244 on sovittaa esim. polynomi LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2KEYrLUYjNiYtRiw2JVEicEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnL0Y6USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGRC8lKXN0cmV0Y2h5R0ZELyUqc3ltbWV0cmljR0ZELyUobGFyZ2VvcEdGRC8lLm1vdmFibGVsaW1pdHNHRkQvJSdhY2NlbnRHRkQvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZTLUkobWZlbmNlZEdGJDYkLUYjNiQtRiw2JVEieEYnRjZGOUZARkBGQC1GPTYtUSI9RidGQEZCRkVGR0ZJRktGTUZPL0ZSUSwwLjI3Nzc3NzhlbUYnL0ZVRlxvLUYjNistSSVtc3ViR0YkNiUtRiw2JVEiYkYnRjZGOS1GIzYkLUkjbW5HRiQ2JFEiMEYnRkBGQC8lL3N1YnNjcmlwdHNoaWZ0R0ZbcC1GPTYtUSIrRidGQEZCRkVGR0ZJRktGTUZPL0ZSUSwwLjIyMjIyMjJlbUYnL0ZVRmJwLUYjNiYtRmFvNiVGY28tRiM2JC1GaW82JFEiMUYnRkBGQEZccC1GPTYtUTEmSW52aXNpYmxlVGltZXM7RidGQEZCRkVGR0ZJRktGTUZPRlFGVEZlbkZARl5wLUYjNictRmFvNiVGY28tRiM2JC1GaW82JFEiMkYnRkBGQEZccEZdcS1GIzYkLUklbXN1cEdGJDYlRmVuRmZxLyUxc3VwZXJzY3JpcHRzaGlmdEdGW3BGQEYrRkBGXnAtRiM2Jy1GYW82JUZjby1GIzYkLUZpbzYkUSIzRidGQEZARlxwRl1xLUYjNiQtRlxyNiVGZW5GZnJGXnJGQEYrRkBGK0ZARitGQEYrRkA= annettuun dataan niin, ett\303\244 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2KC1JK211bmRlcm92ZXJHRiQ2Jy1JI21vR0YkNi1RJiZTdW07RicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y9LyUpc3RyZXRjaHlHUSV0cnVlRicvJSpzeW1tZXRyaWNHRj0vJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGPS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjE2NjY2NjdlbUYnLUYjNiYtRiw2JVEiaUYnLyUnaXRhbGljR0ZCL0Y5USdpdGFsaWNGJy1GNTYtUSI9RidGOEY7Rj4vRkFGPUZDL0ZGRj0vRkhGPUZJL0ZMUSwwLjI3Nzc3NzhlbUYnL0ZPRltvLUkjbW5HRiQ2JFEiMUYnRjhGOC1GLDYlUSJuRidGVkZYRkkvJSxhY2NlbnR1bmRlckdGPUYrLUknbXNwYWNlR0YkNiYvJSdoZWlnaHRHUSYwLjBleEYnLyUmd2lkdGhHUSQ1LjBGJy8lJmRlcHRoR0ZbcC8lKmxpbmVicmVha0dRJWF1dG9GJy1GIzYkLUklbXN1cEdGJDYlLUkobWZlbmNlZEdGJDYkLUYjNictSSVtc3ViR0YkNiUtRiw2JVEieUYnRlZGWC1GIzYkRlNGOC8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLUY1Ni1RKCZtaW51cztGJ0Y4RjtGPkZnbkZDRmhuRmluRkkvRkxRLDAuMjIyMjIyMmVtRicvRk9GXXItRiM2Ji1GLDYlUSJwRidGVkZYLUY1Ni1RMCZBcHBseUZ1bmN0aW9uO0YnRjhGO0Y+RmduRkNGaG5GaW5GSUZLL0ZPRk0tRmpwNiQtRiM2JC1GX3E2JS1GLDYlUSJ4RidGVkZYRmRxRmZxRjhGOEY4RitGOEY4LUZebzYkUSIyRidGOC8lMXN1cGVyc2NyaXB0c2hpZnRHRmhxRjhGK0Y4RitGOA== minimoituu. Yleens\303\244 datapisteiden lkm =n on paljon suurempi kuin polynomin asteluku m (t\303\244ss\303\244 3). Jos A on matriisi , jonka alkiot ovat LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1JJW1zdWJHRiQ2JS1GLDYlUSJhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiYtRiw2JVEiaUYnRjdGOi1JI21vR0YkNi1RIixGJy9GO1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjkvJSlzdHJldGNoeUdGSi8lKnN5bW1ldHJpY0dGSi8lKGxhcmdlb3BHRkovJS5tb3ZhYmxlbGltaXRzR0ZKLyUnYWNjZW50R0ZKLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRiw2JVEiakYnRjdGOkZGLyUvc3Vic2NyaXB0c2hpZnRHUSIwRictRkM2LVEiPUYnRkZGSC9GTEZKRk1GT0ZRRlNGVS9GWFEsMC4yNzc3Nzc4ZW1GJy9GZW5GYm8tSShtc3Vic3VwR0YkNictRiw2JVEieEYnRjdGOi1GIzYkRj9GRkZnbi8lMXN1cGVyc2NyaXB0c2hpZnRHRlxvRmpuRkZGK0ZG , i=1..n, j=0..3, saadaan b-kertoimien m\303\244\303\244r\303\244\303\244miseksi ylim\357\277\275\357\277\275r\357\277\275ytyv\357\277\275 yht\357\277\275l\357\277\275systeemi A b = y, miss\357\277\275 siis A on n x 3 - matriisi. PNS-ratkaisuun p\357\277\275\357\277\275st\357\277\275\357\277\275n kertomalla t\357\277\275m\357\277\275 systeemi A:n transpoosilla, jolloin ratkaistavaksi tulee systeemi (1) LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYzLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEiQUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIlRGJ0YyRjVGMkY1LyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RJyZzZG90O0YnL0Y2USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGSC8lKXN0cmV0Y2h5R0ZILyUqc3ltbWV0cmljR0ZILyUobGFyZ2VvcEdGSC8lLm1vdmFibGVsaW1pdHNHRkgvJSdhY2NlbnRHRkgvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZXRi5GQC1GLzYlUSJiRidGMkY1LUZBNi1RIn5GJ0ZERkZGSUZLRk1GT0ZRRlNGVUZYLUZBNi1RIj1GJ0ZERkZGSUZLRk1GT0ZRRlMvRlZRLDAuMjc3Nzc3OGVtRicvRllGXm9GZ25GK0ZALUYvNiVRInlGJ0YyRjVGZ25GZ25GZ25GZ24tSShtZmVuY2VkR0YkNiQtRiM2LUYrRmduLUYvNiVRI29uRidGMkY1RmduRi4tRkE2LVEiOkYnRkRGRkZJRktGTUZPRlFGU0Zdb0Zfby1GLzYlUSJuRidGMkY1RmduLUYvNiVRK3RyYW5zcG9vc2lGJ0YyRjUtRkE2LVEiLkYnRkRGRkZJRktGTUZPRlFGU0ZVRlhGREZERkQ= (T\303\244m\303\244 johdetaan joko osittaisderivoimalla tai ortogonaaliprojektioon perustuvalla geometrisluonteisella argumentilla.) Jos tuntemattomia kertoimia on yht\303\244 monta kuin datapisteit\303\244 (polynomin asteluku = n-1), niin systeemill\303\244 A b = y on yksik\303\244sitteinen ratkaisu, mik\303\244li x-pisteet ovat erillisi\303\244. T\303\244ll\303\266in on kyse interpolaatiopolynomista. Yll\357\277\275 oleva ratkaisu voidaan Maplella tehd\357\277\275 rakentamalla matriisi A ja ratkaisemalla systeemi (1) . Maplen linalg:ssa on my\357\277\275leastsqrs- funktio, jolle annetaan pelkk\357\277\275 A - matriisi ja datapisteet, joka siis ratkaisee ylim\357\277\275\357\277\275r\357\277\275ytyv\357\277\275n systeemin A b = y PNS-mieless\357\277\275 T\357\277\275m\357\277\275n lis\357\277\275ksi Maplen stats-pakkauksessa on viel\357\277\275 automaattisempi funktio fit[leastsquare], joka tekee kaiken puolestamme (siis ei tarvitse muodostaa edes kerroinmatriisia). Suosimme kahta ensin mainittua tapaa ensinn\357\277\275kin opin kannalta, mutta my\357\277\275s siksi, ett\357\277\275 viimemainittu on syntaksiltaan outo ilmestys. restart:with(LinearAlgebra):with(linalg):with(plots):#with(stats): alias(Tr=Transpose,Van=VandermondeMatrix):

Sovita 2. asteen polynomi annettuun dataan xd:=[8,10,12,16,20,30,40,60,100]; yd:=[0.88,1.22,1.64,2.72,3.96,7.66,11.96,21.56,43.16]; n:=nops(xd); C:=Van(xd,n,3); # tai # Matrix((i,j)->xd[i]^(j-1),n,3); M:=Tr(C).C; b:=Tr(C).Vector(yd); a:=LinearSolve(M,b); p:=add(a[i]*x^(i-1),i=1..3); datakuva:=plot([seq([xd[i],yd[i]],i=1..n)],style=point,symbol=cross,symbolsize=18): display(plot(p,x=7..100,color=blue),datakuva); Tarkistetaan linalg-pakkauksen leastsqrs-funktiolla: aa:=leastsqrs(C,yd);

n:=5: m:=n-2:x:=[$1..n]; C:=Van(x,n,m);cond(Tr(C).C): evalf(%); n:=6: m:=n-2:x:=[$1..n]:C:=Van(x,n,m):cond(Tr(C).C): evalf(%); n:=10: m:=n-2:x:=[$1..n]:C:=Van(x,n,m):cond(Tr(C).C): evalf(%); cond(M) mittaa sensitiivisyytt\357\277\275 virheille (datan virheille ja py\357\277\275ristysvirheille). Normaaliyht\357\277\275l\357\277\275ratkaisun varjopuoli on suuri h\357\277\275ri\357\277\275alttius, jos teht\357\277\275v\357\277\275n koko on v\357\277\275h\357\277\275nk\357\277\275\357\277\275n suuri. n:=6: m:=2:x:=[$1..n]:C:=Van(x,n,m):cond(Tr(C).C): evalf(%); n:=20: m:=2:x:=[$1..n]:C:=Van(x,n,m):cond(Tr(C).C): evalf(%); n:=20: m:=3:x:=[$1..n]:C:=Van(x,n,m):cond(Tr(C).C): evalf(%); Paremmat menetelm\357\277\275t perustuvat ortogonalisointiin tai ns. singulaariarvohajoitelmaan.

Muodosta PNS-suora annettuun pisteist\357\277\275\357\277\275 Plista:=[[2.8,30],[2.9,26],[3,33],[3.1,31],[3.2,33],[3.2,35],[3.2,37],[3.3,36],[3.4,33]]; Matrix(Plista); plot(Plista); gp:=plot(Plista,style=POINT,symbol=circle): xd:=map(lista->lista[1],Plista); n:=nops(xd):yd:=map(lista->lista[2],Plista); A:=<Vector([1$n])| Vector(xd)>; Jos haluttaisiin jatkaa sarakkeittain latomista (ja siis korkeampiasteisten polynomien sovitusta), voitaisiin edet\357\277\275 t\357\277\275h\357\277\275n tapaan: (Alempana on kuitenkin yleisp\357\277\275tev\357\277\275mpi\357\277\275 keinoja.) <Vector([1$n])| Vector(xd) | Vector(map(t->t^2,xd))>; Tehd\357\277\275\357\277\275n aluksi "teoriaa" mukaillen: kertoimet:=LinearSolve(Tr(A).A,Tr(A).Vector(yd)); Suoremmin saadaan n\357\277\275in: leastsqrs(A,Vector(yd)); suora:=innerprod(kertoimet,<1,x>); gs:=plot(suora,x=2.8..3.5,color=blue): display({gs,gp});

Sovita a) 2. asteen PNS-polynomi ja b) interpolaatiopolynomi P:=[[1.4,7400],[1.8,7500],[2.4,7600],[3,7500],[4,7200]]; xd:=map(lis->lis[1],P); yd:=map(lis->lis[2],P); Jos tuo n\357\277\275ytt\357\277\275\357\277\275 oudolta, voitaisiin tehd\357\277\275 funktiot poimieka:=x->x[1]; poimitoka:=x->x[2]; ja suorittaa xd:=map(poimieka,P);yd:=map(poimitoka,P); T\357\277\275ss\357\277\275 on k\357\277\275tev\357\277\275 ja yleisp\357\277\275tev\357\277\275 tapa muodostaa matriisi. C:=Matrix(5,3,(i,j)->xd[i]^(j-1)); Toisaalta yleissivistykseen kuuluu my\357\277\275s tiet\357\277\275\357\277\275, ett\357\277\275 kyseess\357\277\275 on ns. Vandermonden matriisi. Siihenkin on linalg:ssa valmis funktio. Matrix(vandermonde(xd)):%[1..5,1..3]; # Vanhat tietorakenteet muutettava uusiksi (Matrix), jotta k\357\277\275tev\357\277\275t poimimispiirteet yms. toimisivat. K\357\277\275ytet\357\277\275\357\277\275n nyt sit\357\277\275 keskimm\357\277\275ist\357\277\275 "harmaata laatikkoa" kertoimet:=leastsqrs(C,yd); poly2:=innerprod(kertoimet,[1,x,x^2]); display(plot(P,style=point,symbol=circle,color=black),plot(poly2,x=1.3..4.2)); b) Kyseess\357\277\275 on interpolaatioply, jonka asteluku = datapisteiden lkm -1 = 4. No siihen l\357\277\275ytyy taas ihan valmiskin funktio: (Jos ratkaistaisiin LSQ-teht\357\277\275v\357\277\275n\357\277\275, matriisina olisi koko Vandermonde (transposilla kertominen voitaisiin "supistaa pois").) intpol:=interp(xd,yd,x); display(plot(intpol,x=1.3..4.2,color=blue),plot(P,style=point,symbol=circle,color=black),plot(poly2,x=1.3..4.2,color=gray));

LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEocmVzdGFydEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIjpGJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTA== QyYtSSV3aXRoRzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliRzYiNiNJLkxpbmVhckFsZ2VicmFHRiUhIiItSSZhbGlhc0dGJjYjL0kjVHJHRihJKlRyYW5zcG9zZUdGKEYr JSFH xd:=[1940, 1950, 1960, 1970, 1980, 1990]; yd:=[132.165, 151.326, 179.323, 203.302, 226.542, 249.633]; ydv:=Vector(yd): Pisteet:=[seq([xd[i],yd[i]],i=1..nops(xd))]; m:=nops(xd); Tehd\357\277\275\357\277\275n tyhm\357\277\275ll\357\277\275 tavalla sik\357\277\275li, ett\357\277\275 ei skaalata, mutta \357\277\275lykk\357\277\275\357\277\275ll\357\277\275 sik\357\277\275li, ett\357\277\275 k\357\277\275ytet\357\277\275\357\277\275n SVD:t\357\277\275 A:=VandermondeMatrix(xd,m,3); U,S,Vt:=SingularValues(evalf(A),output=['U','S','Vt']); JSFH S; S[1]/S[3]; Hyvin h\357\277\275iri\357\277\275altis. L\357\277\275hes LRV sarakkeet A-matriisissa d:=Tr(U).ydv; z:=seq(d[j]/S[j],j=1..3);z:=Vector([z]); c:=V.z; neliosumma:=add(d[i]^2,i=4..m); q:=<1 | x | x^2>.c; Aseteteaan pienin singulaariarvo 0:ksi, t\357\277\275ll\357\277\275in h\357\277\275iri\357\277\275luku S[1]/S[2]; on jo paljon siedett\357\277\275v\357\277\275mpi. z; T\357\277\275m\357\277\275 merkitsee, ett\357\277\275 kolmas z-komponentti pannaan 0:ksi (se saadaan valita mielivaltaisesti, mutta t\357\277\275m\357\277\275 on tapana). z[3]:=0; z; QyQ+SSJWRzYiLUkjVHJHRiU2I0kjVnRHRiUiIiI= c:=V.z; p:=<1 | x | x^2>.c; Kertoimet muuttuivat ihan toisiksi, mutta funktiot (polynomit) ovat varsin tarkasti samat ko v\357\277\275lill\357\277\275 plots[display](plot(Pisteet,style=point,symbol=circle,color=black),plot([q,p],x=1900..2000,color=[blue,red])); plot(q-p,x=1900..2000); neliosumma:=add(d[i]^2,i=3..m); Neli\357\277\275summa kasvoi hieman, mutta h\357\277\275iri\357\277\275alttius v\357\277\275heni paljon.