interpolacija

library(tidyverse)

Interpolacijski polinom

Primjer 1. Odredite interpolacijski polinom kroz točke \((-2,0)\), \((-1,1)\), \((1,2)\) i \((2,0)\).

podaci <- tibble(x = c(-2, -1, 1, 2), 
                 y = c(0, 1, 2, 0))

Pomoću naredbe lm možemo dobiti koeficijente interpolacijskog polinoma.

model1 <- lm(y ~ 1 + x + I(x^2) + I(x^3), data = podaci)
coef(model1)

## (Intercept)           x      I(x^2)      I(x^3) 
##   2.0000000   0.6666667  -0.5000000  -0.1666667

U općenitijem slučaju je jednostavnije unijeti formulu pomoću naredbe poly. Opcija raw = TRUE daje koeficijente s obzirom na kanonsku bazu, a raw = FALSE daje koeficijente s obzirom na pripadnu ortogonalnu bazu.

model2 <- lm(y ~ poly(x, 3, raw = TRUE), data = podaci)
coef(model2)

##             (Intercept) poly(x, 3, raw = TRUE)1 poly(x, 3, raw = TRUE)2 
##               2.0000000               0.6666667              -0.5000000 
## poly(x, 3, raw = TRUE)3 
##              -0.1666667

Ako tražimo polinom čiji stupanj je veći lili jednak od broja točaka, onda uz veće potencije koeficijenti nisu određeni jer je premalo zadanih točaka.

model3 <- lm(y ~ poly(x, 10, raw = TRUE), data = podaci)
coef(model3)

##               (Intercept)  poly(x, 10, raw = TRUE)1  poly(x, 10, raw = TRUE)2 
##                 2.0000000                 0.6666667                -0.5000000 
##  poly(x, 10, raw = TRUE)3  poly(x, 10, raw = TRUE)4  poly(x, 10, raw = TRUE)5 
##                -0.1666667                        NA                        NA 
##  poly(x, 10, raw = TRUE)6  poly(x, 10, raw = TRUE)7  poly(x, 10, raw = TRUE)8 
##                        NA                        NA                        NA 
##  poly(x, 10, raw = TRUE)9 poly(x, 10, raw = TRUE)10 
##                        NA                        NA

Graf interpolacijskog polinoma

ggplot(podaci, aes(x = x, y = y)) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 3, raw = TRUE), se = FALSE, fullrange = TRUE) +
  geom_point(size = 2.5, color = "red", shape = 21, fill = "yellow") +
  ggtitle("interpolacijski polinom") + xlim(-4, 2.5)

Predikcije pomoću dobivenog polinoma

Vrijednosti u kojima želimo napraviti predikcije.

vr <- c(-3, -1.5, -0.5, 0.5, 1)

Predikcije preko model1

pr1 <- predict(model1, tibble(x = vr))
names(pr1) <- vr
pr1

##            -3          -1.5          -0.5           0.5             1 
## -4.440892e-16  4.375000e-01  1.562500e+00  2.187500e+00  2.000000e+00

Predikcije preko model2

pr2 <- predict(model2, tibble(x = vr))
names(pr2) <- vr
pr2

##            -3          -1.5          -0.5           0.5             1 
## -4.440892e-16  4.375000e-01  1.562500e+00  2.187500e+00  2.000000e+00

Predikcije preko model3. Upozorava nas na neodređene koeficijente uz veće potencije zbog premalo točaka, no predikcije dobro funkcioniraju.

pr3 <- predict(model3, tibble(x = vr))

## Warning in predict.lm(model3, tibble(x = vr)): prediction from a rank-deficient
## fit may be misleading

names(pr3) <- vr
pr3

##            -3          -1.5          -0.5           0.5             1 
## -4.440892e-16  4.375000e-01  1.562500e+00  2.187500e+00  2.000000e+00

Interpolacija polinomom velikog stupnja

Promatramo funkciju \(f(x)=\dfrac{1}{1+x^2}\) na segmentu \([-5,5]\).

zvono <- as_mapper(~1 / (1 + .x^2))

Definiramo ekvidistantne razdiobe segmenta \([-5,5]\) s \(5\), \(10\), \(15\) i \(20\) točaka.

razdioba <- list()
for (i in c(5,10,15,20)) {
  razdioba[[as.character(i)]] <- list(x = seq(from = -5, to = 5, length.out = i),
                                   y = zvono(seq(from = -5, to = 5, length.out = i))
                                   )
}
razdioba <- razdioba %>% map(as_tibble) %>% map_dfr(rbind, .id = "broj") %>%
  mutate(broj = fct_relevel(broj, "5"))

