|
|
| |
11. Numerické řešení rovnice f(x)=0 |
|
|
| |
- Motto: Každý numerický matematik má svou oblíbenou
metodu, která obecně nefunguje.
Důsledek: Následuje seznam oblíbených metod řešení
nelineárních rovnic.
Rada: Svoji oblíbenou metodu pochopíte jedině tak,
že ji budete krokovat v případě, že pro danou rovnici
funguje. Připravte si vhodnou zásobu funkcí pro levé
strany příslušných rovnic.
- Obecně je při řešení rovnice vhodné:
- neopustit definiční obor funkce f(x)
- použít maximálně zadaný počet vyčíslení
funkce f(x)
- končit výpočet při dosažení tolerance pro
f(x) nebo pro přírůstek x
- v případě nesplnění podmínek vracet
hodnotu inf
- Metody řešení nelineárních rovnic buď vylepšují
počáteční odhad a při tom riskují, nebo zužují interval
obsahující alespoň jeden kořen a přitom jsou pomalé.
Krásnou výjimkou je Riddersova metoda na konci
kapitoly.
- Pokud metoda vyžaduje derivace funkce v bodě, je
možné je vyčíslovat buď analyticky nebo numericky.
Vždy však jde o to, jak minimalizovat počet vyčíslení
funkce a jejích derivací.
- Poznámka: Teď jistě pochopíte, co je nudné na for
cyklu. Do těla cyklu stačí napsat jeden důležitý
vzorec, několik větvení s returnem a je hotovo. Pokud
to budete psát jinak, sice se nebudete tolik nudit, ale
bude to komplikovanější.
|
|
|
| |
Příklad 11A.1: |
|
| |
Vytvořte funkci pro realizaci prosté iterační
metody řešení rovnice x = F(x)
Instrukce: Krokujte a pak napište v textovém editoru vzorec pro
prostou iterační metodu. V případě problémů s konvergencí využijte
funkci Lipschitz k volbě lepšího intervalu a počáteční hodnoty
nebo k přeformulování úlohy. |
|
| |
Řešení: |
|
| |
| function [x]=Prosta(jmf,a,b,x0,nmax,epsilon) |
| % Obecne reseni rovnice x=F(x) metodou proste iterace |
| % [x]=Prosta(jmf,a,b,x0,nmax,epsilon); |
| % x ... nalezene reseni rovnice |
| % jmf ... nazev funkce F(x) jako retezec znaku |
| % a ... dolni mez definicniho oboru F(x) |
| % b ... horni mez definicniho oboru F(x) |
| % x0 ... pocatecni odhad reseni a<=x>=b |
| % nmax ... maximalni pocet vycisleni funkce F(x) |
| %epsilon ... tolerance rozdilu mezi levou a pravou stranou rovnice |
| for k=1:nmax |
% bezpecny cykl |
| x=feval(jmf,x0); |
% neprime volani funkce jmenem v bode x0 |
| if x<a | x>b |
% mimo definicni obor? |
| x=inf; |
% trapeni s konvergenci |
| return; |
% neslavny navrat z funkce |
| end |
| if abs(x-x0)<epsilon |
% postacuje presnost |
| return; |
% slavny navrat z funkce |
| end |
| x0=x; |
% zivot jde dal ... |
| end |
| x=inf; |
% problemy s konvergenci |
|
|
| |
| |
Příklad 11A.2: |
|
| |
Vytvořte funkci pro realizaci metody sečen pro
řešení rovnice f(x) = 0.
Instrukce: Krokujte a pak napište v textovém editoru
vzorec pro metodu sečen. Jde o nespolehlivou, ale superlineární
metodu řádu 1.618. |
|
| |
Řešení: |
|
| |
| function [x]=Secen(jmf,a,b,nmax,epsx,epsf) |
| % Obecne reseni rovnice f(x)=0 metodou secen |
| % [x]=Secen(jmf,a,b,nmax,epsx,epsf); |
| % x ... nalezene reseni rovnice |
| % jmf ... nazev funkce f(x) jako retezec znaku |
| % a ... dolni mez definicniho oboru f(x) |
| % b ... horni mez definicniho oboru f(x) |
| % nmax ... maximalni pocet vycisleni funkce f(x) |
| % epsx ... postacujici sirka intervalu <a;b> |
| % epsy ... postacujici hodnota funkce f(x) |
| if a>=b |
% zakladni kontrola |
| x=inf; |
% chybny predpoklad |
| return; |
| end |
| x1=a; |
| x2=b; |
| f1=feval(jmf,x1); |
% levy okraj |
| if abs(f1)<=epsf |
% proc se nezeptat? |
| x=x1; |
% hura!!! |
| return; |
| end |
| f2=feval(jmf,x2); |
% pravy okraj |
| if abs(f2)<=epsf |
% proc se nezeptat? |
| x=x2; |
% hura!!! |
| return; |
| end |
| if f1*f2>0 |
| x=inf; |
% chybny predpoklad |
| return; |
| end |
| for k=3:nmax |
% bezpecny cykl |
| x=(x1*f2-x2*f1)/(f2-f1); |
% kouzlo metody secen |
| if x<a | x>b |
% mimo definicni obor? |
| x=inf; |
% riziko metody secen |
| return; |
| end |
| f=feval(jmf,x); |
% vycislit v novem bode |
| if abs(f)<=epsf |
% proc se nezeptat? |
| return; |
| end |
| if abs(x-x2)<=epsx |
% staci to? |
| return; |
| end |
| x1=x2;f1=f2;x2=x;f2=f; |
% skatule hejbejte se |
| end |
|
|
| |
| |
Příklad 11A.3: |
|
| |
Vytvořte funkci pro realizaci Newtonovy
metody pro řešení rovnice f(x) = 0 za předpokladu, že derivace
je dána analyticky jako druhá funkce.
Instrukce: Krokujte a pak napište v textovém editoru vzorec pro
Newtonovu metodu. Jde o nespolehlivou metodu, ale druhého řádu. |
|
| |
Řešení: |
|
| |
| function [x]=Newton(jmf,jmd,a,b,x0,nmax,epsx,epsf) |
| % Obecne reseni rovnice f(x)=0 Newtonovou metodou |
| % [x]=Newton(jmf,jmd,a,b,x0,nmax,epsx,epsf); |
| % x ... nalezene reseni rovnice |
| % jmf ... nazev funkce f(x) jako retezec znaku |
| % jmd ... nazev funkce f'(x) jako retezec znaku |
| % a ... dolni mez definicniho oboru f(x) |
| % b ... horni mez definicniho oboru f(x) |
| % x0 ... pocatecni odhad reseni |
| % nmax ... maximalni pocet vycisleni funkce f(x) |
| % epsx ... postacujici sirka intervalu <a;b> |
| % epsy ... postacujici hodnota funkce f(x) |
| if a>=b |
% zakladni kontrola |
| x=inf; |
% chybny predpoklad |
| return; |
| end |
| for k=1:nmax |
% az do vycerpani trpelivosti |
| if x0<a | x0>b |
% mimo definicni obor? |
| x=inf; |
% riziko Newtonovy metody |
| return; |
| end |
| f=feval(jmf,x0); |
% vycislit f(x) v novem bode |
| if abs(f)<=epsf |
% proc se nezeptat? |
| x=x0; |
% hura!!! |
| return; |
| end |
| d=feval(jmd,x0); |
% vycislit derivaci v novem bode |
| if d==0 |
% nutny dotaz!!! |
| x=inf; |
% riziko Newtonovy metody |
| return; |
| end |
| x=x0-f/d; |
% tajemstvi Newtonovy metody |
| if abs(x-x0)<=epsx |
% proc se nezeptat? |
| return; |
| end |
| x0=x; |
% skatule hejbejte se |
| end |
|
|
| |
| |
Příklad 11A.4: |
|
| |
Vytvořte funkci pro realizaci Newtonovy metody
pro řešení rovnice f(x) = 0 za předpokladu, že derivace není známa.
Instrukce: Krokujte a pak napište v textovém editoru vzorec pro
Newtonovu metodu s numerickým výpočtem derivace. Jde o
nespolehlivou, ale superlineární metodu řádu 1.414. |
|
| |
Řešení: |
|
| |
| function [x]=Newton2(jmf,a,b,x0,nmax,epsx,epsf) |
| % Obecne reseni rovnice f(x)=0 Newtonovou metodou s numerickou derivaci |
| % [x]=Newton2(jmf,a,b,x0,nmax,epsx,epsf); |
| % x ... nalezene reseni rovnice |
| % jmf ... nazev funkce f(x) jako retezec znaku |
| % a ... dolni mez definicniho oboru f(x) |
| % b ... horni mez definicniho oboru f(x) |
| % x0 ... pocatecni odhad reseni |
| % nmax ... maximalni pocet vycisleni funkce f(x) |
| % epsx ... postacujici sirka intervalu <a;b> |
| % epsy ... postacujici hodnota funkce f(x) |
| if a>=b |
% zakladni kontrola |
| x=inf; |
% chybny predpoklad |
| return; |
| end |
| for k=1:floor(nmax/2) |
% proc asi |
| if x0<a | x0>b |
% mimo definicni obor? |
| x=inf; |
% riziko Newtonovy metody |
| return; |
| end |
| f=feval(jmf,x0); |
% vycislit v novem bode |
| if abs(f)<=epsf |
% proc se nezeptat? |
| x=x0; |
% hura!!! |
| return; |
| end |
| % numericka derivace usporne |
| d=(feval(jmf,x0+epsx)-f)/epsx; |
| if d==0 |
% nutny dotaz!!! |
| x=inf; |
% riziko Newtonovy metody |
| return; |
| end |
| x=x0-f/d; |
% tajemstvi Newtonovy metody |
| if abs(x-x0)<=epsx |
% proc se nezeptat? |
| return; |
| end |
| x0=x; |
% skatule hejbejte se |
| end |
|
|
| |
| |
Příklad 11A.5: |
|
| |
Vytvořte funkci pro realizaci Čebyševovy
metody pro řešení rovnice f(x) = 0 za předpokladu, že první i
druhá derivace jsou dány analyticky jako další funkce.
Instrukce: Krokujte a pak napište v textovém editoru vzorec pro
Čebyševovu metodu. Jde o nespolehlivou metodu, ale třetího řádu. |
|
| |
Řešení: |
|
| |
| function [x]=Cebysev(jmf,jmd,jmd2,a,b,x0,nmax,epsx,epsf) |
| % Obecne reseni rovnice f(x)=0 Cebysevovou metodou |
| % [x]=Cebysev(jmf,jmd,jmd2,a,b,x0,nmax,epsx,epsf); |
| % x ... nalezene reseni rovnice |
| % jmf ... nazev funkce f(x) jako retezec znaku |
| % jmd ... nazev funkce f'(x) jako retezec znaku |
| % jmd2 ... nazev funkce f''(x) jako retezec znaku |
| % a ... dolni mez definicniho oboru f(x) |
| % b ... horni mez definicniho oboru f(x) |
| % x0 ... pocatecni odhad reseni |
| % nmax ... maximalni pocet vycisleni funkce f(x) |
| % epsx ... postacujici sirka intervalu <a;b> |
| % epsy ... postacujici hodnota funkce f(x) |
| if a>=b |
% zakladni kontrola |
| x=inf; |
% chybny predpoklad |
| return; |
| end |
| for k=1:nmax |
% az do vycerpani trpelivosti |
| if x0<a | x0>b |
% mimo definicni obor? |
| x=inf; |
% riziko metody |
| return; |
| end |
| f=feval(jmf,x0); |
% vycislit f(x) v novem bode |
| if abs(f)<=epsf |
% proc se nezeptat? |
| x=x0; |
% hura!!! |
| return; |
| end |
| d=feval(jmd,x0); |
% vycislit derivaci v novem bode |
| if d==0 |
% nutny dotaz!!! |
| x=inf; |
% riziko metody |
| return; |
| end |
| d2=feval(jmd2,x0); |
% vycislit druhou derivaci v novem bode |
| x=x0-f/d-d2*f^2/2/d^3; |
% tajemstvi Cebysevovy metody |
| if abs(x-x0)<=epsx |
% proc se nezeptat? |
| return; |
| end |
| x0=x; |
% skatule hejbejte se |
| end |
|
|
| |
| |
Příklad 11A.6: |
|
| |
Vytvořte funkci pro realizaci
Halley-Richmondovy metody pro řešení rovnice f(x) = 0 za
předpokladu, že první i druhá derivace jsou dány analyticky
jako další funkce.
Instrukce: Krokujte a pak napište v textovém editoru vzorec pro
Halley-Richmondovu metodu. Jde o nespolehlivou metodu, ale
třetího řádu. |
|
| |
Řešení: |
|
| |
| function [x]=Richmond(jmf,jmd,jmd2,a,b,x0,nmax,epsx,epsf) |
| % Obecne reseni rovnice f(x)=0 Halley-Richmondovou metodou |
| % [x]=Richmond(jmf,jmd,jmd2,a,b,x0,nmax,epsx,epsf); |
| % x ... nalezene reseni rovnice |
| % jmf ... nazev funkce f(x) jako retezec znaku |
| % jmd ... nazev funkce f'(x) jako retezec znaku |
| % jmd2 ... nazev funkce f''(x) jako retezec znaku |
| % a ... dolni mez definicniho oboru f(x) |
| % b ... horni mez definicniho oboru f(x) |
| % x0 ... pocatecni odhad reseni |
| % nmax ... maximalni pocet vycisleni funkce f(x) |
| % epsx ... postacujici sirka intervalu <a;b> |
| % epsy ... postacujici hodnota funkce f(x) |
| if a>=b |
% zakladni kontrola |
| x=inf; |
% chybny predpoklad |
| return; |
| end |
| for k=1:nmax |
% az do vycerpani trpelivosti |
| if x0<a | x0>b |
% mimo definicni obor? |
| x=inf; |
% riziko metody |
| return; |
| end |
| f=feval(jmf,x0); |
% vycislit f(x) v novem bode |
| if abs(f)<=epsf |
% proc se nezeptat? |
| x=x0; |
% hura!!! |
| return; |
| end |
| d=feval(jmd,x0); |
% vycislit derivaci v novem bode |
| d2=feval(jmd2,x0); |
% vycislit druhou derivaci v novem bode |
| q=2*d^2-f*d2; |
% jmenovatel Richmondova vzorce |
| if q==0 |
% nutny dotaz!!! |
| x=inf; |
% riziko metody |
| return; |
| end |
| x=x0-2*f*d/q; |
% tajemstvi Richmondovy metody |
| if abs(x-x0)<=epsx |
% proc se nezeptat? |
| return; |
| end |
| x0=x; |
% skatule hejbejte se |
| end |
|
|
| |
| |
Příklad 11A.7: |
|
| |
Vytvořte funkci pro řešení rovnice f(x) = 0
metodou půlení intervalu za předpokladu, že funkce f(x) je
spojitá a a v krajních bodech intervalu mění znaménko.
Instrukce: Krokujte a pak napište v textovém editoru vzorec
pro metodu půlení intervalu. Jde o spolehlivou, ale pomalou metodu
prvního řádu. |
|
| |
Řešení: |
|
| |
| function [x]=Puleni(jmenof,a,b) |
| % Obecne reseni rovnice f(x)=0 metodou puleni intervalu |
| % [x]=Puleni(jmenof,a,b); |
| % x ... nalezene reseni rovnice |
| % jmenof ... nazev funkce f(x) jako retezec znaku |
| % a ... dolni mez intervalu |
| % b ... horni mez intervalu |
| nmax=40; |
% maximalni pocet vycisleni |
| fa=feval(jmenof,a); |
% levy okraj |
| fb=feval(jmenof,b); |
% pravy okraj |
| if fa*fb>0 |
| x=inf; |
% chybny predpoklad uspechu |
| return; |
| end |
| for k=3:nmax |
% bezpecny cykl |
| x=(a+b)/2; |
% stred intervalu |
| f=feval(jmenof,x); |
% vycislit v polovine |
| if fa*f<=0 |
% je neco v intervalu <a;x>? |
| b=x; |
% zahozeni intervalu (x;b> |
| else |
a=x; |
% zahozeni intervalu <a;x) |
| end |
| end |
| x=(a+b)/2; |
% naposled |
|
|
| |
| |
Příklad 11A.8: |
|
| |
Vytvořte funkci pro řešení rovnice
f(x) = 0 metodou Regula Falsi za předpokladu, že
funkce f(x) je spojitá a v krajních bodech intervalu
mění znaménko.
Instrukce: Krokujte a pak napište v textovém editoru
vzorec pro metodu půlení intervalu. Jde o spolehlivou, ale někdy
pomalou metodu prvního řádu. |
|
| |
Řešení: |
|
| |
| function [x]=RegulaFalsi(jmenof,a,b,nmax,epsx,epsf) |
| % Obecne reseni rovnice f(x)=0 metodou Regula Falsi |
| % [x]=RegulaFalsi(jmenof,a,b,nmax,epsx,epsf); |
| % x ... nalezene reseni rovnice |
| % jmenof ... nazev funkce f(x) jako retezec znaku |
| % a ... dolni mez definicniho oboru f(x) |
| % b ... horni mez definicniho oboru f(x) |
| % nmax ... maximalni pocet vycisleni funkce f(x) |
| % epsx ... postacujici sirka intervalu <a;b> |
| % epsy ... postacujici hodnota funkce f(x) |
| fa=feval(jmenof,a); |
% levy okraj |
| if abs(fa)<=epsf |
% proc se nezeptat? |
| x=a; |
% hura!!! |
| return; |
| end |
| fb=feval(jmenof,b); |
% pravy okraj |
| if abs(fb)<=epsf |
% proc se nezeptat? |
| x=b; |
% hura!!! |
| return; |
| end |
| if fa*fb>0 |
| x=inf; |
% chybny predpoklad |
| return; |
| end |
| for k=3:nmax |
% bezpecny cykl |
| x=(a*fb-b*fa)/(fb-fa); |
% kouzlo metody Regula Falsi |
| f=feval(jmenof,x); |
% vycislit v novem bode |
| if abs(f)<=epsf |
% proc se nezeptat? |
| return; |
| end |
| if fa*f<=0 |
% je neco v intervalu <a;x>? |
| b=x; |
% zahozeni intervalu (x;b> |
| fb=f; |
% aktualizace krajni hodnoty |
| else |
| a=x; |
% zahozeni intervalu >a;x) |
| fa=f; |
% aktualizace krajni hodnoty |
| end |
| if abs(b-a)<=epsx |
% uzky interval? |
| x=(a+b)/2; |
% to nic nestoji |
| return; |
| end |
| end |
| x=(a+b)/2; |
% pro jistotu |
|
|
| |
| |
Příklad 11A.9: |
|
| |
Vytvořte funkci pro řešení rovnice f(x) = 0
Riddersovou metodou za předpokladu, že funkce f(x) je spojitá
a v krajních bodech intervalu mění znaménko.
Instrukce: Krokujte a pak napište v textovém editoru vzorce
pro Riddersovu metodu. Jde o spolehlivou a navíc superlineární
metodu řádu 1.414. |
|
| |
Řešení: |
|
| |
| function [x]=Ridders(jmenof,a,b,nmax,epsx,epsf) |
| % Obecne reseni rovnice f(x)=0 Riddersovou metodou |
| % [x]=Ridders(jmenof,a,b,nmax,epsx,epsf); |
| % x ... nalezene reseni rovnice |
| % jmenof ... nazev funkce f(x) jako retezec znaku |
| % a ... dolni mez definicniho oboru f(x) |
| % b ... horni mez definicniho oboru f(x) |
| % nmax ... maximalni pocet vycisleni funkce f(x) |
| % epsx ... postacujici sirka intervalu <a;b> |
| % epsy ... postacujici hodnota funkce f(x) |
| if a>=b |
% zakladni kontrola |
| x=inf; |
% chybny predpoklad |
| return; |
| end |
| fa=feval(jmenof,a); |
% levy okraj |
| if abs(fa)<=epsf |
% proc se nezeptat? |
| x=a; |
% hura!!! |
| return; |
| end |
| fb=feval(jmenof,b); |
% pravy okraj |
| if abs(fb)<=epsf |
% proc se nezeptat? |
| x=b; |
% hura!!! |
| return; |
| end |
| if fa*fb>0 |
| x=inf; |
% chybny predpoklad |
| return; |
| end |
| for k=2:floor(nmax/2) |
% proc? |
| s=(a+b)/2; |
% zacina to obycejnym pulenim |
| fs=feval(jmenof,s); |
% vycislit uprostred |
| if abs(fs)<=epsf |
% proc ne? |
| return; |
| end |
| % tajemstvi Riddersovy metody |
| x=s+(s-a)*sign(fa-fb)*fs/sqrt(fs^2-fa*fb); |
| fx=feval(jmenof,x); |
% vycislit v lepsim bode |
| if abs(fx)<=epsf |
% proc ne? |
| return; |
| end |
| if x<s |
% je odhad x v leve polovine intervalu? |
| if fx*fs<0 |
% kdyz ano |
| a=x;fa=fx;b=s;fb=fs; |
% posun obou mezi |
| else |
| b=x;fb=fx; |
% posun horni meze |
| end |
| else |
| if fx*fs<0 |
% kdyz ne |
| a=s;fa=fs;b=x;fb=fx; |
% posun obou mezi |
| else |
| a=x;fa=fx; |
% posun dolni meze |
| end |
| end |
| if abs(b-a)<=epsx |
% uzky interval? |
| x=(a+b)/2; |
% vylepseni |
| return; |
| end |
| end |
| x=(a+b)/2; |
% pro jistotu |
|
|
| |
| |
Příklady pro samostatné vypracování |
|
| |
Příklad 11B.1: |
|
| |
Pro potřeby metod využívajících
analytický výpočet derivací (Newton, Čebyšev, Richmond)
vytvořte funkci, která z hodnoty x vypočte vektor tvořený
hodnotou funkce a jejich derivací. Pak upravte uvedené metody
tak, aby byly schopné takovou funkci obecně volat. |
|
| |
| |
Příklad 11B.2: |
|
| |
Všechny funkce doplňte tak, aby vracely
také hodnotu funkce a počet vyčíslení funkce. |
|
| |
| |
Příklad 11B.3: |
|
| |
Porovnejte efektivitu jednotlivých metod
řešení rovnic. |
|
| |
| |
Příklad 11B.4: |
|
| |
Vytvořte funkci, která zdůrazní klady metody Regula
Falsi a newtonovy metody. |
|
| |
| |
Příklad 11B.5: |
|
| |
Vytvořte funkci, která zdůrazní klady
metody půlení intervalu. |
|
| |
| |
Příklad 11B.6: |
|
| |
Vytvořte funkci, která zdůrazní klady
Riddersovy metody. |
|
| |
| |
Seznam použitých příkazů |
|
| |
|
|
|
|
|