# Gyakorlat 2019. 04. 24.
close all
output_precision(15)

disp("Gyakorlat 2019. 04. 24.");

disp("\n\
Polinomiaris interpolacio: Lagrange interpolacio")

disp("\n\n\
Alap problema:\n\
- y = f(x) ismeretlen fuggveny\n\
- adott n db meres, azaz (x1, y1) ... (xn, yn) (kulonbozo x-ek)\n\
  melyekre yi = f(xi)")

disp("\n\n\
Lagrange interpolacio: az n meresi ponthoz (xi, yi) meghatarozzuk\n\
azt az (n - 1)-ed foku p(x) polinomot, amire yi = f(xi) = p(xi)")

disp("\n\n\
Tetel: \n\
Letezik pontosan egy ilyen (n-1)-ed foku polinom\n\n\
Ez onnan latszik, hogy az n db egyenlet felirasa egy\n\
linearis egyenletrendszert ad az egyutthatokra:")

A = ["1 " "x1 " "x1^2 " "x1^3 " " ... " "x1^n-1 ";
"1 " "x2 " "x2^2 " "x2^3 " " ... " "x2^n-1 ";
"1 " "... " "... " "... " " ... " "...";
"1 " "xn " "xn^2 " "xn^3 " " ... " "xn^n-1 "]

b = ["y1"; "y2"; "..."; "yn"]

x = ["a0"; "a1"; "..."; "an-1"]

disp("\n\
Megkapjuk az egyutthatokat ha megoldjuk az Ax=b egyenletet\n\
Itt az egyutthatomatrix A egy ugynevezett Vandermonde matrix ami\n\
mindig invertalhato, igy mindig pontosan egy megoldas letezik.")

###########################################################################

disp("\n\n\
1. Feladat\n\
Szamoljuk ki a Lagrange interpolacios polinomot\n\
f(x) = x^2-hez az x=1, 2, 3 meresi pontokkal")

disp("\n\
A Vandermonde matrix:")
A = [1 1 1; 1 2 4; 1 3 9]

disp("\n\
A meresi eredmenyek:")
b = [1; 4; 9]

disp("\n\
A megoldas, azaz az interpolacios polinom egyutthatoi:")
egyutthatok = A \ b

disp("\n\
Az eredmeny: 0 + 0 * x + 1 * x^2\n\
Latjuk, hogy visszakaptuk az eredeti x^2-et, ez nem meglepo, hiszen\n\
egy 3 ponton atmeno masodfoku polinomot kerestunk es ilyenbol csak egy letezik.")

###########################################################################

disp("\n\n\
2. Feladat\n\
Szamoljuk ki a Lagrange interpolacios polinomot\n\
f(x) = 3x^2 - 5x + 4-hez az x=1, 2, 3 meresi pontokkal")

f = @(x) 3*x.^2 - 5*x + 4

disp("\n\
A fuggveny ertekei a meresi pontokon:")
f(1)
f(2)
f(3)

disp("\n\
Tehat a meresi eredmenyek:")
b = [f(1); f(2); f(3)]

disp("\n\
Mivel a meresi pontok megegyeznek az elozo feladateval, igy a\n\
Vandermonde matrix ugyanaz")
A

disp("\n\
A megoldas:")
egyutthatok = A \ b

disp("\n\
Az eredmeny: 4 - 5x + 3x^2\n\
Ami megegyezik a kiindulasi fuggvennyel, az indoklas erre ugyanaz mint fentebb.")

###########################################################################

disp("\n\n\
3. Feladat\n\
Szamoljuk ki es abrazoljuk [-2, 6]-on a Lagrange interpolacios polinomot\n\
f(x) = x^3-hoz az x=1, 2, 3 meresi pontokkal")

disp("\n\
Mivel a meresi pontok megegyeznek az elozo feladateval, igy a\n\
Vandermonde matrix ugyanaz")
A

f = @(x) x.^3

disp("\n\
A fuggveny ertekei a meresi pontokon:")
f(1)
f(2)
f(3)

disp("\n\
Tehat a meresi eredmenyek:")
b = [f(1); f(2); f(3)]

disp("\n\
A megoldas:")
egyutthatok = A \ b

disp("\n\
Az eredmeny: 6x^2 - 11x + 6")
p = @(x) 6 * x.^2 - 11 * x + 6


figure (3)

x = -2 : 0.01 : 6;
plot(x, [f(x); p(x)])