Vidimo da interpolacijski polinom većeg stupnja loše aproksimira funkciju \(f\) blizu rubova segmenta \([-5,5]\).

ggplot(razdioba, aes(x = x, y = y)) +
  stat_function(fun = zvono, color = "red") +
  geom_smooth(method = "lm", formula = y ~ poly(x, 19, raw = TRUE), 
              se = FALSE, fullrange = TRUE, linewidth = 0.7) +
  geom_point(size = 2, color = "red", shape = 21, fill = "yellow") +
  facet_wrap(vars(broj), scales = "free_y")

Čebiševljevi polinomi

Čebiševljev polinom stupnja \(n\) je funkcija \[T_n:[-1,1]\to\mathbb{R},\quad T_n(x)=\cos{(n\arccos{x})},\quad n=0,1,2,\dotsc\] Nultočke od \(T_n\) su dane s \(x_k=\cos{\dfrac{2k+1}{2n}\pi},\quad k=0,1,2,\dotsc,n-1\).

Ako interpoliramo neku funkciju na segmentu \([-1,1]\), za čvorove interpolacije možemo uzeti i nultočke nekog Čebiševljevog polinoma (ne moraju čvorovi uvijek biti dobiveni podjelom segmenta na \(n\) jednakih dijelova, nego to mogu biti bilo koje točke unutar zadanog segmenta). Međutim, ako imamo proizvoljni segment \([a,b]\), tada nultočke Čebiševljevog polinoma možemo transformirati na taj segment pomoću linearne transformacije \[\varphi:[-1,1]\to[a,b],\qquad \varphi(x)=\frac{b-a}{2}x+\frac{b+a}{2}\] i za čvorove interpolacije uzeti točke \[\varphi(x_k),\quad k=0,1,\dotsc,n-1\] gdje su \(x_k\) nultočke Čebiševljevog polinoma \(T_n\).

Definiramo Čebiševljeve razdiobe segmenta \([-5,5]\) s \(5\), \(10\), \(15\) i \(20\) točaka.

cheb_nul <- list()
for (i in c(5,10,15,20)) {
  cheb_nul[[as.character(i)]] <- cos((2*(0:(i-1)) + 1) * pi / (2 * i))
}
a <- -5
b <- 5
cheb_razdioba <- map(cheb_nul, ~(b - a) / 2 * .x + (b + a) / 2)
cheb_razdioba <- cheb_razdioba %>% map(as_tibble) %>% 
  map_dfr(rbind, .id = "broj") %>%  rename(x = value) %>%
  mutate(broj = fct_relevel(broj, "5"), y = zvono(x))

U ovom slučaju interpolacijski polinom većeg stupnja dobro aproksimira funkciju \(f\) na segmentu \([-5,5]\).

ggplot(cheb_razdioba, aes(x = x, y = y)) +
  stat_function(fun = zvono, color = "red") +
  geom_smooth(method = "lm", formula = y ~ poly(x, 19, raw = TRUE), 
              se = FALSE, fullrange = TRUE, linewidth = 0.7) +
  geom_point(size = 2, color = "red", shape = 21, fill = "yellow") +
  facet_wrap(vars(broj), scales = "free_y")

Kubični splajn

Napravimo ekvidistantnu razdiobu segmenta \([-5,5]\) s \(10\) točaka i pripadnim vrijednostima funkcije \(f(x)=\dfrac{1}{1+x^2}\) u tim točkama.

tocke <- tibble(x = seq(from = -5, to = 5, length = 10), 
                y = zvono(x))
tocke

## # A tibble: 10 × 2
##         x      y
##     <dbl>  <dbl>
##  1 -5     0.0385
##  2 -3.89  0.0620
##  3 -2.78  0.115 
##  4 -1.67  0.265 
##  5 -0.556 0.764 
##  6  0.556 0.764 
##  7  1.67  0.265 
##  8  2.78  0.115 
##  9  3.89  0.0620
## 10  5     0.0385

Tražimo prirodni kubični splajn na danim podacima tako da generiramo \(100\) vrijednosti tog splajna na segmentu \([-5,5]\). Prirodni kubični splajn u rubnim točkama segmenta ima druge derivacije jednake nula.

natur <- spline(tocke$x, tocke$y, n = 100, method = "natural") %>% as_tibble()
natur

## # A tibble: 100 × 2
##        x      y
##    <dbl>  <dbl>
##  1 -5    0.0385
##  2 -4.90 0.0398
##  3 -4.80 0.0411
##  4 -4.70 0.0425
##  5 -4.60 0.0441
##  6 -4.49 0.0458
##  7 -4.39 0.0477
##  8 -4.29 0.0499
##  9 -4.19 0.0524
## 10 -4.09 0.0552
## # … with 90 more rows

