Scienti c Computing Tutorial : Interpolation and Curve Fitting

Universiteit van Stellenbosch
Departement Wiskundige Wetenskappe
Afdeling Rekenaarwetenskap
University of Stellenbosch
Department of Mathematical Sciences
Computer Science Division
Scientific Computing 372
Tutorial 6: Interpolation and Curve Fitting
For submission
1. Background: Newton’s method of interpolation involves two steps: Computation of the coefficients,
followed by evaluation of the polynomial. This works wel if the interpolation is carried out repeatedly
at different values of x using the same polynomial. If only point is to be interpolated, a method that
computes the interpolant in a single step, such as Neville’s algorithm, is a better choice.
Let Pk [x i , x i+1 , . . . , x i+k ] denote the polynomial of degree k that passes through the k + 1 data points
(x i , y i ), (x i+1 , y i+1 ), . . . , (x i+k , y i+k ). For a single data point, we have
P0 [x i ] = y i .
(1)
The interpolant based on two data points is
P1 [x i , x i+1 ] =
(x − x i+1 )P0 [x i ] + (x i − x)P0 [x i+1 ]
.
x i − x i+1
It is easily verified that P1 [x i , x i+1 ] passes through the two data points; that is, P1 [x i , x i+1 ] = y i when
x = x i , and P1 [x i , x i+1 ] = y i+1 when x = x i+1 .
The three-point interpolant is
P2 [x i , x i+1 , x i+2 ] =
(x − x i+2 )P1 [x i , x i+1 ] + (x i − x)P1 [x i+1 , x i+2 ]
.
x i − x i+2
To show that this interpolant does intersect the data points, we first substitute x = x i , obtaining
P2 [x i , x i+1 , x i+2 ] = P1 [x i , x i+1 ] = y i .
Similarly, x = x i+2 yields
P2 [x i , x i+1 , x i+2 ] = P1 [x i+1 , x i+2 ] = y i+2 .
Finally, when x = x i+1 , we have
P1 [x i , x i+1 ] = P1 [x i+1 , x i+2 ] = y i+1 ,
so that
P2 [x i , x i+1 , x i+2 ] =
(x i+1 − x i+2 )y i+1 − (x i − x i+1 )y i+1
= y i+1 .
x i − x i+2
Having established the pattern, we can now deduce the general recursive formula:
Pk [x i , x i+1 , . . . , x i+k ] =
(x − x i+k )Pk−1 [x i , x i+1 , . . . , x i+k−1 ] + (x i − x)Pk−1 [x i+1 , x i+2 , . . . , x i+k ]
.
x i − x i+k
(2)
Assignment: Implement Neville’s method in a Python starts with the signature given below, and save
your work to a file called neville.py.
Scientific Computing 372
Tutorial 6
Stellenbosch University, 2014
def neville(x_data, y_data, x):
"""p = neville(x_data, y_data, x)
Evaluate the polynomial interpolant p(x) that passes through the
specified data points by Neville's method."""
Questions for test preparation
2. (a) Use the natural cubic spline to determine y at x = 1.5. The data points are as follows:
x
y
1 2 3 4 5
0 1 0 1 0
(b) Find the zero of y(x) from the following data:
y
x
0
0.5
1
1.5
2
2.5
3
1.8421 2.4694 2.4921 1.9047 0.8509 −0.4112 −1.5727
(c) Use Neville’s method to compute y at x = π/4 from the following data points:
x
y
0
0.5
1
1.5
2
−1.00 1.75 4.00 5.75 7.00
(d) Given the the following data, find y at x = π/4 and x = π/2. Use the method that you consider to
be most convenient.
x
y
0
0.5
1
1.5
2
−0.7854 0.6529 1.7390 2.2071 1.9425
(e) The table shows the drag coefficient c D of a sphere as a function of the Reynolds number Re. Use
the natural cubic spline to find c D at Re = 5, 50, 500, and 5000. Hint: Use the log-log scale.
Re
cD
0.2
2
20
200 2000 20 000
103 13.9 2.72 0.800 0.401 0.433
3. (a) Fit a straight line to the following data, and compute the standard deviation.
x
y
0.0 1.0 2.0 2.5 3.0
2.9 3.7 4.1 4.4 5.0
(b) The relative density ρ of air was measured at various altitudes h (in km). The results are:
h
ρ
0 1.525 3.050 4.575
6.10
7.625 9.150
1 0.8617 0.7385 0.6292 0.5328 0.4481 0.3741
Use a quadratic least-squares fit to determine the relative air density at h = 10.5 km.
Page 2
Scientific Computing 372
Tutorial 6
Stellenbosch University, 2014
(c) Fit a straight line and quadratic to the follow data.
x
y
1.0
2.5
3.5
4.0
1.1
1.8
2.2
3.7
6.008 15.722 27.130 33.772 5.257 9.549 11.098 28.828
(d) The intensity of radiation of a radioactive substance was measured at half-year intervals. The
results are
t
γ
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
5.5
1.000 0.994 0.990 0.985 0.979 0.977 0.972 0.969 0.967 0.960 0.956 0.952
where t is the time in years and γ is the relative intensity of radiation. Knowing that radioactivity
decays exponentially with time, γ(t) = ae −bt , estimate the radioactive half-life of the substance.
Page 3

Download Report