Hatodik labor

Contents

Gyökkereső eljárások

Intervallumfelezés

function gyok=felezo(f,a,b,lepes)
% Függvény használata: felezo(f,a,b,lepes) ahol f egy egyváltozós függvény, 
% amelyre az [a,b] intervallumon teljesül, hogy f(a)f(b)<0. 
% Meghívásakor a függvény lepes darab iterációt végez.
if f(a)*f(b)>0
    error('Nem megfelelő kezdőintervallum')
elseif f(a)==0
    gyok=a;
    return
elseif f(b)==0
    gyok=b;
    return
end

c=a+(b-a)/2;
for i=1:lepes
    if f(c)==0
        gyok=c;
        return
    elseif f(c)*f(a)<0
        b=c;
    else
        a=c;
    end
    c=a+(b-a)/2;
end
gyok=c;

Példa az eljárás használatára:

gyok1=felezo(@sin,3,4, 10)
gyok1 =

    3.1411

pi-gyok1
ans =

   4.7937e-04

Newton-módszer

function gyok=newtonIt(f,df,x0,lepes)
% Függvény használata: newtonIt(f,df,x0) ahol f egy egyváltozós függvény, 
% df pedig a deriváltja.
% Meghívásakor a függvény lepes darab iterációt végez.
gyok=x0;
for i=1:lepes
    gyok=gyok-f(gyok)/df(gyok);
end
    

Példa az eljárás használatára:

gyok2=newtonIt(@sin,@cos,3,10)
gyok2 =

    3.1416

pi-gyok2
ans =

     0

Órai feladatok

Vizsgáljuk meg a konvergenciasebességet a két függvénynél. Próbáljuk ki, hogy a fenti két példában tényleg lineáris, illetve négyzetes-e a konvergencia. Mi lesz C értéke az $e_{k+1}<Ce_k$ becslésben a futás alapján?

A hibabecsléshez:

for i=1:10
    gyok1=felezo(@sin,3,4, i);
    gyok2=newtonIt(@sin,@cos,3,i);
    fprintf('Az %d-edik lépésben a hiba intervallumfelezéssel: %d, Newton-módszerrel: %d \n',i, gyok1-pi, gyok2-pi)
end
Az 1-edik lépésben a hiba intervallumfelezéssel: 1.084073e-01, Newton-módszerrel: 9.538895e-04 
Az 2-edik lépésben a hiba intervallumfelezéssel: -1.659265e-02, Newton-módszerrel: -2.893161e-10 
Az 3-edik lépésben a hiba intervallumfelezéssel: 4.590735e-02, Newton-módszerrel: 0 
Az 4-edik lépésben a hiba intervallumfelezéssel: 1.465735e-02, Newton-módszerrel: 0 
Az 5-edik lépésben a hiba intervallumfelezéssel: -9.676536e-04, Newton-módszerrel: 0 
Az 6-edik lépésben a hiba intervallumfelezéssel: 6.844846e-03, Newton-módszerrel: 0 
Az 7-edik lépésben a hiba intervallumfelezéssel: 2.938596e-03, Newton-módszerrel: 0 
Az 8-edik lépésben a hiba intervallumfelezéssel: 9.854714e-04, Newton-módszerrel: 0 
Az 9-edik lépésben a hiba intervallumfelezéssel: 8.908910e-06, Newton-módszerrel: 0 
Az 10-edik lépésben a hiba intervallumfelezéssel: -4.793723e-04, Newton-módszerrel: 0 

Órai feladat

Próbáljuk ki az x^3+3x^2-4 függvény gyökeinek keresésére az x0=-3 pontból indítva a Newton-módszert. Mi történt?

Megoldás

for i=1:10
    gyok=newtonIt(@(x)x^3+3*x^2-4,@(x)3*x^2+6*x,-3,i);
    fprintf('Az %d-edik lépésben a közelítő gyök: %d \n',i,gyok)
end
Az 1-edik lépésben a közelítő gyök: -2.555556e+00 
Az 2-edik lépésben a közelítő gyök: -2.297907e+00 
Az 3-edik lépésben a közelítő gyök: -2.155390e+00 
Az 4-edik lépésben a közelítő gyök: -2.079562e+00 
Az 5-edik lépésben a közelítő gyök: -2.040288e+00 
Az 6-edik lépésben a közelítő gyök: -2.020277e+00 
Az 7-edik lépésben a közelítő gyök: -2.010172e+00 
Az 8-edik lépésben a közelítő gyök: -2.005095e+00 
Az 9-edik lépésben a közelítő gyök: -2.002550e+00 
Az 10-edik lépésben a közelítő gyök: -2.001275e+00 

Látjuk, hogy jóval lassabb a konvergencia négyzetesnél. Ennek oka, hogy a -2 többszörös gyök, ami például az alábbi ábrából is kiderül:

t=linspace(-3,2);
plot(t,t.^3+3*t.^2-4,t, zeros(size(t)))

Ennek megfelelően az 1 közelében tényleg négyzetes a konvergencia

for i=1:10
    gyok=newtonIt(@(x)x^3+3*x^2-4,@(x)3*x^2+6*x,0.5,i);
    fprintf('Az %d-edik lépésben a közelítő gyök: %d \n',i,gyok)
end
Az 1-edik lépésben a közelítő gyök: 1.333333e+00 
Az 2-edik lépésben a közelítő gyök: 1.055556e+00 
Az 3-edik lépésben a közelítő gyök: 1.001949e+00 
Az 4-edik lépésben a közelítő gyök: 1.000003e+00 
Az 5-edik lépésben a közelítő gyök: 1.000000e+00 
Az 6-edik lépésben a közelítő gyök: 1 
Az 7-edik lépésben a közelítő gyök: 1 
Az 8-edik lépésben a közelítő gyök: 1 
Az 9-edik lépésben a közelítő gyök: 1 
Az 10-edik lépésben a közelítő gyök: 1