# Gyakorlat 2019. 03. 20.
close all
output_precision(5)

########
# Az Octave grafikai motorjat lehet itt kivalasztani
# ha gondok vannak az abrazolassal erdemes probalni mas
# beallitast, hatha azzal sikerul
########
# graphics_toolkit qt
# graphics_toolkit fltk
 graphics_toolkit gnuplot

disp("Gyakorlat 2019. 03. 20.");

disp("\n\
Nemlinaris egyenlet numerikus megoldasa: f(x) = 0, x = ?")

disp("\n\
1. Pelda: cos(x) + x = 0")

disp("\n\
Erdemes abrazolni a fuggvenyt es ebbol\n\
levonni valamilyen kovetkeztetest,\n\
sejtest arra, hogy mi a megoldas.")


disp("\n\
Fuggveny abrazolasnal klasszikusan egy\n\
intervallumon abrazolunk, legyen pl [-2, 2]\n\
s ezen az intervallumon veszunk fel egy rakas\n\
pontot amiken mind kiertekeljuk a fuggvenyt.\n\
Az igy kapott pontokat a szamitogep valamilyen\n\
logika alapjan osszekoti.")

disp("\n\
x -2 es 2 kozotti intervallum")
x = [-2 : 1 : 2]

disp("\n\
Az intervallumot 1-es lepessel 5 pontra bontottuk fel, ezen\n\
az 5 ponton fogjuk kiertekelni a fuggvenyt")
y = cos(x) + x

disp("\n\
Majd abrazoljuk a plot paranccsal: Figure 1")

# uj kep nyitasa, nem kotelezo csak itt most
# tobb kepet fogunk csinalni es nem akarjuk felulirni
figure (1);

# adatpontok kirajzolasa es osszekotese
plot(x,y);

# adunk egy cimet a kepnek
title("cos(x) + x abrazolasa [-2, 2]-n 5 ponttal")

disp("\n\
Lathatoan egyenes vonalakkal koti ossze az adatpontokat.")

disp("\n\
Hogy \"simabban\", pontosabban lassuk a fuggvenyt,\n\
vegyunk fel tobb adatpontot: x = [-2 : 0.01 : 2];\n\
Itt mar tegyunk ;-t a sor vegere, hogy ne irja ki\n\
az Octave az x erteket a kepernyore mert az felesleges\n\
es lassit is. Abra: Figure 2")
x = [-2 : 0.01 : 2];
y = cos(x) + x;

figure(2);
plot(x,y);
title("cos(x) + x abrazolasa [-2 : 0.01 : 2]-n")

disp("\n\
a \"whos x\" parancs mutatja, hogy 401 adatpontot vettunk fel\n")

whos x

disp("\n\
Latszik, hogy valahol -1 es 0 kozott\n\
van egy zerushelye a fuggvenynek.")

disp("\n\
Az Octave az \"fsolve\" eljarast ajanlja alapbol az f(x) = 0\n\
egyenlet megoldasara.\n\n\
Ennek szintaxisa: fsolve(f, x0, opciok)\n\n\
Reszletes leirasert lasd: \"help fsolve\" illetve\n\
Google kereses: Octave fsolve")

disp("\n\
Mi most csak az alapokkal ismerkedunk meg,\n\
f a kerdeses fuggveny, x0 a kezdeti becsles.\n\
Az opciokban lehet finomhangolni a gyokkeresest\n\
amire gyakorlati alkalmazasoknal gyakran szukseg is van.")

disp("\n\
Elso feladatunk tehat a fuggveny definialasa.")

disp("\n\
Erre tobb modot nezunk meg gyorsan.\n\
a)  kulon file: File -> New -> New Function ->\n\
                itt adunk neki egy nevet, pl Pelda1_a")
               
disp("\n\n\
    A lenyegi resze a megnyilo uj fuggvenyfilenak:\n\n\
    function retval = Pelda1_a (input1, input2)\n\
\n\
    endfunction")
    
