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 typu f (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]. 1
Řešit rovnici (1) numericky znamená navrhnout algoritmus přibližného řešení, při kterém získáme posloupnost aproximací x0 , x1 , ..., xk , ...
(3)
lim xk = ξ
(4)
takovou, že k→∞
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í aproximace x0 . 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 < a0 , b0 >, ve kterém se nachází hledaný kořen. Pro krajní body intervalu musí platit f (a0 )f (b0 ) < 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 : x1 =
a0 + b0 2
(6)
Poté testujeme, zda x1 neodpovídá hledanému řešení: |f (x1 )| < p 2
(7)
kde p je požadovaná přesnost. Není-li tato podmínka splněna, testujeme zda f (a0 )f (x1 ) < 0
(8)
V takovém případě se kořen nachází v intervalu < a0 , x1 >. 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 f (b) − f (a) y= (x − a) + f (a) (9) b−a Aproximace řešení xk se získá z rovnice (9), položíme-li y = 0: xk =
ak−1 f (bk−1 ) − bk−1 f (ak−1 ) f (bk−1 ) − f (ak−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 neobejde bez „dobréhoÿ odhadu počáteční aproximace x0 . 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 ) Prostá iterační metoda nemusí vždy konvergovat.
3
(12)
2.4
Newtonova metoda
Jde opět o metodu iterační vyžadující odhad počáteční aproximace x0 . Aproximace ř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 0 (xk )
(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ě xk nahradíme přímkou procházející body xk a xk−1 . Je-li x0 počáteční aproximací, hodnotu x1 určíme podle vztahu x1 = x0 (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 0 (x) nahradíme výrazem f (xk ) − f (xk−1 ) f 0 (x) ≈ (15) xk − xk−1 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: 8x3 − 6x − 1 = 0 π x − sin x − = 0 4
(16) (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 Newtonovou metodou. U Newtonovy metody jsme zkoumali dvě verze – s analytickým výpočtem derivace f 0 (x) a s jeho numerickým přiblížením podle vzorce (15). Přesnost výpočtu byla nastavena na hodnotu ² = 10−10 a výpis programů v jazyce Fortran je v Příloze. 4
3.1
Metoda půlení intervalu
Metodou půlení intervalu byly nalezeny tyto kořeny: Funkce f1 (x) = 8x3 − 6x − 1 = 0 poč. interval < −1, −0.5 > < −0.5, 0 > < 0.5, 1.5 >
x −0.766044443123974 −0.173648177675204 0.939692620784626
Funkce f2 (x) = x − sin x − poč. interval < 0, 3.14 >
π 4
f1 (x) −4.03855113413620 · 10−11 4.36562151446412 · 10−11 −1.94774608767989 · 10−11
počet iterací 32 32 36
=0
x 1.76634028665372
f2 (x) 6.07372557738690 · 10−11
počet iterací 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) = 8x3 − 6x − 1 = 0 poč. interval < −1, −0.5 > < −0.5, 0 > < 0.5, 1.5 >
x −0.766044443110925 −0.173648177668642 0.939692620779913
Funkce f2 (x) = x − sin x − poč. interval < 0, 3.14 >
π 4
f1 (x) 6.50997606516102 · 10−11 9.03027197289574 · 10−12 −9.10867358101808 · 10−11
počet iterací 25 11 37
=0
x 1.76634028652257
f2 (x) −9.58941248890149 · 10−11 5
počet iterací 21
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 f1 (x) = 8x3 − 6x − 1 = 0 Obě funkce jsme upravili podle vzorce (11). U funkce f1 (x) jsme testovali dva tvary rovnice (11): 8x3 − 1 6 1√ 3 x = 6x + 1 2
x =
(18) (19)
S použitím funkce (18) jsme došli k následujícímu řešení: poč. odhad −0.25
x −0.173648177676812
f1 (x) 5.21391611526048 · 10−11
počet iterací 12
Použijeme-li funkci (19), metoda konvergovat nebude. Funkce f2 (x) = x − sin x −
π 4
=0
Funkci f2 (x) jsme upravili na tvar x = sin(x) +
π 4
(20)
Výsledek výpočtu: poč. odhad 1.5
x 1.76634028665244
6
f2 (x) 5.92004918394465 · 10−11
počet iterací 14
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 f1 (x) = 8x3 − 6x − 1 = 0 poč. odhad −1 −0.25 0.75
x −0.766044443118984 −0.173648177666930 0.939692620785911
Funkce f2 (x) = x − sin x − poč. odhad 1.5
3.5
π 4
f1 (x) −4.91854820761084 · 10−14 5.32343266690383 · 10−17 3.94601885889134 · 10−14
počet iterací 5 4 5
=0
x 1.76634028660287
f2 (x) 8.96000938374608 · 10−15
počet iterací 4
Newtonova metoda s numerickou derivací
Po úpravě vzorce (13) pomocí vztahu (15) jsme dostali následující výsledky: Funkce f1 (x) = 8x3 − 6x − 1 = 0 poč. odhad −1 −0.25 0.75
x −0.766044443118978 −0.173648177666931 0.939692620785909
Funkce f2 (x) = x − sin x − poč. odhad 1.5
π 4
f1 (x) −3.41393580072236 · 10−15 3.42152521592975 · 10−15 6.66133814775094 · 10−16
počet iterací 7 5 7
=0
x 1.76634028660001
f2 (x) −3.40559947863833 · 10−12 7
počet iterací 5
Shrnutí Newtonova metoda opět nemusí vždy konvergovat. Pokud však konverguje, konverguje velmi rychle (šlo o nejrychlejší metodu v testu). Nevýhodou klasické Newtonovy metody je nutnost analyticky derivovat vstupní funkci. Tuto nevýhodu však lze eliminovat náhradou derivace f 0 (x) jejím numerickým přiblíž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 pravdě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 studovaných rovnic. Jak ilustrují výše uvedené tabulky, úspěch při numerickém řešení rovnice f (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).
8
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 integer, intent(in):: metoda real(8), intent(out) :: vysl integer, intent(out) :: err
! ! ! ! ! ! ! !
poc. interval pouzita metoda, 0=puleni, 1=secny nalez. koren 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 9
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. 10
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 integer :: citac_it integer :: citac_tisk 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 11
! --- tisk hlavicky --if(TISK > 0) print *,"c.kroku, x, f(x)" 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 print *,"chyba: prekrocen max. pocet iteraci!" err = 1 exit end if ! --- dalsi krok --x0 = g_funkce(x0) end do vysl = x0 print *,"Pocet iteraci:",citac_it end subroutine iterace
12
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 integer :: citac_it integer :: citac_tisk
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 ! --- tisk hlavicky --if(TISK > 0) print *,"c.kroku, x, f(x)" do ! --- odhad korene --x = x0 - funkce(x0)/derivace(x0) citac_it = citac_it+1 citac_tisk = citac_tisk+1
13
! --- je nalezen koren? --fx = funkce(x) if(abs(fx) < PRESNOST) exit
! --- tisk mezivysledku --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 print *,"chyba: prekrocen max. pocet iteraci!" 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 real(8),intent(out) :: vysl integer,intent(out) :: err
14
! poc. odhad korene ! nalezeny koren ! chyba (0 .. OK, ! 1 .. prekrocen max. pocet
! kroku) real(8) :: x,x0,x1,fx,fx0,fx1,n_derivace integer :: citac_it integer :: citac_tisk
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.)
! --- inicializace --x0 = x_vstup x1 = x_vstup + KROK citac_it = 0; err = 0 ! --- tisk hlavicky --if(TISK > 0) print *,"c.kroku, x, f(x)" 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) x = x0 - funkce(x0)/n_derivace citac_it = citac_it+1 citac_tisk = citac_tisk+1 15
! odhad derivace
! --- je nalezen koren? --fx = funkce(x) if(abs(fx) < PRESNOST) exit
! --- tisk mezivysledku --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 print *,"chyba: prekrocen max. pocet iteraci!" 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
16
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 print print print
* *,"S pouzitim jake metody?" * *,"[1] separace - puleni intervalu" *,"[2] separace - metoda secen" 17
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) if(vst_funkce == 1) then call newton(f,df,odhad,vysledek,chyba) else call newton(g,dg,odhad,vysledek,chyba) 18
end if case(4) if(vst_funkce == 1) then call newton2(f,odhad,vysledek,chyba) else call newton2(g,odhad,vysledek,chyba) end if case(5) if(vst_funkce == 1) then 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) 19
implicit none real(8),intent(in) :: x real(8) :: f f = 8.0 * x**3.0 - 6.0 * x end function f
- 1.0
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 end function f2
- 1.0) / 6.0
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 20
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
21