Graf prirodnog kubičnog splajna i funkcije \(f(x)=\frac{1}{1+x^2}\)

ggplot(natur, aes(x, y)) +
  stat_function(fun = zvono, n = 180, aes(color = "zvono")) +
  geom_line(aes(color = "prirodni splajn")) +
  geom_point(data = tocke, size = 2, color = "red", shape = 21, fill = "yellow") +
  scale_color_manual("Graf", values = c("zvono" = "red", "prirodni splajn" = "blue"))

fmm metoda daje kubični splajn koji u rubnim točkama segmenta ima treće derivacije jednake trećim derivacijama kubičnih polinoma koje dobijemo tako da tražimo interpolacijske polinome kroz prve četiri točke razdiobe segmenta i kroz zadnje četiri točke razdiobe segmenta.

fmm <- spline(tocke$x, tocke$y, n = 100, method = "fmm") %>% as_tibble()
fmm

## # A tibble: 100 × 2
##        x      y
##    <dbl>  <dbl>
##  1 -5    0.0385
##  2 -4.90 0.0400
##  3 -4.80 0.0415
##  4 -4.70 0.0430
##  5 -4.60 0.0446
##  6 -4.49 0.0463
##  7 -4.39 0.0482
##  8 -4.29 0.0503
##  9 -4.19 0.0527
## 10 -4.09 0.0554
## # … with 90 more rows

Graf fmm kubičnog splajna i funkcije \(f(x)=\frac{1}{1+x^2}\)

ggplot(fmm, aes(x, y)) +
  stat_function(fun = zvono, n = 180, aes(color = "zvono")) +
  geom_line(aes(color = "fmm splajn")) +
  geom_point(data = tocke, size = 2, color = "red", shape = 21, fill = "yellow") +
  scale_color_manual("Graf", values = c("zvono" = "red", "fmm splajn" = "blue"))

Iako sa slika ne vidimo veliku razliku između prirodnog i fmm splajna, ipak to nisu iste funkcije ako pogledamo njihove vrijednosti u gornjim tablicama. Možemo definirati splajnove kao funkcije i računati vrijednosti u željenim točkama.

spl1 <- splinefun(tocke$x, tocke$y, method = "natural")
spl2 <- splinefun(tocke$x, tocke$y, method = "fmm")

Možemo vidjeti da na primjer u točki \(4.9\) ova dva splajna nemaju iste vrijednosti. Razlika tih vrijednosti je jako mala pa ih na slici nije lako uočiti osim ako jako zumiramo sliku.

c(spl1(4.9), spl2(4.9))

## [1] 0.03975041 0.03996341

Zumirana slika

ggplot() +
  geom_line(data = fmm, aes(x, y, color = "fmm splajn")) +
  geom_line(data = natur, aes(x, y, color = "prirodni splajn")) +
  scale_color_manual("Graf", values = c("prirodni splajn" = "red", "fmm splajn" = "blue")) +
  xlim(4.19,5) + ylim(0.035,0.055)

Što je veći broj točaka u razdiobi segmenta \([-5,5]\), kubični splajn bolje aproksimira funkciju \(f(x)=\dfrac{1}{1+x^2}\). Pogledajmo to na primjeru prirodnog splajna.

splajn <- razdioba %>% group_split(broj) %>% map(~splinefun(.$x, .$y, method = "natural"))
names(splajn) <- levels(razdioba$broj)

x_spl <- rep(with(razdioba, seq(min(x), max(x), length = 100)), times = length(splajn))
broj_spl <- rep(names(splajn), each = 100)
y_spl <- map2_dbl(x_spl, broj_spl, ~splajn[[.y]](.x))

tab_spl <- tibble(broj = fct_relevel(broj_spl, "5"), x = x_spl, y = y_spl)

ggplot(razdioba, aes(x = x, y = y)) +
  stat_function(fun = zvono, color = "red") +
  geom_line(data = tab_spl, aes(x = x, y = y), color = "blue") +
  geom_point(size = 1.5, color = "red", shape = 21, fill = "yellow") +
  facet_wrap(vars(broj), scales = "free_y") +
  ggtitle("Interpolacija s prirodnim splajnom")

---

Primjer 2. Odredite prirodni kubični splajn za podatke dane u donjoj tablici.

x	0	2.0	2.5	3.0	3.5	4.0	4.5	5.0	6
y	0	2.9	3.5	3.8	3.5	3.5	3.5	2.6	0

Zadane točke i pripadni prirodni kubični splajn

tocke1 <- tibble(x = c(0, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6), 
                 y = c(0, 2.9, 3.5, 3.8, 3.5, 3.5, 3.5, 2.6, 0))