disp("\n\n\
    Ezt a reszt a kovetkezo modon toltjuk ki:\n\n\
    function retval = Pelda1_a (x)\n\
      retval = cos(x) + x;\n\
    endfunction")

disp("\n\
    Ezt a file-t most letre kell hozni kulonben\n\
    ennek a scriptnek a futasa megszakad.")

disp("\n\
    Pelda1_a tesztelese:")
Pelda1_a(3)

disp("\n\n\
b)  kulon file helyett, ugyanez a fuggveny deklaracio\n\
    itt sorfolytonosan is megtortenhet:\n\n\
    function retval = Pelda1_b (x)\n\
      retval = cos(x) + x;\n\
    endfunction")
function retval = Pelda1_b (x)
  retval = cos(x) + x;
endfunction

disp("\n\
    Pelda1_b tesztelese:")
Pelda1_b(3)

disp("\n\n\
c)  Egyszerubb fuggvenyekre meg rovidebb megoldas az\n\
    ugynevezett lambda deklaracio.\n\
    Kicsit csunyanak nez ki de sokat lehet vele sporolni:\n\
    Pelda1_c = @(x) cos(x) + x\n\n\
    Az @(x) fejezi ki, hogy fuggvenyrol\n\
    van szo aminek x a valtozoja\n\
    Ezt koveti szokoz utan a fuggveny keplete.")
Pelda1_c = @(x) cos(x) + x;
    
disp("\n\
    Pelda1_c tesztelese:")

Pelda1_c(3)

disp("\n\
=================\n\
Az igy definialt fuggvenyt is lehet grafikon\n\
rajzolasara hasznalni: Figure 3")

# uj kep
figure(3);

# a kep 1x3 abrabol all, az 1. az aktiv
subplot(1,3,1);

plot(x, Pelda1_a(x));
title("Pelda1_a");

subplot(1,3,2);

plot(x, Pelda1_b(x));
title("Pelda1_b");

subplot(1,3,3);

plot(x, Pelda1_c(x));
title("Pelda1_c");

disp("\n\
===========\n\
A megoldashoz barmelyik modon megadott fuggvenyt hasznalhatjuk.\n\n\
A grafikon kirajzolasa egy gyokot sejtet -1 korul, \n\
legyen, hat a kezdo becslesunk -1.\n\n\
Az a) es b) modon megadott fuggvenyek es a c) modon\n\
megadott verzio kozott egy lenyeges kulonbseget\n\
azert ki kell emelni:\n\
- mashogy lehet rajuk hivatkozni.\n\n\
Az a) es b) modon megadott fuggvenyekre idezojelbe tett\n\
nevvel tudunk hivatkozni, mig a c)-ben a fuggvenyt\n\
tulajdonkeppen egy valtozoba helyeztuk el\n\
    (Pelda1_c = @(x) cos(x) + x => Pelda1_c valtozo erteke a fgv)\n\
igy itt a valtozo neve hasznalando:\n\n\
fsolve( \"Pelda1_a\", -1 )")
fsolve( "Pelda1_a", -1)

disp("\n\
fsolve( \"Pelda1_b\", -1 )")
fsolve( "Pelda1_b", -1)

disp("\n\
fsolve( Pelda1_c, -1 )")
fsolve( "Pelda1_b", -1)

disp("\n\
Az fsolve a gyokot alapbeallitasokkal \n\
a -0.73909 ertekkel kozeliti, nezzuk erre mit ad a fgv:\n\
Pelda1_c(-0.73909) = ")
Pelda1_c(-0.73909)

disp("\n\
==============\n\n\
2. Pelda - Technika laborrol")

disp("\n\
10q - 10 = q^1.2124 - 1      => q = ?\n\
(nem a q = 1 az erdekes)\n\n\
Pelda2 = @(x) 10 * x - 10 - x.^1.2124 + 1")
# a .^ arra utal, hogy ha x matrix/vektor akkor
# elemenkent hatvanyozunk.
#
# erre akkor van elsosorban szukseg mikor az itt
# definialt fuggvenyt abrazoljuk, hiszen akkor
# a fuggveny bemenete egy x = [a : lepes : b]
# tartomany ami egy sorvektor es ebbol akarunk
# kapni egy masik sorvektort ami pontonkent
# tartalmazza a fgv ertekeket
Pelda2 = @(x) 10 * x - 10 - x.^1.2124 + 1

