Epälineaariset yhtälöt ja optimointi

Päivitys: 26.01.2014 HA
[Edellinen] [Seuraava] [Alkusivu]
Viitteitä: Halutkaamme tutkia ja ratkaista funktion $f(x)= 9\sin(x)-x$ nollakohtia.

Piirretään ja valitaan hiirellä klikkaamalla alkuarvaus nollakohdalle.

Tiedostossa fzero_johdantoesim.m on koodi: 
  close all    % Suljetaan grafiikkaikkunat
  f=@(x) 9*sin(x)-x
  fplot(f,[0 10]); grid on
  title('valitse hiirellä piste nollakohdan läheltä')
  [x0,y] = ginput(1);
  hold on
  plot(x0,0,'o');
  x00=fzero(f,x0);
  plot(x00,f(x00),'*')
  shg      % show graphics
  gtext('Alkuarvaus')
  gtext('Nollakohta') 
  gtext('fzero:lla laskettuna')
Tässä istunnon ajo kuvineen (yllä olevan m-tiedoston valinnalla FILE->PUBLISH).

Yhtälön numeerisen ratkaisun perusteita

Moler: NCM, zeros and roots - dokumentti on hyvä lähde opiskella perusteita tarkemmin.

Bisektio - välin puolitus

Perusasia: Välillä [a,b] jatkuva funktio f, jolla on päätepisteissä erimerkkiset arvot, saa arvon 0 jossain välin pisteessä. Liukuluvuilla laskettaessa tarkkaa nollakohtaa ei yleensä saavuteta, tavoite on asetettava näin:
Määritä hyvin pieni väli [a,b] (parhaassa tapauksessa 2 peräkkäistä liukulukua) siten, että f(a)f(b) < 0.
Bisektiomenetelmää havainnollistava koodi, mukalma C.van Loanin ShowBisect.m-koodista.
bisektio-funktio. (Moleria mukaillen)

Tehtävä:

  1. Kopioi koodi omaan m-tiedostoosi ja kokeile kumpaakin joihinkin funktioihin, vaikkapa johdantoesimerkin funktioon ja seuraaviin: ...
  2. Harjoittele debuggerin käyttöä: Katso m-tiedoston debug-valikon valintoja.
    Aseta "breakpoint" vaikka while-lauseen eteen ja suorita "continue" (F5). Katso muuttujia x,X,a,b työtilassa ja/tai m-tiedostossa pitämällä osoitinta kevyesti muuttujannimen päällä. Tee erilaisia kokeiluja, niistä voi olla arvaamaton hyöty "tosi paikan tullen".

Suppenemisnopeus

Bisektio-menetelmä suppenee hitaasti, mutta varmasti jotain nollakohtaa kohti, kunhan alkuväli on valittu oikein (ja funktio on jatkuva). Suppenemisnopeus ei riipu lainkaan tarkasteltavasta funktiosta. Joka askeleella virhe puolittuu, joten jos pyritään tarkkuuteen tol ja alkuväli on [a,b], niin n on valittava niin, että

(b-a)/2n < tol, eli n > log2((b-a)/tol).
Jos välin pituus on 1 ja tol=eps, niin 
>> n=log2(1/eps)  % Välin pituus = 1
n =
    52
>> n=log2(10/eps) % Välin pituus = 10
n =
   55.3219        % (Välin pituus ei juuri näyttele osaa.)

Newtonin menetelmä

Edellä ohjelmointi-osassa esiintyi yhden rivin "nuolinäppäinNewton", newt1. Siellä myös annettiin ohjelmointitehtävä, itse asiassa sama koodi on Molerin kirjan s. 119. (Tosin laskuri k ei ole välttämätön.)
Kirjoitin ajotiedoston (skriptin) newtonajo, jossa kutsutaan for-silmukassa newt1-funktiota. Samalla havainnollistetaan prosessia graafisesti.

Tehtävä

  1. Kopioi newtonajo itsellesi, testaile ja muokkaa oman makusi mukaan.
  2. Laske samoja laskuja myös omanewton-funktiolla (tehokkaampaa ja numeerisesti tarkempaa).

Eräs ongelmatapaus

Moler: NCM, 4.3 Zeros and roots s. 121 "A Perverse Example"

Ehtona kehälle on: $x_{n+1} - a = - (x_n - a).$
Se pätee, jos $$ \frac{f'(x)}{f(x)} = \frac{1}{2(x-a)},$$ josta integroimalla: $$f(x) = \text{sgn}(x-a)\sqrt{|x-a|}$$

Tehtävä Tutustu, aja ja muokkaa tarvittaessa tätä: esim/newtonpaha.m
Yleiset suppenemisehdot eivät päde, koska derivaatta menee äärettömyyteen O:ssa.

Newtonin ongelmat

Tehtävä:

Muokkaa Newtonin menetelmän ohjelmaasi Käyttämällä "symbolic toolboxia" derivaatan symboliseen muodostamiseen.
Ei ihan helppo. Kannattaa ensialkuun siirtää Maple-puolelle, tai tehdä Matlabin symbolipuolella (Mupad). Ei siis ota onnistuakseen tyyliin
syms u; f=@sin; df=@(u)diff(f(u),u)

fzero Yhdistelmä: sekanttimenetelmä, käänteinen kvadraattinen interpolaatio, bisektio (Perustuu: Richard Brent:n algoritmiin "zeroin" 1960-luvulla, Molerin viite [12])

Lue lisää:

  • Moler: NCM, zeros and roots

    Tässä ovat Molerin kirjan m-tiedostot Hae fzerogui.m Siinä on hieno havainnollistus fzero:n toiminnasta. Suorita: 

         fzerogui(@(x)besselj(0,x),[0 3.83])
    Seuraa vaiheita Molerin kirjasta.
    Lisätehtäviä s. 135 alkaen.

    Yksi sopiva: Hae 10 ensimmäistä nollakohtaa yhtälölle tan x = x. (Vrt. tangentin piirtämisteht.) 

    

    

    

    

    

    

    

    

    

    



    
[Edellinen]
[Seuraava]
[Alkusivu]