http://math.aalto.fi/~apiola/matlab/opas/lyhyt/polynomit.html Päivitetty 5.2.2019 HA

Polynomit, interpolaatio, käyrän sovitus

Polynomeista

Viitteitä:

Mathworks-doku, selkeä valaisevin käyttöesimerkein
Solmu-kirjoitus polynomeista ja interpolaatiosta, perusteellinen, tarkoitettu lukiolaisille.
"Moler: Numerical Computing with MATLAB, Interpolation"
Numeerisen analyysin oppikirjat, esim. Burden-Faires tai Cheney-Kincaid, ym.
Kahaner-Moler-Nash: Numerical methods and software.
Forsythe-Malcolm-Moler: Computer methods for mathematical computation
Van Loan: Intro to scientific computing ... with Matlab

Polynomin arvot, polyval

Tarkastellaan polynomia $$ p(x)=2 x^4 - x^2 + 5\, x +17$$
Polynomin määrää kerroinvektori [2,0,-1,5,17] .

Matlabissa polynomi esitetään kertoimien vektorina korkeimman asteen kertoimesta alkaen.
Jos x on jokin vektori, vaikkapa välin [-1,1] 10 jakopistettä jaettaessa 9:ään osaväliin, niin polynomin arvo vektorin x koordinaattipisteissä voitaisiin laskea tässä tapauksessa näin:

 >> x=linspace(-1,1,10);
 >> y=2*x.^4-x.^2+5*x+17;

Tämä on hankalaa ja tehotonta. Polynomin arvon laskemiseen on Malabissa funktio polyval.
Esim:

>> c=[2 0 -1 5 17];

>> polyval(c,0)
ans =
    17
>> y=polyval(c,-1)
y =
    13
>> y=polyval(c,1)
y =
    23
>> x=linspace(-1,1,10);
>> y=polyval(c,x)
y =
13.0000   13.2381   14.1041   15.2469   16.4324   17.5435   18.5802   19.6597   21.0159   23.0000

Yleinen n-asteinen polynomi $$p(x)=c_1 x^n + c_2 x^{n-1} + \ldots + c_n x + c_{n+1}$$ esitetään kerroinvektorina c. Siis korkeimman asteisesta kertoimesta alkaen. Huomaa, että Matlabissa vektorin indeksi alkaa aina 1:stä, ts. indeksointiyritys c(0) johtaa virheeseen. Siksi joudutaan hiukan tinkimään eleganssista, jonka indeksi 0 tarjoaisi.

Funktio polyval ei laske alussa esitetyllä "tyhmällä" tavalla, vaan aritmeettisia operaatioita säästävällä Hornerin menetelmällä, eli tyyliin:

p(x)=(... (c₁ x+ c₂)x + ...+ c_n)x + c_n+1

Tehtävä 1) Piirrä yllä tarkastellun (tai jonkin muun) polynomin kuvaaja vaikkapa välillä [-5,5] (tai jollain muulla välillä).
Neuvo: Polynomin arvot kannattaa laskea vektorissa x, jonka pituus on luokkaa 100. Aseta myös plot-komennon jälkeen grid on .

Polynomin juuret, roots, poly

Polynomin nollakohtien laskeminen on matematiikan historiassa saanut runsaasti huomiota osakseen. Käytännön laskennan kannalta siitä on jäänyt käteen vain toisen asteen yhtälön ratkaisukaava (hiukan kärjistäen). Kolmannen ja neljännen asteen yhtälöiden ratkaisukaavat ovat toivottoman mutkikkaita, ja korkeammille ei ratkaisukaavoja edes voida johtaa. Joka tapauksessa polynomiyhtälöiden ratkaisukaavojen kehittely (ja epäolemassaolo) on kiehtova kappale matematiikan historiaa. Solmu-lehdessä asiaa valaisee hauskalla tavalla Eero Saksman
kirjoituksellaan 3. asteen yhtälön ratkaisukaavan vaiheista

Numeeriikan käytännön kannalta polynomiyhtälö ei ole enempää eikä vähempää kuin eräs epälineaarinen yhtälö kaltaistensa joukossa.
Yleisiä epälineaarisia yhtälöitä käsitellään myöhemmin. Polynomiyhtälöt ovat tarkemmin tarkasteltuna kuitenkin siinä määrin erikoisia, että niille on voitu kehittää luotettavia menetelmiä, joilla saadaan approksimaatiot kaikille (kompleksi)juurille.

Matlab turvautuu vahvaan lineaarialgebrataustaansa ja hoitaa asian "mutkan kautta" muodostamalla ensin polynomin "kumppanimatriisin", "companion matrix", eli matriisin A, jonka karakterisitinen polynomi meillä on käsissämme. Sitten se laskee (numeerisesti) sen ominaisvektorit ja ominaisarvot. Huomaa, että numeerisesti edetään päinvastoin kuin lineaarialgebran harjoitus- ja tenttitehtävissä. Parempia ja tehokkaampiakin menetelmiä on kehitetty, mutta Matlabissa on tyydytty tähän, joka saadaan tärkeästä ominaisarvolaskennasta helppona (ja ovelana) sivutuotteena.