natur1 <- spline(tocke1$x, tocke1$y, n = 100, method = "natural") %>% as_tibble()

Graf pripadnog prirodnog kubičnog splajna

ggplot(natur1, aes(x, y)) +
  geom_line(color = "blue") +
  geom_point(data = tocke1, size = 2, color = "red", shape = 21, fill = "yellow")

---

Primjer 3. Odredite prirodni kubični splajn za podatke dane u donjoj tablici.

x	0	1	2	3	4
y	0	0	1	0	0

Zadane točke i pripadni prirodni kubični splajn

tocke2 <- tibble(x = c(0, 1, 2, 3, 4), 
                 y = c(0, 0, 1, 0, 0))
natur2 <- spline(tocke2$x, tocke2$y, n = 100, method = "natural") %>% as_tibble()

Graf pripadnog prirodnog kubičnog splajna

ggplot(natur2, aes(x, y)) +
  geom_line(color = "blue") +
  geom_point(data = tocke2, size = 2, color = "red", shape = 21, fill = "yellow")

Hermite interpolacija

Kubični splajn u čvorovima ima neprekidne derivacije do drugog reda, dok po dijelovima kubična Hermiteova interpolacija u čvorovima ima samo neprekidnu derivaciju prvog reda. Uz zadane čvorove i funkcijske vrijednosti moramo zadati i vrijednosti derivacije u čvorovima.

Pogledajmo sve na primjeru funkcije \(f(x)=\dfrac{1}{1+x^2}\) na segmentu \([-5,5]\) sa \(10\) ekvidistantnih čvorova.

der_zvono <- D(expression(1 / (1 + x^2)), "x")
tocke <- tocke %>% mutate(dy = eval(der_zvono))
herm <- splinefunH(tocke$x, tocke$y, tocke$dy)

herm_tocke <- tibble(x = seq(-5, 5, length = 150),
                     y = herm(x),
                     d1 = herm(x, deriv = 1),
                     d2 = herm(x, deriv = 2))

ggplot(herm_tocke, aes(x, y)) +
  stat_function(fun = zvono, n = 180, aes(color = "zvono")) +
  geom_line(aes(color = "Hermite")) +
  geom_point(data = tocke, size = 2, color = "red", shape = 21, fill = "yellow") +
  scale_color_manual("Graf", values = c("zvono" = "red", "Hermite" = "forestgreen"))

Usporedba prirodnog splajna i Hermiteove interpolacije za navedene čvorove

ggplot() +
  stat_function(fun = zvono, n = 180, aes(color = "zvono")) +
  geom_line(data = herm_tocke, aes(x, y, color = "Hermite")) +
  geom_line(data = natur, aes(x, y, color = "prirodni splajn")) +
  geom_point(data = tocke, aes(x, y), size = 2, color = "red", shape = 21, fill = "yellow") +
  scale_color_manual("Graf", values = c("prirodni splajn" = "blue", 
                                        "Hermite" = "forestgreen",
                                        "zvono" = "red"))

Pogledajmo grafove prirodnog splajne i njegove prve i druge derivacije. Vidimo da su obje derivacije neprekidne funkcije.

natur <- natur %>% mutate(d1 = spl1(x, deriv = 1),
                          d2 = spl1(x, deriv = 2))

ggplot(natur, aes(x, y)) +
  geom_line(aes(color = "prirodni splajn")) +
  geom_line(aes(x, d1, color = "1. derivacija")) +
  geom_line(aes(x, d2, color = "2. derivacija")) +
  scale_color_manual("Graf", values = c("prirodni splajn" = "red", 
                                        "1. derivacija" = "blue",
                                        "2. derivacija" = "forestgreen")) +
  ggtitle("Prirodni splajn i njegove derivacije")

Pogledajmo grafove Hermiteovog “splajna” i njegove prve i druge derivacije. Uočite kako prva derivacija nije glatka u čvorovima pa u tim točkama ne postoji druga derivacija, odnosno druga derivacija više nije neprekidna u čvorovima. Preciznije, druga derivacija u čvorovima se ne može dodefinirati tako da postane neprekidna funkcija. Pazite, vertikalne zelene linije nisu dio grafa druge derivacije, u okolini čvorova druga derivacija ima nagle promjene (skokove).

ggplot(herm_tocke, aes(x, y)) +
  geom_line(aes(color = "Hermite")) +
  geom_line(aes(x, d1, color = "1. derivacija")) +
  geom_line(aes(x, d2, color = "2. derivacija")) +
  scale_color_manual("Graf", values = c("Hermite" = "red", 
                                        "1. derivacija" = "blue",
                                        "2. derivacija" = "forestgreen")) +
  ggtitle("Hermite 'splajn' i njegove derivacije")