Numerické řešení rovnice f(x) = 0

(1)

Numerické řešení rovnice f (x) = 0

Přemysl Vihan 9.10.2003

Katedra fyziky, Pedagogická fakulta Univerzity J.E. Purkyně v Ústí n.L.

2. ročník, PMVT-mag.

Abstrakt

Seminární práce se zabývá numerickým řešením rovnice f(x) = 0. Jsou uvedeny metody, které vyžadují separaci (lokalizaci) kořenů:

metoda půlení intervalu a metoda regula falsi (sečen) a metody, které vyžadují ”dobrý” odhad počáteční aproximace: prostá iterační metoda a Newtonova metoda. Použití metod, jejich výhody a nevýhody jsou demonstrovány na několik příkladech.

1 Úvod

V praxi se často setkáváme s úlohou řešit rovnici

f(x) = 0 (1)

kde f(x) je reálná funkce proměnné x. Řešit rovnici (1) znamená nalézt všechna ξ, pro která platí

f(ξ) = 0 (2)

ξ nazýváme kořenem či řešením rovnice (1).

Jen málo rovnic typuf(x) = 0 lze řešit analyticky (např. algebraické rovnice do 3. řádu). Většinou však musíme použít numerické (přibližné) metody pro řešení rovnic f(x) = 0 [1].

(2)

Řešit rovnici (1) numericky znamená navrhnout algoritmus přibližného řešení, při kterém získáme posloupnost aproximací

x₀, x₁, ..., x_k, ... (3)

takovou, že

k→∞lim x_k =ξ (4)

2 Teorie

Numerické metody pro řešení rovnice f(x) = 0 lze rozdělit na metody, které vyžadují separaci (lokalizaci) kořenů a metody, které vyžadují ”dobrý” odhad počáteční aproximacex0. Separovat kořeny znamená nalézt interval< a, b >, pro který platí:

• f(a)f(b)<0,

• funkce f(x) je v intervalu spojitá a

• v intervalu leží jen jeden kořen ξ.

2.1 Metoda půlení intervalu

Tato metoda vyžaduje separaci kořenů, t.j. znalost intervalu < a₀, b₀ >, ve kterém se nachází hledaný kořen. Pro krajní body intervalu musí platit

f(a₀)f(b₀)<0 (5)

t.j. funkce v krajních bodech intervalu nabývá opačných znamének. V prvním kroku rozdělíme interval na dvě části - střední hodnotu označme x1:

x₁ = a₀+b₀

2 (6)

Poté testujeme, zda x1 neodpovídá hledanému řešení:

|f(x₁)|< p (7)

(3)

kde pje požadovaná přesnost. Není-li tato podmínka splněna, testujeme zda

f(a₀)f(x₁)<0 (8)

V takovém případě se kořen nachází v intervalu < a₀, x₁ >. V opačném případě kořen leží v intervalu< x1, b0 >. Tento postup opakujeme tak dlouho, dokud není řešení nalezeno s požadovanou přesností.

2.2 Metoda regula falsi (sečen)

Tato metoda též patří mezi metody vyžadující separaci kořenů. Je podobná metodě půlení intervalu - opět vymezíme interval < a, b >, ale bod x získá- váme jako průsečík sečny spojující oba krajní body s osou x. Rovnice této sečny je

y= f(b)−f(a)

b−a (x−a) +f(a) (9)

Aproximace řešení x_k se získá z rovnice (9), položíme-li y= 0:

x_k= a_k−1f(b_k−1)−b_k−1f(a_k−1)

f(b_k−1)−f(a_k−1) (10)

Metoda sečen obvykle konverguje rychleji, než metoda půlení intervalu.

2.3 Prostá iterační metoda

Prostá iterační metoda nevyžaduje separaci kořene, na druhou stranu se ne- obejde bez „dobréhoÿ odhadu počáteční aproximace x₀. Funkci f(x) je dále nutné upravit na tvar

x=ϕ(x) (11)

V k-tém kroku iterace potom provedeme odhad kořene jako

xk=ϕ(xk−1) (12)

Prostá iterační metoda nemusí vždy konvergovat.

(4)

2.4 Newtonova metoda

Jde opět o metoduiterační vyžadující odhad počáteční aproximacex₀. Apro- ximace řešení xk+1 se získá jako průsečík tečny v bodě [xk, f(xk)] s osou x

xk+1 =xk− f(xk)

f⁰(x_k) (13)