Aloitetaan 2. asteen polynomista, vaikkapa: p(x)=x² -x -1 .

>> c=[1 -1 -1];
>> z=roots(c)
z =
   -0.6180
    1.6180
>> poly(z)
ans =
    1.0000   -1.0000   -1.0000

Funktio poly on roots-funktion käänteisfunktio, se palauttaa polynomin kertoimet, kun juuret tunnetaan (ts. kertoo auki tulot (z-z₁)...(x-z_n) ja poimii kertoimet.).

Tehtävä: Laske kaikki ykkösen n:nnet juuret muutamalla n:n arvolla.
Piirrä yksikköympyrä ja nuo juuret rinkulalla (tms.) kompleksitasoon.
Ohje Tarkastele polynomiyhtälöä zⁿ - 1 = 0.

Funktiot polyder,polyint, conv,deconv

Polynomin derivointi- ja integrointi ovat sopivia pieniä Matlab-ohjelmointiharjoituksia. Saat miettiä, mutta ehkä annetaan varsinaisesti tehtäväksi ohjelmoinnin yhteydessä. (Ratkaisut ovat valmiina polyder ja polyint-funktioina, tosin erikoistapauskäsittelyt vaikeuttavat hiukan koodin lukemista.)
Polynomien kertominen on hiukan "haastavampaa", mutta aivan tehtävissä ja varsin opettavainenkin. Ei nyt kuitenkaan upoteta siihen aikaresurssiamme. Sekä kertominen että jakaminen hoituvat conv ja deconv-funktioilla, jotka ovat keskeisiä työkaluja signaalinkäsittelyssä.

Polynomin sovitus dataan

Olkoon annettu datapisteet:

      x₁ x₂ ...  x_m 
      y₁ y₂ ...  y_m

Yleinen polynomin sovitustehtävä on valita korkeintaan astetta $m-1$ oleva polynomi $p$, jolle $$p(x_k) \approx y_k, k=1, \ldots, m$$ Huom: Yleensä kirjallisuudessa datan indeksointi alkaa nollasta, jolloin valitaan "korkeintaan astetta m" oleva polynomi. Noudatetaan tässä Matlabin indeksointityylissä, joka ei salli indeksiä 0.

Jos vaaditaan, että yllä oleva ehto toteutuu tarkasti, on kyseessä interpolaatiotehtävä. Jos x_k-pisteet ovat erilliset, tällä on yksikäsitteinen ratkaisu. Interpolaatioehdot voidaan kirjoittaa yhtälöryhmäksi, jonka matriisiksi tulee annettuun x-dataan liittyvä Vandermonden matriisi. Kts. vander. (Matlab:ssa: help vander, doc vander).

Tehtävä: Kirjoita funktio

    function [kertoimet,condnr]=vandinterp(xdata,ydata)
    % Funktio laskee interpolaatiopolynomin kertoimet Vandermonden
    % systeemin ratkaisulla ja palauttaa myös cond-luvun.
    %  Esim:
    %  xdata=0:5; ydata=xdata.*sin(xdata);
    %  [c,cnr]=vandinterp(xdata,ydata);

Laske vaikkapa kommenttiesimerkin tapaus ja vertaa polyfit-funktion antamiin kertoimiin. Piirrä data ja interpolaatiopolynomi. Käytä arvojen laskentaan polyval-funktiota.

Vandermonden matriisi on ei-sigulaarinen, mutta sen "häiriöalttius" kasvaa nopeasti m:n kasvaessa. Interpolaatiotehtävälle voidaan johtaa tyylikkäämpiä ja numeerisesti vähemmän häiriöalttiita ratkaisutapoja. Menetelmien isät ovat herrat Lagrange ja Newton.
Kts. HA:n luentomatskua interpolaatiosta, ja käyrän sovituksesta (englanniksi)

Molerin vapaasti verkosta saatavassa kirjassa Numerical Computing with MATLAB on luvut "Interpolation" ja "Least Squares" Edellisessä on esitetty Lagrangen menetelmä ja sen toteuttaminen vähäeleisen yksinkertaiseksi Matlab-koodiksi nimeltään polyinterp. Se ei kuitenkaan sisälly Matlabiin, vaan palvelee opetustarkoitusta. Tähän voidaan palata ohjelmoinnin yhteydessä. Kannattaa huomata, että polyinterp toimii hyvin myös symboliselle muuttujalle (kun Symbolic toolbox on käytössä).

`polyfit`

