Download Linear Interpolation: Determining Function Values between Tabulated Data and more Schemes and Mind Maps Physics in PDF only on Docsity!
4 Interpolation
4.1 Linear Interpolation
A common computational problem in physics involves determining the value of a particular function at one or more points of interest from a tabulation of that function. For instance, we may wish to calculate the index of refraction of a type of glass at a particular wavelength, but be faced with the problem that that particular wavelength is not explicitly in the tabulation. In such cases, we need to be able to interpolate in the table to find the value of the function at the point of interest. Let us take a particular example. BK-7 is a type of common optical crown glass. Its index of refraction n varies as a function of wavelength; for shorter wavelengths n is larger than for longer wavelengths, and thus violet light is refracted more strongly than red light, leading to the phenomenon of dispersion. The index of refraction is tabulated in the following table:
Refractive Index for BK7 Glass λ (˚A) n λ (˚A) n λ (˚A) n 3511 1.53894 4965 1.52165 8210 1. 3638 1.53648 5017 1.52130 8300 1. 4047 1.53024 5145 1.52049 8521 1. 4358 1.52669 5320 1.51947 9040 1. 4416 1.52611 5461 1.51872 10140 1. 4579 1.52462 5876 1.51680 10600 1. 4658 1.52395 5893 1.51673 13000 1. 4727 1.52339 6328 1.51509 15000 1. 4765 1.52310 6438 1.51472 15500 1. 4800 1.52283 6563 1.51432 19701 1. 4861 1.52238 6943 1.51322 23254 1. 4880 1.52224 7860 1.
Let us suppose that we wish to find the index of refraction at a wavelength of 5000˚A. Unfortunately, that wavelength is not found in the table, and so we must estimate it from the values in the table. We must make some assumption about how n varies between the tabular values. Presumably it
varies in a smooth sort of way and does not take wild excursions between the tabulated values. The simplest and quite often an entirely adequate assumption to make is that the actual function varies linearly between the tabulated values. This is the basis of linear interpolation.
Exercise 4.1: Determine, by hand, the value of the index of refraction of BK7 at 5000˚A using linear interpolation.
How do we carry out linear interpolation on the computer? Let us sup- pose that the function is tabulated at N points and takes on the values y 1 , y 2 , y 3 ...yN at the points x 1 , x 2 , x 3 ...xN , and that we want to find the value of the function y at a point x that lies someplace in the interval between x 1 and xN. The first thing that we must do is to bracket x, that is we must find a j such that xj < x ≤ xj+1. This can be accomplished by the following code fragment:
for(i=1;i<N;i++) { if(xn[i] < x && xn[i+1] >= x) { j = i; break; } }
where the xn’s are the tabulated points. When the if statement is satisfied, j is assigned the value of i and the procedure drops out of the loop. Please note that this is not the most efficient way to accomplish this task, especially if N is very large. We will look at a more efficient way later on. Once we have bracketed x, we can find the equation of the line between the points (xj , yj ) and (xj+1, yj+1). This equation will be of the form y = mx + b where m is the slope and b is the y-intercept. As we all know, the slope is given by
m =
yj+1 − yj xj+1 − xj
4.2 Polynomial Interpolation
Linear interpolation is good enough for government work, and it is even better than that. Because it is simple and makes the simplest possible assumption about the data, it should be employed in all cases except where it is manifestly inadequate. There are such cases. Sometimes the function being interpolated is very non-linear or has been tabulated at such wide intervals that linear interpolation would lead to large errors. Some applications demand more than simply the functional values at the interpolated points; sometimes the derivative of the function is required as well. With linear interpolation, the derivative is a constant between the tabulated points, and may actually be undefined at the tabulated points! For such applications it may be best to interpolate using a polynomial interpolating function or functions. It can be shown that the following poly- nomial P (x) of degree N −1 will exactly pass through the N tabulated points of the function y = f (x):
P (x) =
(x − x 2 )(x − x 3 )... (x − xN ) (x 1 − x 2 )(x 1 − x 3 )... (x 1 − xN )
y 1
(x − x 1 )(x − x 3 )... (x − xN ) (x 2 − x 1 )(x 2 − x 3 )... (x 2 − xN )
y 2
(x − x 1 )(x − x 2 )...(x − xN − 1 ) (xN − x 1 )(xN − x 2 )...(xN − xN − 1 )
yN
The problem with the direct application of the polynomial P (x) is that for tabulations with many points, it can lead to very high degree polynomi- als. For instance, if a function is tabulated at 100 points, the above equation would yield a polynomial of degree 99! Such a polynomial could potentially fluctuate wildly between the tabulated points and thus not be a good rep- resentation of the actual function. P (x) is more usually applied to subsets of the tabulated points. For instance, if the polynomial is applied to subsets of 3 points, it yields parabolic interpolation, which can be much superior to linear interpolation if the function has a number of minima and maxima. One problem with parabolic 3-point interpolation is that when it comes to bracketing x, there is an ambiguity – does one bracket between x 1 and x 2 or between x 2 and x 3? This is one reason why cubic 4-point interpolation is more commonly practiced – the bracketing is then between x 2 and x 3 with no
ambiguity. Such interpolation is also called Lagrangian 4-point interpolation. The Lagrangian 4-point interpolation equation can be written:
L(x) =
(x − x 2 )(x − x 3 )(x − x 4 ) (x 1 − x 2 )(x 1 − x 3 )(x 1 − x 4 )
y 1 +
(x − x 1 )(x − x 3 )(x − x 4 ) (x 2 − x 1 )(x 2 − x 3 )(x 2 − x 4 )
y 2 +
(x − x 1 )(x − x 2 )(x − x 4 ) (x 3 − x 1 )(x 3 − x 2 )(x 3 − x 4 )
y 3 +
(x − x 1 )(x − x 2 )(x − x 3 ) (x 4 − x 1 )(x 4 − x 2 )(x 4 − x 3 )
y 4
which the student can easily verify by reference to the equation for P (x). To apply this interpolation equation, the user should first bracket x be- tween xj and xj+1 as before, but now identify xj with x 2 in the above equation and xj+1 with x 3. It then follows that x 1 will be xj− 1 and x 4 will be xj+2. The perceptive student will see that this will lead to a problem at the end- points. For instance, if x is situated between the first two tabulated points, xj− 1 will be undefined. Likewise, if x is situated between the last two tabu- lated points, xj+2 will be undefined. Thus in those intervals, the user must either interpolate linearly, or, in the first case, use the polynomial that would be used for an x bracketed between the 2nd and 3rd tabulated points, and similarly for the last case.
Exercise 4.4: Write a C-function that will implement the 4-point La- grangian interpolation formula above. For the endpoints, use the polynomial that would have been defined for the adjacent interval as described above. Modify the driver program in Exercise 4.2 (interpolation in a table of the wavelength and the index of refraction for BK-7 glass) to use this new C- function. Compare the results between the two programs.
Exercise 4.5: Write a C-function that will implement the 4-point La- grangian interpolation formula above. For the endpoints, use the polynomial that would have been defined for the adjacent interval as described above. Modify the driver program in Exercise 4.3 (interpolation in a table of atmo- spheric pressure and the boiling point of water) to use this new C-function. Compare the results between the two programs.
independent variables. The functional value y is tabulated on a Cartesian grid, and so the first thing the programmer must do is to bracket the desired point (x 1 , x 2 ) in this grid (see figure below):
x
x
x
x
(x 1 ,i, x 2 ,j ) (x 1 ,i+1, x 2 ,j )
(x 1 ,i, x 2 ,j+1) (x 1 ,i+1, x 2 ,j+1)
h(x 1 , x 2 )
This bracketing can be done by bracketing in the two dimensions one at at time, using the technique discussed earlier. The simplest interpolation technique in two dimensions is bilinear inter- polation. If we define
y 1 = y(x 1 ,i, x 2 ,j ) y 2 = y(x 1 ,i+1, x 2 ,j ) y 3 = y(x 1 ,i+1, x 2 ,j+1) y 4 = y(x 1 ,i, x 2 ,j+1)
i.e., working our way counterclockwise around the above figure, then the interpolation formulae are:
t = (x 1 − x 1 ,i)/(x 1 ,i+1 − x 1 ,i) u = (x 2 − x 2 ,j )/(x 2 ,j+1 − x 2 ,j )
which make both t and u lie between 0 and 1. Then,
y(x 1 , x 2 ) = (1 − t)(1 − u)y 1 + t(1 − u)y 2 + tuy 3 + (1 − t)uy 4
where (x 1 , x 2 ) are the coordinates of the desired point.
Exercise 4.6: A good example of the need to carry out interpolation in 2 dimensions is found in the calculation of partition functions. In certain plasma-physics contexts, it is necessary to calculate the partition functions of atoms and ions. A partition function, U, is essentially an overall “statistical weight” for the atom, calculated by carrying out a weighted sum – weighted according to the populations of the levels – of the statistical weights for all of the energy levels in the atom. At low temperatures, the partition function is simply the statistical weight of the ground level of the atom, but at higher temperatures, the partition function becomes larger, as the populations in the excited levels become significant. The partition function is also a function of density, as at high densities in the plasma, the outermost energy levels are effectively “stripped off”, resulting in a lowering of the ionization energy, usually denoted as ∆E. Thus, U = U(T, ∆E). Because the calculation for a partition function can be very complex, they are usually calculated and tabulated so that users need not carry out the full calculation. The following table is a tabulation for the partition function for the hydrogen atom:
The Hydrogen Partition Function T(K) ∆E = 0. 10 ∆E = 0. 50 ∆E = 1. 00 ∆E = 2. 00 3250 2.000 2.000 2.000 2. 10083 2.000 2.000 2.000 2. 14188 2.025 2.006 2.005 2. 15643 2.068 2.016 2.012 2. 17246 2.168 2.037 2.027 2. 19014 2.384 2.080 2.058 2. 20963 2.814 2.162 2.114 2. 23111 3.610 2.308 2.213 2. 25480 4.991 2.551 2.377 2.
Write a function that will perform a 2-D interpolation in this table, and use it to determine the partition function for Hydrogen at the following points (T, ∆E): (16000, 0.25), (18500, 1.50), (19000, 0.15), (25023, 1.99).
Exercise 4.7: The study of plasmas is an important field of physics, and is essential in the understanding of the interiors of stars and the functioning of fusion reactors. Hydrogen becomes increasingly ionized (hydrogen loses