disp("\n\
A tort hatvany miatt nemnegativ x-ekkel fogunk dolgozni.\n\n\
Nezzuk a fuggvenyt [0, 100]-on: Figure 4")
figure(4);
x = [0 : 0.01 : 100];
plot(x, Pelda2(x));
title("Pelda2 [0 : 0.01 : 100]-on")

disp("\n\
Nem latunk masik gyokot, mint 0 korul valamit (ez az 1)\n\
vizsgaljunk joval nagyobb intervallumot: [0, 30000]: Figure 5\n\n\
(turelem kicsit lassu lehet windowson)")
figure(5);
x = [0 : 0.01 : 30000];
plot(x, Pelda2(x));
title("Pelda2 [0 : 0.01 : 30000]-en")

disp("\n\
Latjuk, hogy gorbul vissza a fuggveny, menjunk meg\n\
joval tovabb: [0, 60000]: Figure 6\n\n\
(turelem kicsit lassu lehet windowson)")
figure(6);
x = [0 : 0.01 : 60000];
plot(x, Pelda2(x));
title("Pelda2 [0 : 0.01 : 60000]-en")

disp("\n\
Latszik, hogy 50000 korul is van egy zerohely.\n\
Ugy kepzeljunk el egy f(x) = 0 -t megoldo modszert\n\
mint ami valahogy kihasznalva a fgv lejteset halad\n\
a zerushely fele.\n\n\
Az fsolve hasznalatanal igy jol kell megvalasztani\n\
a kezdeti becslest. Ha a \"pukli\" baloldalarol indulunk\n\
jo esellyel az x = 1 -es megoldashoz jutunk,\n\
mig ha a jobb oldalrol indulunk akkor erkezunk\n\
az x ~ 50000-es zerohelyhez:")

disp("\n\
10000-es kezdeti becslessel:\n\n\
fsolve(Pelda2, 10000) (... turelem)")
fsolve(Pelda2, 10000)

disp("\n\
Megjegyzes: erdekes, hogy komplex szamkent kapjuk\n\
 a kozelitest. Ennek oka, hogy a tort hatvany fuggveny negativ\n\
 szamokon komplex erteket vesz fel es az fsolve altal\n\
 hasznalt modszer jo esellyel at megy neha\n\
 a negativ szamok koze is.\n\n\
 Mindenesetre a becsles kepzetes resze elhanyagolhato,\n\
 latszik, hogy ez kb 1 + 0i.")

disp("\n\
30000-es kezdeti becslessel:\n\n\
fsolve(Pelda2, 30000) (... turelem)")
fsolve(Pelda2, 30000)

disp("\n\
Nezzuk mit is talaltunk:\n\
Pelda2(51057.77512) = ")
Pelda2(51057.77512)

disp("\n\n\
============================\n")

disp("\n\
Bracketing (bezaro, kozrefogo) algoritmusok (bracket: [])")

disp("\n\
Ket algoritmusrol ejtunk szot:\n\
- Felezo algoritmus (bisection)\n\
- hu'r modszer (regula falsi)")

disp("\n\
Bracketing modszerek kozos jellemzok:\n\
- kell egy kiindulasi intervallum [a, b]\n\
  melyre teljesul, hogy f(a) es f(b) elojele\n\
  kulonbozo. Innen folytonos f fgv-re\n\
  a Bolzano tetelbol tudjuk, hogy [a, b]-n\n\
  van zerohelye, azaz letezik x ertek [a, b]-n,\n\
  hogy f(x) = 0\n\
- valamilyen modszerrel valasztunk [a, b]-n\n\
  egy koztes pontot: c\n\
- ekkor vagy f(c) = 0 es talaltunk egy megoldast\n\
- vagy f(c) > 0 vagy f(c) < 0\n\
  ekkor azt a vegpontot kicsereljuk c-re amelyikben\n\
  az f fgv elojele megegyezik f(c) elojelevel.\n\n\
  Azaz pl ha f(a) > 0, f(b) < 0 es f(c) > 0 akkor\n\
  az uj keresesi intervallumunk [c, b],\n\
  mig ha pl f(a) > 0, f(b) < 0 es f(c) < 0 akkor\n\
  az uj keresesi intervallumunk [a, c]")
  