# cimet adunk a kepnek
title("Lagrange interpolacio x^3-hoz")

# tovabbi dolgokat rajzolunk ugyanerre a kepre
hold on

# meresi pontok kirajozlasa: 'ks' jelentese k <-> fekete (blac_k_), s <-> negyzet (_s_quare)
plot([1 2 3], [1 8 27], 'ks')

# legend kiirasa ez azonositja, hogy mi micsoda a kepen
legend("x^3", "interpolacio", "meresek")

# befejeztuk a munkat ezzel a keppel
hold off

###########################################################################

disp("\n\n\
4. Feladat\n\
Szamoljuk ki es abrazoljuk egyutt az elozo eredmennyel [-2, 6]-on a\n\
Lagrange interpolacios polinomot\n\
f(x) = x^3-hoz az x=1/2, 0, -1 meresi pontokkal")

disp("\n\
A Vandermonde matrix valtozik, hiszen masok a meresi pontok:")
A = [1 1/2 1/4; 1 0 0; 1 -1 1]

disp("\n\
A meresi ertekek:")
b = [f(1/2); f(0); f(-1)]

disp("\n\
A megoldas:")
egyutthatok = A \ b

disp("\n\
Igy az interpolacios polinom:")
p4 = @(x) -1 / 2 * x.^2 + 1 / 2 * x

# hasznaljuk az elozoleg kiszamolt p(x) polinomot a 3-as feladatbol
figure (4)

plot(x, [f(x); p(x); p4(x)])

# cimet adunk a kepnek
title("Lagrange interpolacio x^3-hoz")

# tovabbi dolgokat rajzolunk ugyanerre a kepre
hold on

# 3. feladat meresi pontjainak kirajozlasa
plot([1 2 3], [1 8 27], 'ks')

# 4. feladat meresi pontjainak kirajozlasa: b* jelentese: b <-> kek (_b_lue), * <-> csillag
plot([1/2 0 -1], [f(1/2) f(0) f(-1)], 'b*')

# legend kiirasa ez azonositja, hogy mi micsoda a kepen
legend("x^3", "[0 1 2] interpolacio", "[1/2 0 -1] interpolacio", "meresek (3)", "meresek (4)")

# befejeztuk a munkat ezzel a keppel
hold off

disp("\n\
Latszik, hogy a ket interpolacio kulonbozo helyeken simul ra az x^3 gorbere")

###########################################################################

disp("\n\n\
5. Feladat\n\
Szamoljuk ki es abrazoljuk [-2, 6]-on a\n\
Lagrange interpolacios polinomot\n\
f(x) = e^x-hez az x=0.1, 0.9, 1.5 meresi pontokkal")

disp("\n\
Meresi pontok:")
xAdat = [0.1 0.9 1.5]

disp("\n\
A Vandermonde matrix eleg hires ahhoz, hogy legyen sajat parancsa:\n\
vander(xAdat) legyartja az x adatsorbol a negyzetes Vandermonde matrixot")
A = vander(xAdat)

# FONTOS, OCTAVE AZ ELOZOEKHEZ KEPEST FORDITVA KEPEZI A MATRIXOT
disp("\n\
Figyeljuk meg, hogy a matrixot forditott oszlopsorrenddel kapjuk")

disp("\n\
Igy a fentebb targyalt linearis egyenletrendszer igy modosul")

A = ["x1^n-1 " " ... " "x1^3 " "x1^2 " "x1 " "1 ";
"x2^n-1 " " ... " "x2^3 " "x2^2 " "x2 " "1 ";
"1 " "... " "... " "... " " ... " "...";
"xn^n-1 " " ... " "xn^3 " "xn^2 " "xn " "1 "]

b = ["y1"; "y2"; "..."; "yn"]

x = ["an-1"; "an-2"; "..."; "a0"]

# vissza a szamokhoz
disp("\n\
Vissza a feladathoz:")
A = vander(xAdat)

# e^x formulaja exp(x)     '_exp_onencialis fuggveny'
disp("\n\
A meresi ertekek:")
b = [exp(0.1); exp(0.9); exp(1.5)]

disp("\n\
A megoldas:")
egyutthatok = A \ b

disp("\n\
Az ebben a formaban kapott egyutthato vektorbol polinomot kepezhetunk es\n\
kiertekelhetjuk a polyval parancs segitsegevel\n\n\
FIGYELEM: az egyutthato vektor jelentese kotott:\n\
- az elso eleme a legnagyobb hatvany kitevoje\n\
- az utolso eleme a konstans tag\n\
azaz eh vektor = [an-1; an-2; ...; a1; a0]")