Někdy může být problémem nutnost analytického výpočtu derivace. V tom případě lze použít upravenou verzi Newtonovy metody, která používá derivaci numerickou. Tečnu v boděx_k nahradíme přímkou procházející body x_k ax_k−1. Je-li x₀ počáteční aproximací, hodnotu x₁ určíme podle vztahu

x₁ =x₀(1 +²) (14)

kde ² je dostatečně malé číslo (např. 0.01).

V dalších krocích (k = 2,3, ...) provede odhad podle vztahu (13) (stejně jako u klasické Newtonovy metody) s tím, že hodnotu f⁰(x) nahradíme vý- razem

f⁰(x)≈ f(x_k)−f(x_k−1)

x_k−x_k−1 (15)

Ve srovnání s prostou iterační metodou konverguje Newtonova metoda rychleji. Nemusí konvergovat vždy, ale často funguje i v případě, kdy prostá iterační metoda selhává.

3 Výsledky a diskuze

Numericky jsme řešili rovnice:

8x³−6x−1 = 0 (16)

x−sinx−π

4 = 0 (17)

Výpočet jsme provedli pomocí výše zmíněných metod: metodou půlení intervalu, metodou sečen (regula falsi), prostou iterační metodou a Newtono- vou metodou. U Newtonovy metody jsme zkoumali dvě verze – s analytickým výpočtem derivace f⁰(x) a s jeho numerickým přiblížením podle vzorce (15).

Přesnost výpočtu byla nastavena na hodnotu ²= 10⁻¹⁰ a výpis programů v jazyce Fortran je v Příloze.

(5)

3.1 Metoda půlení intervalu

Metodou půlení intervalu byly nalezeny tyto kořeny:

Funkce f₁(x) = 8x³ −6x−1 = 0

poč. interval x f1(x) počet iterací

<−1,−0.5> −0.766044443123974 −4.03855113413620·10⁻¹¹ 32

<−0.5,0> −0.173648177675204 4.36562151446412·10⁻¹¹ 32

<0.5,1.5> 0.939692620784626 −1.94774608767989·10⁻¹¹ 36 Funkce f₂(x) = x−sinx− ^π₄ = 0

poč. interval x f₂(x) počet iterací

<0,3.14> 1.76634028665372 6.07372557738690·10⁻¹¹ 34 Shrnutí

Metoda půlení intervalu je velmi jednoduchá metoda s relativně pomalou konvergencí k hledanému kořeni – v našem testu byla metoda půlení intervalu nejpomalejší. Při správně provedené separaci kořene však vždy konverguje k přesnému řešení.

3.2 Metoda regula falsi (sečen)

Metodou regula falsi byly nalezeny tyto kořeny:

Funkce f1(x) = 8x³ −6x−1 = 0

poč. interval x f₁(x) počet iterací

<−1,−0.5> −0.766044443110925 6.50997606516102·10⁻¹¹ 25

<−0.5,0> −0.173648177668642 9.03027197289574·10⁻¹² 11

<0.5,1.5> 0.939692620779913 −9.10867358101808·10⁻¹¹ 37 Funkce f₂(x) = x−sinx− ^π₄ = 0

poč. interval x f₂(x) počet iterací

<0,3.14> 1.76634028652257 −9.58941248890149·10⁻¹¹ 21

(6)

Shrnutí

Metoda regula falsi konverguje ve srovnání s metodou půlení intervalu o poznání rychleji. Určitý problém však může nastat ve vzorci (10), objeví-li se ve jmenovateli příliš malé číslo a může dojít k eventuálnímu dělení nulou.

3.3 Prostá iterační metoda

Funkce f₁(x) = 8x³ −6x−1 = 0

Obě funkce jsme upravili podle vzorce (11). U funkce f1(x) jsme testovali dva tvary rovnice (11):

x = 8x³−1

6 (18)

x = 1

2

√3

6x+ 1 (19)

S použitím funkce (18) jsme došli k následujícímu řešení:

poč. odhad x f1(x) počet iterací

−0.25 −0.173648177676812 5.21391611526048·10⁻¹¹ 12 Použijeme-li funkci (19), metoda konvergovat nebude.

Funkce f2(x) = x−sinx− ^π₄ = 0 Funkci f₂(x) jsme upravili na tvar

x= sin(x) + π

4 (20)

Výsledek výpočtu:

poč. odhad x f₂(x) počet iterací