disp("\n\
A lenyeg, hogy az uj, kisebb intervallumra\n\
ugyanugy teljesul, hogy a vegpontokon\n\
az elojelek kulonbozoek.\n\n\
Mivel az intervallumok osszehuzodnak es\n\
a Bolzano tetel feltetelei vegig teljesulnek,\n\
igy garantaltan kozelitunk egy gyokot.\n\n\
(Ha f nem folytonos, lehet, hogy szakadasi pontot,\n\
 polust kozelitunk, erre vigyazni kell.)")

disp("\n\
Az elteres a bracketing modszerek kozott\n\
abban nyilvanul meg, hogy milyen logika alapjan\n\
valasztja az algoritmus ezt a koztes c pontot [a, b]-n:\n\n\
1) Bisection / Felezo modszer\n\
c az [a, b] felezopontja, azaz c = (a + b) / 2\n\n\
-> elony:   intervallumok merete folyamatosan felezodik\n\
            igy relative gyorsan kozeliti a gyokot\n\n\n\
2) Regula Falsi / hu'r modszer\n\
c az (a, f(a)) es (b, f(b)) pontokat osszekoto hu'r\n\
metszespontja az x-tengellyel.\n\n\
-> elony:   altalaban gyorsabb mint a felezo modszer\n\
-> hatrany: neha egy vegpont \"be tud ragadni\", ekkor\n\
            az algoritmus kulonfele modositasai\n\
            valtoztatnak c kivalasztasan:\n\
            - pl beiktatnak egy felezo lepest\n\
            - vagy a beragadt pont fgv erteket\n\
              \"lesulyozzak\" (Illinois modszer)")

disp("\n\
===============\n\
3. Pelda: 5x^2 + 9x - 80\n\
          polinomnak van gyoke [3, 4]-en.\n\
          Keressuk meg felezo modszerrel!")
          
disp("Deklaraljuk a fuggvenyt:\n\
Pelda3 = @(x) 5 * x.^2 + 9 * x - 80")
Pelda3 = @(x) 5 * x.^2 + 9 * x - 80

disp("\n\
Abrazolas [3, 4]-en: Figure 7")

figure(7);
x = [3 : 0.01 : 4];
plot(x, Pelda3(x))
title("Pelda3 [3, 4]-en");

disp("\n\
Valoban Pelda3(3) < 0, hisz")
Pelda3(3)

disp("es Pelda3(4) > 0")
Pelda3(4)

disp("\n\
Igy a [3, 4]-en van zerohely.")

a = 3
b = 4

disp("\n\
Felezo modszer: [3, 4] felezopontja (3+4)/2=3.5")
c = (a + b) / 2

disp("\n\
Itt a fuggveny ertek:")
Pelda3(c)

disp("\n\
Ez pozitiv, igy b-t csereljuk le (mert f(b) > 0 volt): [3, 3.5]\n\
az uj bezaro intervallumunk.")
a
b = c

disp("\n\
Felezo modszer: [3, 3.5] felezopontja (3+3.5)/2=3.25")
c = (a + b) / 2

disp("\n\
Itt a fuggveny ertek:")
Pelda3(c)

disp("\n\
Ez pozitiv, igy b-t csereljuk le: [3, 3.25]\n\
az uj bezaro intervallumunk.")
a
b = c

disp("\n\
Felezo modszer: [3, 3.25] felezopontja (3+3.25)/2=3.125")
c = (a + b) / 2

disp("\n\
Itt a fuggveny ertek:")
Pelda3(c)