disp("\n\
Interpolacio 0-ban:")
polyval(egyutthatok, 0)

disp("\n\
Interpolacio 1-ben:")
polyval(egyutthatok, 1)

figure (5)

x = -2 : 0.01 : 6;
plot(x, [exp(x); polyval(egyutthatok, x)])

hold on

yAdat = exp(xAdat)
plot(xAdat, yAdat, 'k*')

title("e^x interpolacioja")

legend("e^x", "interpolacio", "meresek")

hold off

###########################################################################

disp("\n\
A fentebbi modszer problemaja, hogy a Vandermonde matrixok tipikusan rosszul\n\
kondicionaltak.")

A

disp("\n\
Kondicios szama:")
cond(A)

disp("\n\
Ez 3x3-as Vandermonde matrix volt a 16 bar nem nagyon nagy de\n\
a matrix meretehez kepest jelentos")

disp("\n\
vander([1 2 3 4 5 6])")
vander([1 2 3 4 5 6])

disp("\n\
Kondicios szama:")
cond(ans)

disp("\n\
Ez mar objektiven sok.")

disp("\n\
Emlekezzunk vissza, ha Ax=b megoldasanal A kondicios szama nagy, akkor\n\
kis valtozasok b-ben nagy valtozasokat eredmenyeznek x-ben.\n\
Azaz a modell nem megbizhato es a kapott x ertekeket\n\
fenntartasokkal kell kezelni.\n\n\
Igy ennek a modszernek a jelentosege inkabb elmeleti, gyakorlatban\n\
az elozo oran elkeszult lix.m es lagrange.m fuggvenyeket hasznaljuk,\n\
melyek kiszamoljak pontosan a Lagrange interpolacios polinom erteket\n\
es nem az egyutthatokat hatarozzak meg.")

# 5. feladat abraja a lagrange.m es lix.m felhasznalasaval
figure (6)

xAdat = [0.1 0.9 1.5]
yAdat = exp(xAdat)

x = -2 : 0.01 : 6;

plot(x, [exp(x); lagrange(xAdat, yAdat, x)]) #a lagrange.m belul hasznalja a lix.m-et

hold on

plot(xAdat, yAdat, 'k*')

title("e^x interpolacioja lagrange.m hasznalataval")

legend("e^x", "interpolacio", "meresek")

hold off

###########################################################################

disp("\n\n\
6. Feladat\n\
A viz fajhojerol a kovetkezot tudjuk:\n\
- 22  fokon 4181 (J / kg K)\n\
- 42  fokon 4179 (J / kg K)\n\
- 52  fokon 4186 (J / kg K)\n\
- 82  fokon 4199 (J / kg K)\n\
- 100 fokon 4217 (J / kg K)\n\n\
Becsuljuk meg a viz fajhojet 61 fokon kobos spline interpolacioval")
xAdat = [22 42 52 82 100]
yAdat = [4181 4179 4186 4199 4217]

disp("\n\
Az Octave spline(xAdat, yAdat) parancsa kiszamitja a kobos spline interpolaciot.")
spline(xAdat, yAdat)

disp("\n\
az eredmeny formatuma ugynevezett pp ~ _p_iecewise _p_olynomial\n\
azaz szakaszonkent polinom. Termeszetes, hiszen a kobos spline\n\
interpolacio meresi pontok kozott kulonbozo 3-adfoku polinomokkal\n\
kozeliti az ismeretlen fuggvenyt")

disp("\n\
Adjunk az eredmenynek nevet, hogy tudjunk vele dolgozni:")
kobosSpline = spline(xAdat, yAdat)

disp("\n\
a pp strukturat a ppval fuggvennyel tudjuk kiertekelni")

disp("\n\
A fajho becsult erteke 61 fokon = ppval(kobosSpline, 61):")
ppval(kobosSpline, 61)

disp("\n\
Abrazoljuk a kobos spline es a 4-ed foku Lagrange interpolaciot:")

figure (7)

x = 10 : 0.1 : 110;

plot(x, [lagrange(xAdat, yAdat, x); ppval(kobosSpline, x)])

hold on

plot(xAdat, yAdat, 'ks')

title("Viz fajhoje")

legend("4-ed foku Lagrange", "Kobos spline", "Meresek")

hold off