1.5 1.76634028665244 5.92004918394465·10⁻¹¹ 14

(7)

Shrnutí

Jde o jednoduchou metodu, která ne vždy bude konvergovat ke správnému řešení. V našem případě se u funkce f1(x), nepodařilo najít takový odhad, aby metoda konvergovala k jinému než k výše uvedenému kořeni.

3.4 Newtonova metoda

Při použití Newtonovy metody byly nalezeny tyto kořeny:

Funkce f₁(x) = 8x³ −6x−1 = 0

poč. odhad x f₁(x) počet iterací

−1 −0.766044443118984 −4.91854820761084·10⁻¹⁴ 5

−0.25 −0.173648177666930 5.32343266690383·10⁻¹⁷ 4 0.75 0.939692620785911 3.94601885889134·10⁻¹⁴ 5 Funkce f₂(x) = x−sinx− ^π₄ = 0

poč. odhad x f₂(x) počet iterací

1.5 1.76634028660287 8.96000938374608·10⁻¹⁵ 4

3.5 Newtonova metoda s numerickou derivací

Po úpravě vzorce (13) pomocí vztahu (15) jsme dostali následující výsledky:

Funkce f₁(x) = 8x³ −6x−1 = 0

poč. odhad x f₁(x) počet iterací

−1 −0.766044443118978 −3.41393580072236·10⁻¹⁵ 7

−0.25 −0.173648177666931 3.42152521592975·10⁻¹⁵ 5 0.75 0.939692620785909 6.66133814775094·10⁻¹⁶ 7 Funkce f₂(x) = x−sinx− ^π₄ = 0

poč. odhad x f2(x) počet iterací

1.5 1.76634028660001 −3.40559947863833·10⁻¹² 5

(8)

Shrnutí

Newtonova metoda opět nemusí vždy konvergovat. Pokud však konverguje, konverguje velmi rychle (šlo o nejrychlejší metodu v testu). Nevýhodou kla- sické Newtonovy metody je nutnost analyticky derivovat vstupní funkci. Tuto nevýhodu však lze eliminovat náhradou derivacef⁰(x) jejím numerickým při- blížením – tato verze Newtonovy metody potřebuje k dosažení výsledku jen o několik (v našem případě o 1 až 2) iterace více, což v praxi nebude prav- děpodobně hrát významnou roli.

4 Závěr

S využitím výše popsaných metod se podařilo najít všechny kořeny studo- vaných rovnic. Jak ilustrují výše uvedené tabulky, úspěch při numerickém řešení rovnicef(x) = 0 závisí na dobré separaci kořenů či na dobrém odhadu počáteční iterace. Rozdíl v počtu provedených kroků se v praxi (vzhledem k výkonu počítače) neprojevil. Nejrychleji byl však výpočet proveden pomocí Newtonovy metody. Tato metoda se také ukázala – v našem konkrétním pří- padě – jako nejspolehlivější: s její pomocí jsme s požadovanou přesností určili všechny kořeny studovaných rovnic.

Reference

[1] M. Vicher, Numerická matematika(2003).

(9)

Příloha: Programy v jazyce Fortran

Separační metody (metoda půlení intervalu a metoda sečen)

subroutine separace(funkce,a0,b0,metoda,vysl,err) external funkce

real(8) :: funkce

real(8), intent(in) :: a0, b0 ! poc. interval integer, intent(in):: metoda ! pouzita metoda,

! 0=puleni, 1=secny real(8), intent(out) :: vysl ! nalez. koren integer, intent(out) :: err ! 0=vypocet OK

! 1=prekrocen max.

! pocet it.

! 2=deleni nulou real(8) :: a,b,x,fx,fa,fb

integer :: citac_it integer :: citac_tisk

real(8),parameter::PRESNOST = 1e-10 integer,parameter::MAX_POCET_IT = 10000

! max. pocet iter.

integer,parameter::TISK = 1

! --- inicializace --- a = a0; b = b0

citac_it = 0; err = 0 x = 0.

! --- tisk hlavicky ---

if(TISK > 0) print *,"c.kroku, x, f(x)"

! --- hl. cyklus --- do

(10)

citac_it = citac_it+1 citac_tisk = citac_tisk+1 fa = funkce(a)

fb = funkce(b)

! --- deleni intervalu podle zvol. metody --- if(metoda == 1) then ! puleni

x = (a + b) / 2.0 else ! secny

