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

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

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)\n\
Kerdes:\n\
- milyen becslest tudunk adni f(x)-re, ahol x nem xi egyik i-re sem\n\n\
Ha x < min(x1, x2, ..., xn) vagy x > max(x1, x2, ..., xn),\n\
akkor extrapolaciorol beszelunk, ellenkezo esetben interpolaciorol.\n\n\
A most targyalando polinomos kozelitesek nem alkalmasak extrapolaciora,\n\
a meresunk altal lefedett x tartomanyon kivul gyakorlatilag\n\
hasznalhatatlan eredmenyeket kapunk.")

disp("\n\
Egy k-ad foku polinom x-ben: a0 + a1 x + a2 x^2 + ... + ak x^k\n\
(k + 1) db egyutthatoval rendelkezik: a0, a1, ..., ak\n\n\
Ha tehat szeretnenk a meresi adatunkat egy k-ad foku polinommal modellezni,\n\
akkor (k + 1) egyenletre lenne szuksegunk. Nekunk n db van: yi = f(xi),\n\
tehat legfeljebb (n - 1)-ed foku polinommal tudjuk modellezni a jelenseget.")

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\
Lagrange modszere:\n\
- Segedpolinomokat konstrualunk: l_i(x)\n\
  olyan (n-1)-ed foku polinom ami x = xi eseten 1, x = xj eseten 0\n\
                                                   ha i nem = j")
                                                   
disp("\n\
Konnyu latni, hogy\n\
l_i(x) = (x - x_1)/(x_i - x_1) * ... * (x - x_i-1)/(x_i - x_i-1) * \n\
         (x - x_i+1)/(x_i - x_i+1) * ... * (x - x_n)/(x_i - x_n)")

disp("\n\
Ha kiszamitottuk az osszes l_i(x) polinomot, a Lagrange interpolacios\n\
polinom p(x) = y1 * l_1(x) + ... + yn * l_n(x)\n\
formaban trivialisan adodik, hiszen\n\
p(xi) = y1 * 0 + ... yi-1 * 0 + yi * 1 + yi+1 * 0 + ... + yn * 0 = yi")
         
disp("\n\
Toltsuk le a lix.m es lagrange.m fileokat!\n\
- lix(i, xAdat, x) kiszamlja l_i(x)-et\n\
  az xAdat-ban megadott meresi pontokra\n\
- lagrange(xAdat, yAdat, x) kiszamolja az interpolacios polinomot x-ben")

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

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

xAdat = [1 2 3]
yAdat = xAdat .^2

figure (1);
x = -4 : 0.01 : 8;
plot(x, [x .^2; lagrange(xAdat, yAdat, x)])
title("x^2 es a Lagrange interpolacios kozelites 3 meresi ponttal")

disp("Nyilvan 3 pontra egy (3-1)=2 foku polinomot illesztunk,\n\
s mivel az eredeti modellunk is 2 foku, igy a kozelites megegyezik az\n\
eredeti fuggvennyel.")

disp("\n\
Kezi szamitashoz:\n\
- l_1(x) = (x - x2) (x - x3) / (x1 - x2) (x1 - x3) =\n\
           (x^2 - 5x + 6) / (-1 * -2) = (x^2 - 5x + 6) / 2\n\
- l_2(x) = (x - x1) (x - x3) / (x2 - x1) (x2 - x3) =\n\
           (x^2 - 4x + 3) / (1 * -1) = (x^2 - 4x + 3) / -1\n\
- l_3(x) = (x - x1) (x - x2) / (x3 - x1) (x3 - x2) =\n\
           (x^2 - 3x + 2) / (1 * 2) = (x^2 - 3x + 2) / 2\n\n\
Tehat p(x) = 1 * l_1(x) + 4 * l_2(x) + 9 * l_3(x) = \n\
             [1 * (x^2 - 5x + 6) - 4 * 2(x^2 - 4x + 3) + 9 * (x^2 - 3x + 2)] / 2\n\
             (2 x^2 + 0 * x + 0) / 2 = x^2")
###########################################################################

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

xAdat = [1 2 3]
yAdat = xAdat .^3

figure (2);
x = -4 : 0.01 : 8;
plot(x, [x .^3; lagrange(xAdat, yAdat, x)])
title("x^3 es a Lagrange interpolacios kozelites 3 meresi ponttal")