disp("\n\
Ez negativ, igy a-t csereljuk le (mert f(a) < 0 volt): [3.125, 3.25]\n\
az uj bezaro intervallumunk.")
a = c
b

disp("\n\
Felezo modszer: [3.125, 3.25] felezopontja (3.125+3.25)/2=3.1875")
c = (a + b) / 2

disp("\n\
Itt a fuggveny ertek:")
Pelda3(c)

disp("\n\
Ez negativ, igy a-t csereljuk le : [3.1875, 3.25]\n\
az uj bezaro intervallumunk.")
a = c
b

disp("\n\
stb...\n\n\
Kerdezzuk meg az Octave-ot a megoldasrol, \n\
amit mi kaptunk az jo is kezdeti becslesnek:\n\
fsolve(Pelda3, 3.1875)")
fsolve(Pelda3, 3.1875)

disp("\n\
===============\n\
Visszateres a 2. Peldahoz\n\n\
Toltsuk le a weboldalrol a Bisection.m es a RegulaFalsi.m\n\
programokat.")

disp("\n\
Emlekeztetoul:\n\
Pelda2 = @(x) 10 * x - 10 - x.^1.2124 + 1\n\n\
Aminek grafikonja [0, 60000]-en egy \"pukli\",\n\
ket gyoke van: 1 es kb 51057.8")

disp("\n\n\
Mindket letoltott program hivasa a kovetkezo:\n\
...( f, a, b, hosszTol, ertekTol, maxIter )\n\n\
ahol\n\
- f    :    a fuggveny\n\
- a, b :    az [a, b] intervallum ket vegpontja amire\n\
            f(a) es f(b) elojele kulonbozo\n\
- hosszTol: ha ennel rovidebbe valik a bezaro\n\
            intervallum a futas soran megallunk\n\
- ertekTol: ha ennel kisebbe valik az |f(c)|\n\
            ertek a futas soran megallunk\n\
- maxIter:  maximum ennyi lepest teszunk")

disp("\n\
Felezo modszer futtatasa:\n")

disp("\n\
[0, 100] -on van gyok:")
Pelda2(0)
Pelda2(100)

disp("\n\
Bisection(Pelda2, 0, 100, 0.01, 0.001, 4)")
Bisection(Pelda2, 0, 100, 0.01, 0.001, 4)

disp("\n\
Latszik, hogy a maximalis 4 lepes miatt alltunk meg, \n\
engedjuk tovabb:\n\
Bisection(Pelda2, 0, 100, 0.01, 0.001, 100)")
Bisection(Pelda2, 0, 100, 0.01, 0.001, 100)

disp("\n\
Latszik, hogy azert alltunk meg mert a bezaro\n\
intervallum hossza kisebb lett mint 0.01\n\n\
a visszaadott c becsult zerohelyen a fuggveny\n\
erteke -0.04 tehat nem kisebb abszolut ertekben mint\n\
0.001. Engedjuk kicsit tovabb futni az algoritmust:\n\
Bisection(Pelda2, 0, 100, 0.000001, 0.001, 100)")
Bisection(Pelda2, 0, 100, 0.000001, 0.001, 100)

disp("\n\
Itt mar a konkluzio, hogy azert allunk meg mert az\n\
utolso c ponton a fgv ertek abszolut ertekben kisebb\n\
mint az elore megadott tolerancia a 0.001.")

disp("\n\
=================\n\
[40000, 60000]-en van gyok:")
Pelda2(40000)
Pelda2(60000)

disp("\n\
Bisection(Pelda2, 40000, 60000, 0.0000001, 0.001, 100)")
Bisection(Pelda2, 40000, 60000, 0.0000001, 0.001, 100)

disp("\n\
=================\n\
Ugyanez a hu'r modszerrel:\n\
RegulaFalsi(Pelda2, 40000, 60000, 0.0000001, 0.001, 100)")
RegulaFalsi(Pelda2, 40000, 60000, 0.0000001, 0.001, 100)

disp("\n\
Itt eppen azt latjuk, hogy a hu'r modszer gyorsabban\n\
er veget.")