if(abs(fb - fa) < PRESNOST) then

! ve jmenovateli je 0 print *,"chyba: deleni nulou!"

err = 2 exit end if

x = (a * fb - b * fa) / (fb - fa) end if

fx = funkce(x)

! --- je nalezen koren? --- if(abs(fx) < PRESNOST) exit

! --- nastaveni mezi pro dalsi krok --- if(fx*fa < 0.0) then

b = x else

a = x end if

! --- tisk mezivysledku (pokud je nastaveno) --- if(TISK > 0 .and. citac_tisk == TISK) then

print *,citac_it,x,fx ! tisk mezivysledku citac_tisk = 0

end if

! --- ukonceni vyp., je-li presazen max. pocet it. --- if(citac_it == MAX_POCET_IT) then ! chyba: max. pocet it.

(11)

print *,"chyba: prekrocen max. pocet iteraci!"

err = 1 exit end if

! --- dalsi krok --- end do

vysl = x

print *,"Pocet iteraci:",citac_it end subroutine separace

Prostá iterační metoda

subroutine iterace(funkce,g_funkce,x_vstup,vysl,err) external funkce

real(8) :: funkce

external g_funkce ! upravena funkce real(8) :: g_funkce

real(8),intent(in) :: x_vstup real(8),intent(out) :: vysl integer,intent(out) :: err real(8) :: x,x0,fx

real(8),parameter::PRESNOST = 1e-10 integer,parameter::MAX_POCET_IT = 10000 integer,parameter::TISK = 1

! --- inicializace --- x0 = x_vstup

citac_it = 0; err = 0

(12)

do

citac_it = citac_it+1 citac_tisk = citac_tisk+1

! --- je nalezen koren? --- fx = funkce(x0)

if(abs(fx) < PRESNOST) exit

! --- tisk mezivysledku ---

if(TISK > 0 .and. citac_tisk == TISK) then print *,citac_it,x0,fx ! tisk mezivysledku citac_tisk = 0

end if

! --- ukonceni vyp., je-li presazen max. pocet it. --- if(citac_it == MAX_POCET_IT) then

err = 1 exit end if

! --- dalsi krok --- x0 = g_funkce(x0) end do

vysl = x0

print *,"Pocet iteraci:",citac_it end subroutine iterace

(13)

Newtonova metoda s analytickou derivací

subroutine newton(funkce,derivace,x_vstup,vysl,err) external funkce

real(8) :: funkce

external derivace ! f’(x) real(8) :: derivace

real(8),intent(in) :: x_vstup real(8),intent(out) :: vysl integer,intent(out) :: err real(8) :: x,x0,fx

real(8),parameter::PRESNOST = 1e-10 integer,parameter::MAX_POCET_IT = 10000 integer,parameter::TISK = 1

citac_it = 0; err = 0

do

! --- odhad korene ---

x = x0 - funkce(x0)/derivace(x0) citac_it = citac_it+1

citac_tisk = citac_tisk+1

(14)

! --- je nalezen koren? --- fx = funkce(x)

if(TISK > 0 .and. citac_tisk == TISK) then print *,citac_it,x,fx ! tisk mezivysledku citac_tisk = 0

end if

err = 1 exit end if

! --- dalsi krok --- x0 = x

end do vysl = x

print *,"Pocet iteraci:",citac_it end subroutine newton

Newtonova metoda s numerickou derivací

subroutine newton2(funkce,x_vstup,vysl,err)

external funkce ! funkce f(x), jejiz koreny

real(8) :: funkce ! hledam

real(8),intent(in) :: x_vstup ! poc. odhad korene real(8),intent(out) :: vysl ! nalezeny koren integer,intent(out) :: err ! chyba (0 .. OK,

! 1 .. prekrocen max. pocet

(15)

! kroku) real(8) :: x,x0,x1,fx,fx0,fx1,n_derivace

real(8),parameter::PRESNOST = 1e-10 ! pozad. presnost real(8),parameter::KROK = 0.01 ! x1 - x0

integer,parameter::MAX_POCET_IT = 10000 ! po kolika krocich se

! ma automat. ukoncit vypocet integer,parameter::TISK = 1 ! po kolika krocich tisk

! vysledku (0 = bez tisku

! mezivysl.)

x1 = x_vstup + KROK citac_it = 0; err = 0

do ! --- provadim k-ty krok ---

! x ... x(k)

! x0 ... x(k-1)