disp("Nyilvan 3 pontra egy (3-1)=2 foku polinomot illesztunk,\n\
s mivel az eredeti modellunk 3 foku igy ertheto a nagy elteres\n\
az interpolacios intervallumon [1, 3] kivul.")

disp("\n\
Kezi szamitashoz:\n\
- l_1(x), l_2(x), l_3(x) ugyanaz mint elobb, hisz ugyanazok az x meresi pontok\n\n\
Tehat p(x) = 1 * l_1(x) + 8 * l_2(x) + 27 * l_3(x) = \n\
             [1 * (x^2 - 5x + 6) - 8 * 2(x^2 - 4x + 3) + 27 * (x^2 - 3x + 2)] / 2\n\
             (12 x^2 - 22 * x + 12) / 2 = 6 x^2 - 11 x + 6")
             
figure (3);
x = -4 : 0.01 : 8;
plot(x, [x .^3; 6 * x.^2 - 11 * x + 6])
title("x^3 es a Lagrange interpolacios kozelites 3 meresi ponttal")

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

disp("\n\n\
3. 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\n\
- a) 1. foku Lagrange interpolacioval\n\
- b) 2. foku Lagrange interpolacioval\n\
- c) 3. foku Lagrange interpolacioval\n\
- d) 4. foku Lagrange interpolacioval")

disp("\n\n\
Osszesitjuk majd a vegen az eredmenyt, igy minden interpolaciot abrazolni fogunk")
x = 10 : 0.1 : 110;

disp("\n\n\
a) Elsofoku interpolaciohoz ket meresi pontot hasznalhatunk fel ( 2 - 1 = 1 )\n\
   Mivel 61 fokon keressuk a fajhot, igy celszeru a ket legkozelebbi meresi\n\
   pontot hasznalni, azaz az 52 es 82 fokon mert adatokat.")
xAdat = [52 82]
yAdat = [4186 4199]
lagrange(xAdat, yAdat, 61)

xLagrangeElsofoku = lagrange(xAdat, yAdat, x);

disp("\n\n\
b) Masodfoku interpolaciohoz harom meresi pontot hasznalhatunk fel ( 3 - 1 = 2 )\n\
   Mivel 61 fokon keressuk a fajhot, igy celszeru a harom legkozelebbi meresi\n\
   pontot hasznalni, azaz a 42, 52 es 82 fokon mert adatokat.")
xAdat = [42 52 82]
yAdat = [4179 4186 4199]
lagrange(xAdat, yAdat, 61)

xLagrangeMasodfoku = lagrange(xAdat, yAdat, x);

disp("\n\n\
c) Harmadfoku interpolaciohoz negy meresi pontot hasznalhatunk fel ( 4 - 1 = 3 )\n\
   Mivel 61 fokon keressuk a fajhot, igy celszeru a negy legkozelebbi meresi\n\
   pontot hasznalni, azaz a 22, 42, 52 es 82 fokon mert adatokat\n\
   (22 helyett a 100 is megfelel, hisz 61 - 22 = 100 - 61 = 39).")
xAdat = [22 42 52 82]
yAdat = [4181 4179 4186 4199]
lagrange(xAdat, yAdat, 61)

xLagrangeHarmadfoku = lagrange(xAdat, yAdat, x);

disp("\n\n\
d) Negyedfoku interpolaciohoz ot meresi pontot hasznalhatunk fel ( 5 - 1 = 4 )\n\
   Igy mindegyik meresi adatot hasznaljuk.")
xAdat = [22 42 52 82 100]
yAdat = [4181 4179 4186 4199 4217]
lagrange(xAdat, yAdat, 61)

xLagrangeNegyedfoku = lagrange(xAdat, yAdat, x);

figure (4);
plot(x, [xLagrangeElsofoku; xLagrangeMasodfoku; xLagrangeHarmadfoku; xLagrangeNegyedfoku]);
title("Kulonbozo foku Lagrange interpolaciok")
legend("1. foku", "2. foku", "3. foku", "4. foku");
grid on;

# kirajzoljuk a meresi pontokat is
hold on;

plot(xAdat, yAdat, 'ks'); #a ks opciobol k: fekete s: negyzetek

hold off;