-
Notifications
You must be signed in to change notification settings - Fork 57
/
cs2_lab01.Rmd
156 lines (132 loc) · 6.6 KB
/
cs2_lab01.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
---
number: 1
course: Computer Science 2
material: Instruction 1
author: B. Górecki, rev. W. Gryglas, eng. T.Wacławczyk
---
# Interpolation
Files for the present tutorial can be found using following links:
- [Header file inter.h](http://ccfd.github.io/courses/code/info2/inter.h)
- [Source file inter.cpp](http://ccfd.github.io/courses/code/info2/inter.cpp)
## 1. Interpolation in engineering practice
Today, we will deal with the issue of interpolation.
Suppose that from a tank that can be placed at different heights, water flows out of the pipe until it is empty.
We made the experiment by placing the tank at several different heights (eg 2, 3 and 4 m) and for each of these heights we measure the time after which the tank was empty. Note, the pipe outlet remains on the same level.
In the design process (optimization?), we usually want to be able to predict the emptying time for any tank
placement height. This height is between the set extreme values, i.e. at an altitude of 3.27 m.
Therefore, we must be able to use our measurement data to carry out a curve that will reliably approximate this relationship.
We are talking about interpolation, i.e. interpolating a set of discrete data (interpolation nodes) into a continuous relationship defined for each argument lying between the points (interpolation nodes) on which we base interpolation.
### Meritum
Interpolation is the issue of leading a curve through all points from the data set so as to obtain the course of the assumed dependence between the measurement points.
The number of degrees of freedom of the fitted curve must be equal to the number of points in the data set (the number of equations must be equal to the number of unknowns).
## 2. Polynomial interpolation
Let us condider a Lagrange interpolation polynomial of order $n$ defined as:
$$
w_n(x) = \sum_{i=0}^{n} y_i \prod_{j=0 \land j\neq i}^{n} \frac{x - x_j}{x_i - x_j}
$$
where $(x_i, y_i)$ denote the $i-th$ coordinate of single measurement point (interpolation node).
### Exercises
1. Write a program generating n points (interpolation nodes) using following definitions:
$$
(x_i,y_i) \rightarrow \left( x_i, \exp(-x_i^2) \right)
$$
where
$$
x_i = a + i \cdot h,\qquad y_i=\exp(-x_i^2)
$$
$$
h = \frac{b - a}{n-1}
$$
$$
i = 0, \ldots n-1
$$
and $a = -2$ i $b = 2$.
**Repetition:** usage of `malloc` function for dynamic mememory allocation of a 1D array:
```c++
double *x;
x = (double *) malloc(n * sizeof(double));
if (x == NULL) {
fprintf(stderr, "malloc: can not allocate x.\n");
exit(1);
}
/* free memory x */
free(x);
```
2. Compute value of the Lagrange interpolation polynomial in points with cooridinates:
$$
t_i = a + i \cdot \frac{h}{4}, \; i = 0, \ldots 4n - 4.
$$
Using function `lagrange(double *x, double *y, int n, double xx)` compute the value of the interpolation polynomial in the arbitrary point `xx` located in between the interpolation nodes. An example of basic usage of `lagrange()` function is given below:
```c++
int main() {
double x[3] = {0, 1, 2};
double y[3] = {0, 1, 4};
double f, t = 1.5
/* Compute value of the Lagrange polynomial of order 2 at t = 1.5 */
f = lagrange(x, y, 3, t);
}
```
3. Print obtained results (for each $t_i$) on the screen.
4. Display result on the diagram using function `scale(double x0, double y0, double x1, double y1)` defined in `winbgi2` library. Experiment with different `x0,y0,x1,y1` parameters. An example below shows how to use it to display sin(x) function:
```c++
int main() {
double pi = 4. * atan(1.);
graphics(600, 400);
scale(0, -1.2, 7, 1.2); // xmin = 0, ymin = -1.2, x_max = 7, ymax = 1.2
double x = 0;
while(x < 2 * pi) {
point(x, sin(x));
x += 0.01;
}
wait();
}
```
5. Write results to the file.
6. Compare results for different number of interpolation nodes $n$.
7. For a selected $n$ calculate and display maximal error of the interpolation.
8. Repeat your results for different set of interpolation nodes (using different function to create values `y[i]` in interpolation nodes).
Look for the results obtained using function $|x|$.
Can you explain why they look like that (see definition of the intepolation error)?
### Note
Parameters of the function `lagrange(double *x, double *y, int n, double xx)` denote:
- `x`, `y` - pointers to n-element arrays holding coordinates of the interpolation nodes,
- `n` - the size of the vector (number of its elements),
- `xx` - current value of the interpolating polynomial argument for which we compute the value of the Lagrange polynomial.
## 3. Runge effect
Let us check how to improve quality of the interpolation by introduction of a more
clever, non-uniform interpolation nodes distribution. We will use function $|x|$ in this example to generate the set of interpolation nodes $(x_i,|x_i|)$.
### Chebyschev distribution of interpolation nodes
Up to now, interpolation nodes were distributed uniformly in the interpolation domain.
Let us now to choose the distribution in such a way that they are denser close to the
ends of the interpolation domain.
This can be done as follows (see Lecture notes for theoretical explanation).
We want to distribute interpolation nodes in the range $x = [-1, 1]$.
Now, let us plot half-circle with center in the middle of this range $x = 0$ and radius $R = 1$.
The uniform division of an resulting arc (measuring in local arc coordinate) on $n - 1$ fragments
allows to generate abscisas of the interpolation nodes given be the formula:
$$
x_i = -\cos\left( \frac{i \cdot \pi}{n-1} \right), \; i = 0, \ldots n - 1
$$
In the more general case of the range $x = [a, b]$ the following formula is valid:
$$
x_i = -\frac{b-a}{2} \cos\left( \frac{i \cdot \pi}{n-1} \right) + \frac{a+b}{2}, \; i = 0, \ldots n - 1
$$
```{r, echo = FALSE}
pl_n <- 200
pl_x <- seq(-1, 1, len = pl_n)
pl_y <- sqrt(1 - pl_x**2)
node_N <- 9
node_i <- seq(0, node_N-1)
node_x <- -cos(pi * node_i / (node_N - 1))
node_y <- sqrt(1 - node_x**2)
nodes <- cbind(node_x, node_y)
plot(pl_x, pl_y, xlab = "", ylab = "", xlim = c(-1, 1), ylim = c(0, 1),
asp = 1, type = "l", lwd = 2, col = "blue")
points(cbind(node_x, 0), pch = 19, col = "blue")
invisible(sapply(0:node_N, function(i) lines(c(0, nodes[i, 1]), c(0, nodes[i, 2]), lty = 2)))
invisible(sapply(1:(node_N-1), function(i) lines(rep(nodes[i, 1], 2), c(0, nodes[i, 2]), lty = 2, col = "blue")))
```
### Exercise
Modify your program in such a way, that interpolation nodes are generated according to the above defined formula for $x_i$.
Check,how interpolation errors change when functions: $|x|$ or $1/(1+10x^2)$ are used to generate the values $y_i$ in interpolation nodes distributed (non-)uniformly.
Can you see a difference ?