! x1 ... x(k-2)

! --- odhad korene --- fx0 = funkce(x0)

fx1 = funkce(x1)

n_derivace = (fx1-fx0)/(x1-x0) ! odhad derivace x = x0 - funkce(x0)/n_derivace

citac_it = citac_it+1 citac_tisk = citac_tisk+1

(16)

! --- je nalezen koren? --- fx = funkce(x)

if(TISK > 0 .and. citac_tisk == TISK) then print *,citac_it,x,fx ! tisk mezivysledku citac_tisk = 0

end if

err = 1 exit end if

! --- dalsi krok --- x0 = x1

x1 = x end do vysl = x

print *,"Pocet iteraci:",citac_it end subroutine newton2

Hlavní program

program num_reseni implicit none

(17)

integer :: vst_funkce,vst_metoda character :: vst_pokrac

logical :: pokracuj integer :: chyba

real(8) :: vysledek,h_mez,d_mez,odhad external f

real(8) :: f external g real(8) :: g external f2 real(8) :: f external g2 real(8) :: g external df real(8) :: df external dg real(8) :: dg real(8) :: PI PI = acos(-1.0) vysledek = 0.0 pokracuj = .true.

do while(pokracuj == .true.)

print *,"Kterou funkci chcete resit?"

print *

print *,"[1] f(x)=8x^3-6x-1"

print *,"[2] f(x)=x-sin(x)-pi/4"

print *

read *,vst_funkce print *

print *,"S pouzitim jake metody?"

print *

print *,"[1] separace - puleni intervalu"

print *,"[2] separace - metoda secen"

(18)

print *,"[3] Newtonova metoda"

print *,"[4] Newtonova metoda s numerickou derivaci"

print *,"[5] jednoducha iteracni metoda"

print *

read *,vst_metoda

! --- vstup poc. hodnot ---

if(vst_metoda == 1 .or. vst_metoda == 2) then print *

print *,"Dolni mez intervalu ?"

read *,d_mez print *

print *,"Horni mez intervalu ?"

read *,h_mez else

print *

print *,"Pocatecni odhad reseni ?"

read *,odhad end if

print *

! --- vypocet ---

select case(vst_metoda) case(1,2)

if(vst_funkce == 1) then

call separace(f,d_mez,h_mez,vst_metoda,vysledek,chyba) else

call separace(g,d_mez,h_mez,vst_metoda,vysledek,chyba) end if

case(3)

call newton(f,df,odhad,vysledek,chyba) else

call newton(g,dg,odhad,vysledek,chyba)

(19)

end if case(4)

call newton2(f,odhad,vysledek,chyba) else

call newton2(g,odhad,vysledek,chyba) end if

case(5)

call iterace(f,f2,odhad,vysledek,chyba) else

call iterace(g,g2,odhad,vysledek,chyba) end if

case default

print *,"Spatne zadane cislo metody!"

end select

! --- vypis nalezeneho vysledku --- if(chyba == 0) then

print *,"nalezeno x = ",vysledek else

print *,"Koren nenalezen."

end if

! --- pokracovat ? --- print *

print *,"Novy vypocet ? (a/n)"

read *,vst_pokrac

if(vst_pokrac == "N" .or. vst_pokrac == "n") pokracuj = .false.

end do

end program num_reseni

Zadané funkce

function f(x)

(20)

implicit none

real(8),intent(in) :: x real(8) :: f

f = 8.0 * x**3.0 - 6.0 * x - 1.0 end function f

function g(x) implicit none real(8) :: PI

real(8),intent(in) :: x real(8) :: g

PI = acos(-1.0)

g = x - sin(x) - PI/4.0 end function g

! --- upravy funkci pro prostou it. metodu --- function f2(x)

implicit none

real(8),intent(in) :: x real(8) :: f2

f2 = (8.0 * x**3.0 - 1.0) / 6.0 end function f2

function g2(x) implicit none real(8) :: PI

real(8),intent(in) :: x real(8) :: g2

PI = acos(-1.0) g2 = sin(x) + PI/4.0

(21)

end function g2

! --- derivace funkci pro Newtonovu metodu --- function df(x)

implicit none

real(8),intent(in) :: x real(8) :: df

df = 24 * x**2.0 - 6.0 end function df

function dg(x) implicit none real(8) :: PI

real(8),intent(in) :: x real(8) :: dg

PI = acos(-1.0) dg = 1 - cos(x) end function dg