Funktio polyfit määrittää interpolaatiopolynomin kertoimet. Se ratkaisee lisäksi ns. PNS-apprioksimaatiotehtävän tapauksessa, jossa haetaan alemman asteista polynomia kuin m-1 (missä siis m=datapisteiden lukumäärä).

Molerin kirjan edellä jo mainitussa LSQ-luvussa mainittu "design matrix" on polynomiapproksimaation tapauksessa Vandermonden matriisi, jossa on vähemmän sarakkeita kuin rivejä. Tämän tehtävän polyfit laskee "takakenolla", johon on sisään rakennettuna PNS-ratkaisu, kuten edellä jo puhuttiin.
Itse asiassa tämä ei ole sen maagisempaa kuin ns. normaaliyhtälöiden muodostaminen ja ratkaisu. Jos V on tuo "kapea" Vandermonde, niin PNS-kertoimet c saadaan yhtälösysteemistä $$(V^T V) c = V^T y,$$ missä y on data-y-vektori.

Tehtävä: Tee pieni modifikaatio yllä olevaan vandinterp-funktioon muodostaaksesi funktion:

function kertoimet=vanderfit(xdata,ydata,n)
% Astetta n olevan PNS-polynomin sovitus
% vandinterp on tämän erikoistapaus.

Tässä ei tarvitse edes muodostaa (V^TV)-matriisia, sillä "takakeno" suorittaa ylimääräytyvän systeemin PNS-ratkaisun. Kyseessä on toiminnallisesti sama kuin polyfit.

Palapolynomit, splinit

Jos annetut datapisteet yhdistetään niiden kautta kulkevilla janoilla, kysessä on paloittain lineaarinen interpolaatio. MATLAB:n plot tekee sen. Kokeile vaikka:

>> x=0:5;
>> y=rand(1,6);
>> plot(x,y)        % Lineaarinen interpolaatio
>> ylim([-1 2])     % Katsotaan hiukan kauempaa.
>> hold on
>> plot(x,y,'o','MarkerSize', 10)  % Merkitään rinkuloilla pelkät datapisteet.

Kun halutaan interpoloida polynomeilla niin, että käyrä näyttää sileältä, mutta asteluku pysyy alhaisena, on suosituimmaksi muodostunut astetta 3 olevista polynomipaloista muodostuva funktio, joka liitoskohdissa interpoloi dataa (on siis jatkuva) ja jonka 2 ensimmäistä derivaattaa ovat jatkuvia liitoskohdissa. Tällä tavoin vältytään yllä esitellyistä polynomi-interpolaatio-ongelmista.

Jos vaikka interpoloidaan funktiosta (1/(x + (1-x)²)) välillä [-2,2] valittua 20 tasavälisen pisteen muodostamaa dataa, toimitaan näin:

xdata=linspace(-2,2,20);
ydata=1./(xdata+(1-xdata).^2);
x=linspace(-2,2,100);    % Tiheä pisteistö, jossa arvot lasketaan
y=spline(xdata,ydata,x);
plot(xdata,ydata,'o',x,y)
title('Splinisovitus')

Koko istunto kuvineen on tässä

Tehtävä: Tee yllä olevaan dataan 3:n asteen PNS-polynomisovitus ja piirrä samaan kuvaan vaikka eri värillä tai viivatyypillä. Käytä legend-komentoa selitteiden lisäämiseksi. Myös grid voi selventää.

Muita interpolointifunktioita

Palapolynomeja voidaan muodostaa monella muullakin tavalla, astelukua voidaan muutella, ehtoja voidaan asettaa derivaatoille, ns. reunaehtoja voidaan vaihdella. MATLAB:ssa on spline toolbox erilaisiin splinitarpeisiin. Katsotaan vielä yhtä, suhteellisen tuoretta perus-MATLAB:iin kuuluvaa: pchip-funktiota. ("Piecewise Hermite interpolation polynomial"). Tässä ei vaadita kolmannen asteen polynomipalojen toisten derivaattojen jatkuvuutta liitoskohdissa. Tuloksena on yleensä vähemmän sileä käyrä, jolla on se etu puolellaan, että se säilyttää datan muodon ja monotonisuuden. Ts. pchip tuottaa funktion, joka on monotoninen niillä väleillä, joilla data on monotoninen ja sillä on samoilla väleillä lokaali ääriarvo kuin datalla. Tavallinen kuutiollinen splini puolestaan pakotetaan sileytensä takia heilahtelemaan, vaikka data ei heilahtelisi.

Tämä istunto ja kuva puhukoot puolestaan

Tehtäviä

Hae Molerin m-tiedostoista interpgui ja censusgui. Kirjan luvussa 3 s. 108 ja 116 on edelliseen liittyvää selvitystä ja tehtäviä, luvussa 5 s. 144 ja 159 jälkimmäiseen.
Molerin m-tiedostoista Näissä luvuissa hyviä lisätehtäviä.

[Edellinen] [Seuraava] [Alkusivu]