Rešitve nalog, 1. teden

1. naloga

(a) Za vsak \(x_i\) imamo enačbo \(p(x_i) = f(x_i)\) oziroma \(a_n x_i^n + a_{n-1} x_i^{n-1} + \cdots + a_0 = f(x_i)\). V sistemu \(A \mathbf{x} = \mathbf{b}\) je torej \(\mathbf{b} = [f(x_0), \ldots, f(x_k)]^\mathsf{T}\), \[ A = \begin{bmatrix} x_0^n & x_0^{n-1} & \cdots & x_0 & 1 \\x_1^n & x_1^{n-1} & \cdots & x_1 & 1 \\ x_2^n & x_2^{n-1} & \cdots & x_2 & 1 \\ \vdots & \vdots & & \vdots & \vdots \\ x_k^n & x_k^{n-1} & \cdots & x_k & 1\end{bmatrix},\] vektor neznank pa \(\mathbf{x} = [a_n, a_{n-1}, \ldots, a_1, a_0]^\mathsf{T}\).

(b) Polinom stopnje \(0\) je oblike \(p(x) = a_0\), pripadajoča matrika \(A\) je torej \[ A = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}, \] desna stran pa \(\mathbf{b} = [f(-1), f(0), f(1)]^\mathsf{T} = [\frac{1}{2},0,\frac{1}{2}]^\mathsf{T}\). Rešitev po linearni metodi najmanjših kvadratov lahko poiščemo tako, da poiščemo rešitev pripadajočega normalnega sistema \[ A^\mathsf{T} A \mathbf{x} = A^\mathsf{T} \mathbf{b}. \] V tem primeru je \(A^\mathsf{T} A = 3\), \(A^\mathsf{T} \mathbf{b} = 1\), torej \(3a_0 = 1\) ali \(p(x) = a_0 = \frac{1}{3}\). Do rešitev za aproksimaciji s polinomom stopnje \(1\) oziroma \(2\) pridemo analogno.

(c) Pripravimo si podatke. V vrstico xp shranimo 21 ekvidistantnih točk na intervalu \([-1,1]\):

octave:1> xp = linspace(-1, 1, 21);

V teh točkah izračunamo funkcijske vrednosti \(f(x) = \frac{x^2}{1+x^2}\):

octave:2> yp = xp.^2./(1+xp.^2);

Poiščimo najprej kar aproksimacijo s polinomom stopnje 2. Pripravimo matriko sistema \(A\):

octave:3> A = [xp'.^2, xp', ones(21, 1)];

in desno stran \(\mathbf{b}\):

octave:4> b = yp';

Koeficienti polinoma so

octave:5> a = A\b;

Narišimo še pripadajoče grafe. V ta namen si pripravimo nabor 100 točk na intervalu \([-1,1]\):

octave:6> x = linspace(-1, 1);

V teh točkah poračunamo funkcijske vrednosti \(f\):

octave:7> y = x.^2./(1+x.^2);

in vrednosti polinomske aproksimacije:

octave:8> ya = polyval(a, x);

Na isto sliko narišemo graf originalne funkcije in graf polinoma ter označimo položaje točk, ki smo jih uporabili za izračun aproksimacije:

octave:9> plot(x, y)
octave:10> hold on
octave:12> plot(x, ya)
octave:11> plot(xp, yp, 'ro')

Poiščimo še polinom stopnje 20, ki se najbolje prilega danim 21 točkam. To bo v resnici kar interpolacijski polinom, ki se točno prilega danim točkam. Pripravimo novo matriko \(A\):

octave:12> A = xp'.^(20:-1:0); % a:t:b vrne vrstico stevil od a do b s korakom t

Desna stran \(\mathbf{b}\) je enaka kot prej, koeficienti interpolacijskega polinoma pa so:

octave:13> a = A\b;

Pripravimo še dovolj gost nabor vrednosti tega polinoma in ga narišimo na isto sliko:

octave:14> ya = polyval(a, x);
octave:15> plot(x, ya)

Kako se polinomske aproksimacije in interpolacije obnašajo z realnimi podatki? Pripravimo najprej podatke z umetno dodano napako:

octave:16> yp = xp.^2./(1+xp.^2) + 0.001*rand(1, 21);

Uporabili smo funkcijo rand; klic rand(m, n) vrne \(m \times n\) matriko naključnih vrednosti med \(0\) in \(1\). Sedaj lahko na teh podatkih ponovimo klice od vrstice 3 naprej. Naredi to in primerjaj rezultate. 

Ponovi nalogo še s funkcijo \(g(x) = \frac{1}{1+25x^2}\).

2. naloga

(a) Če odštejemo dve zaporedni enačbi oblike \((x-p_i)^2 + (y-q_i)^2 = d_i^2\) dobimo sistem enačb \[ 2(p_i - p_{i+1}) x + 2(q_i - q_{i+1})y = d_i^2-d_{i+1}^2-p_i^2+p_{i+1}^2-q_i^2+q_{i+1}^2\] za vse \(i = 1,2,\ldots,n-1\).

(b) V matriko \(A\) shranimo koeficiente, ki nastopajo pri neznankah \(x\) in \(y\) v zgornjem sistemu enačb: \[ A = 2 \begin{bmatrix} p_1- p_{2} & q_1 - q_{2} \\ p_2- p_{3} & q_2 - q_{3}  \\ \vdots & \vdots \\ p_{n-1}- p_{n} & q_{n-1} - q_{n}\end{bmatrix},\] v vektor \(\mathbf{b}\) pa koeficinte na desni strani: \[ \mathbf{b} = \begin{bmatrix} d_1^2-d_{2}^2-p_1^2+p_{2}^2-q_1^2+q_{2}^2 \\ d_2^2-d_{3}^2-p_2^2+p_{3}^2-q_2^2+q_{3}^2 \\ \vdots \\ d_{n-1}^2-d_{n}^2-p_{n-1}^2+p_{n}^2-q_{n-1}^2+q_{n}^2 \end{bmatrix}.\]

(c) Funkcija sprejemnik. 

(d) Datoteki oddajniki.txt in razdalje.txt shranimo v mapo, v kateri poganjamo octave. Podatke shranimo:

octave:1> P = load("oddajniki.txt")
octave:2> d = load("razdalje.txt")

in poiščemo (npr.) položaj prvega sprejemnika

octave:3> sprejemnik(P, d(:, 1))