Cubic interpolation online. Spline interpolation. Spline interpolation: an example of building a spline in the STATISTICA program

20.06.2022 diets

MINISTRY OF EDUCATION AND SCIENCE OF THE RUSSIAN FEDERATION

Federal State Autonomous Educational Institution

higher professional education

"Ural Federal University named after the first President of Russia B.N. Yeltsin"

Institute of Radio Electronics and Information Technologies - RTF

department Automation and information technology

Spline interpolation

METHODOLOGICAL INSTRUCTIONS For laboratory work ON THE DISCIPLINE "Numerical methods"

Compiled by I.A. Selivanova, senior lecturer.

SPLINE INTERPOLATION: Guidelines for practical exercises in the discipline "Numerical Methods"

The instructions are intended for students of all forms of education in the direction 230100 - "Informatics and Computer Engineering".

Ó FSAEI HPE "UrFU named after the first President of Russia B.N. Yeltsin", 2011

1. INTERPOLATION BY SPLINES. four

1.1. Cubic splines. four

1.2. A special form of writing a spline. 5

1.3. Quadratic splines. 13

1.4. Practice assignment. eighteen

1.5. Task options. 19

References 21

1. Interpolation by splines.

In cases where the interval [ a,b], on which it is required to replace the function f(x) large, you can apply spline interpolation.

1.1. Cubic splines.

Interpolation splines 3rd order are functions consisting of pieces of polynomials 3 th order. The conjugation nodes ensure the continuity of the function, its first and second derivatives. The approximating function is composed of separate polynomials, as a rule, of equally small degree, each defined on its own part of the segment.

Let on the segment [ a, b] real axis x the grid is given, in the nodes of which the values are defined
functions f(x). It is required to build on the segment [ a, b] continuous spline function S(x), which satisfies the following conditions:

To construct the desired spline, you need to find the coefficients
polynomials
,i=1,… n, i.e. 4 n unknown coefficients that satisfy 4 n-2 equations (1), (2), (3). In order for the system of equations to have a solution, two more additional (boundary) conditions are added. Three types of boundary conditions are used:

Conditions (1), (2), (3) and one of the conditions (4), (5), (6) form a SLAE of order 4 n. The system can be solved using the Gauss method. However, by choosing a special form of writing a cubic polynomial, one can significantly reduce the order of the system of equations being solved.

1.2. A special form of writing a spline.

Consider the segment
. Let us introduce the following notation for variables:

Here
- segment length
,

,
- auxiliary variables,

x- intermediate point on the segment
.

When x runs through all values in the interval
, variable changes from 0 to 1, and
changes from 1 to 0.

Let the cubic polynomial
on the segment
looks like:

Variables and
are determined in relation to a specific segment of interpolation.

Find the value of the spline
at the ends of the segment
. Dot
is the initial for the segment
, that's why =0,
=1 and according to (3.8):
.

At the end of the segment
=1,
=0 and
.

For interval
dot
is final, so =1,
=0 and from formula (9) we get:
. Thus, the continuity condition for the function is satisfied S(x) at the junction points of cubic polynomials, regardless of the choice of numbers  i .

To determine the coefficients  i , i=0,… n we differentiate (8) twice as a complex function of x. Then

Define the second derivatives of the spline
and
:

For polynomial
dot is the beginning of the interpolation segment and =0,
=1, so

From (15) and (16) it follows that on the segment [ a,b]spline function, "glued" from pieces of polynomials of the 3rd order, has a continuous derivative of the 2nd order.

To obtain the continuity of the first derivative of a function S(x), we require the following conditions to be fulfilled in the internal nodes of interpolation:

For natural cubic spline
, therefore, the system of equations will look like:

and the system of equations (17) will look like:

Example.

Initial data:

Replace function
interpolating cubic spline, the values of which at the given nodal points (see table) coincide with the values of the function at the same points. Consider different boundary conditions.

Let's calculate the value of the function at the nodal points. To do this, we substitute the values from the table into the given function.

For different boundary conditions (4), (5), (6) we find the coefficients of cubic splines.

Consider the first boundary conditions.

In our case n=3,
,
,
. To find
we use the system of equations (3.18):

Compute and , using formulas (7) and (11):

We substitute the obtained values into the system of equations:

System solution:

Taking into account the first boundary conditions, the spline coefficients:

Consider the definition of the spline coefficients taking into account the boundary conditions (3.5):

Let's find the derivative of the function
:

Compute
and
:

Let us substitute into the system of equations (21) the values and :

Using formula (20), we determine  0 and  3:

Given specific values:

and coefficient vector:

Let's calculate the values of the cubic spline S(x) at the midpoints of the interpolation segments.

Middle sections:

To calculate the value of the cubic spline at the midpoints of the interpolation segments, we use formulas (7) and (9).

3.1.

Let's find and
:

In formula (3.9) we substitute the coefficients

3.2.

Let's find and
:

, for boundary conditions (4), (5), (6):

3.3.

Let's find and
:

In formula (9) we substitute the coefficients
, for boundary conditions (4), (5), (6):

Let's make a table:

(1 cr. condition.)	(2 cr. conditions)	(3 cr. conditions)

Curves and surfaces encountered in practical problems often have a rather complex shape, which does not allow a universal analytical specification as a whole with the help of elementary functions. Therefore, they are assembled from relatively simple smooth fragments - segments (curves) or cuts (surfaces), each of which can be quite satisfactorily described using elementary functions of one or two variables. In this case, it is quite natural to require that smooth functions that are used to construct partial curves or surfaces have a similar nature, for example, be polynomials of the same degree. And in order for the resulting curve or surface to be sufficiently smooth, it is necessary to be especially careful at the junctions of the corresponding fragments. The degree of polynomials is chosen from simple geometric considerations and, as a rule, is small. For a smooth change of the tangent along the entire compound curve, it is sufficient to describe the joining curves using polynomials of the third degree, cubic polynomials. The coefficients of such polynomials can always be chosen so that the curvature of the corresponding composite curve is continuous. Cubic splines that arise in solving one-dimensional problems can be adapted to the shaping of fragments of composite surfaces. And here, quite naturally, there appear bicubic splines, described by polynomials of the third degree in each of the two variables. Working with such splines requires much more calculations. But a properly organized process will allow taking into account the continuously growing capabilities of computer technology to the maximum extent. Spline functions Let on the segment , that is, Remark. The index (t) of the numbers a^ indicates that. that the set of coefficients by which the function S(x) is determined on each partial segment D is its own. On each of the segments D1, the spline 5(x) is a polynomial of degree p and is determined on this segment by p + 1 coefficient. Total partial segments - then. Hence, in order to completely determine the spline, it is necessary to find (p + 1) then numbers. Condition) means the continuity of the function S(x) and its derivatives at all internal grid nodes w. The number of such nodes is m - 1. Thus, to find the coefficients of all polynomials, p(m - 1) conditions (equations) are obtained. For a complete definition of the spline, there are not enough (conditions (equations). The choice of additional conditions is determined by the nature of the problem under consideration, and sometimes simply by the desire of the user. SPLINE THEORY Examples of solutions Most often, interpolation and smoothing problems are considered, when it is required to build one or another spline from a given array of points on a plane In interpolation problems, it is required that the spline graph pass through points, which imposes m + 1 additional conditions (equations) on its coefficients. The remaining p - 1 conditions (equations) for the unique construction of a spline are most often set in the form of values of the lower derivatives of the spline at the ends of the segment under consideration [a, 6] - boundary (boundary) conditions. The ability to select different boundary conditions allows you to build splines with a variety of properties. In smoothing problems, a spline is built so that its graph passes near the points (i "" Y "), * = 0, 1, ..., m, and not through them. The measure of this closeness can be defined in different ways, which leads to a significant variety of smoothing splines. The described options for choosing when constructing spline functions are far from exhausting their diversity. And if initially only piecewise polynomial spline functions were considered, then as the scope of their applications expanded, splines began to appear, "glued" from other elementary functions as well. Interpolation cubic splines Statement of the interpolation problem Let the grid w be given on the interval [a, 6) Consider a set of numbers Problem. Construct a function that is smooth on the segment (a, 6] and takes on the given values at the nodes of the grid o, i.e. "By imposing additional conditions on the function being constructed, one can achieve the necessary uniqueness. In applications, it often becomes necessary to approximate a function given analytically by means of a function with prescribed sufficiently good properties. For example, in cases where the calculation of the values of a given function f(x) at points segment [a, 6] is associated with significant difficulties and/or the given function f(x) does not have the required smoothness, it is convenient to use another function that would approximate the given function well enough and would be devoid of its shortcomings. [a, 6] a smooth function a(x) coinciding at the grid nodes w with the given function /(X). Definition of an interpolating cubic spline An interpolating cubic spline S(x) on a mesh w is a function that 1) on each of the segments is a polynomial of the third degree, 2) is twice continuously differentiable on the segment [a, b], that is, belongs to the class C2[ a, 6], and 3) satisfies the conditions On each of the segments, the spline S(x) is a polynomial of degree three and is determined on this segment by four coefficients. The total number of segments is m. This means that in order to completely determine the spline, it is necessary to find 4m numbers. The condition means the continuity of the function S (x) and its derivatives S "(x) and 5" (x) at all internal grid nodes w. The number of such nodes is m - 1. Thus, to find the coefficients of all polynomials, 3 (m - 1) more conditions (equations) are obtained. Together with conditions (2), conditions (equations) are obtained. Boundary (boundary) conditions Two missing conditions are specified as restrictions on the values of the spline and/or its derivatives at the ends of the interval [a, 6]. When constructing an interpolating cubic spline, the boundary conditions of the following four types are most often used. A. Boundary conditions of the 1st type. - at the end of the interval [a, b], the values of the first derivative of the desired function are given. B. Boundary conditions of the 2nd type. - at the end of the interval (a, 6) the values of the second derivative of the desired function are set. B. Boundary conditions of the 3rd type. are called periodic. It is natural to require the fulfillment of these conditions in cases where the interpolated function is periodic with period T = b-a. D. Boundary conditions of the 4th type. require special comment. Comment. At internal sepsi nodes, the third derivative of the function S(x) is, generally speaking, discontinuous. However, the number of discontinuities of the third derivative can be reduced by using conditions of the 4th type. In this case, the constructed spline will be continuously differentiable three times on intervals. Construction of an interpolating cubic spline Let us describe a method for calculating the coefficients of a cubic spline, in which the number of quantities to be determined is equal. On each of the intervals, the interpolation spline function is sought in the following form For boundary conditions of the 1st and 2nd types, this system has the following form where Coefficients depend on the choice of boundary conditions. Boundary conditions of the 1st type: Boundary conditions of the 2nd type: In the case of boundary conditions of the 3rd type, the system for determining numbers is written as follows. For boundary conditions of the 4th type, the system for determining numbers has the form The matrices of all three linear algebraic systems are matrices with diagonal dominance. These matrices are not degenerate, and therefore each of these systems has a unique solution. Theorem. An interpolation cubic spline that satisfies conditions (2) and a boundary condition of one of the listed four types exists and is unique. Thus, to construct an interpolating cubic spline means to find its coefficients. When the coefficients of the spline are found, the value of the spline S(x) at an arbitrary point of the segment [a, b] can be found using formula (3). However, for practical calculations, the following algorithm for finding the quantity S(x) is more suitable. Let x 6 [x", First, the values A and B are calculated according to the formulas and then the value 5(x) is found: The use of this algorithm significantly reduces the computational costs for determining the value Advice to the user The choice of boundary (boundary) conditions and interpolation nodes allows to a certain extent to control properties of interpolation splines. A. Choice of boundary (boundary) conditions. The choice of boundary conditions is one of the central problems in the interpolation of functions. It acquires special importance in the case when it is necessary to ensure high accuracy of the approximation of the function f(x) by the spline 5(g) near the ends of the segment [a, 6]. The boundary values have a noticeable effect on the behavior of the spline 5(g) near the points a and b, and this effect rapidly weakens as we move away from them. The choice of boundary conditions is often determined by the availability of additional information about the behavior of the function f(x) being approximated. If the values of the first derivative f "(x) are known at the ends of the segment (a, 6), then it is natural to use the boundary conditions of the 1st type. If the values of the second derivative f "(x) are known at the ends of the segment [a, 6], then it is natural use boundary conditions of the 2nd type. If it is possible to choose between the boundary conditions of the 1st and 2nd types, then preference should be given to the conditions of the 1st type. If f(x) is a periodic function, then we should stop at the boundary conditions of the 3rd type. If there is no additional information about the behavior of the function being approximated, the so-called natural boundary conditions are often used. However, it should be borne in mind that with such a choice of boundary conditions, the accuracy of approximating the function f (x) by the spline S (x) near the ends of the segment (a, ft] decreases sharply. Sometimes, boundary conditions of the 1st or 2nd type are used, but not with the exact values of the corresponding derivatives, but with their difference approximations. The accuracy of this approach is low. Practical experience of calculations shows that in the situation under consideration, the most appropriate choice is boundary conditions of the 4th type. B. Choice of interpolation nodes. If the third derivative f""(x) of the function suffers a discontinuity at some points of the segment [a, b], then to improve the quality of the approximation, these points should be included in the number of interpolation nodes. If the second derivative /"(x) is discontinuous, then in order to avoid oscillation of the spline near the discontinuity points, special measures must be taken. Usually, the interpolation nodes are chosen so that the discontinuity points of the second derivative fall inside the interval \xif), such that. The value and can be chosen by numerical experiment (it is often sufficient to set a = 0.01). There is a set of recipes for overcoming the difficulties that arise when the first derivative f "(x) is discontinuous. As one of the simplest, we can propose this: divide the approximation segment into intervals where the derivative is continuous, and build a spline on each of these intervals. Choice of interpolation function (pluses and minuses) Approach 1st. Lagrange interpolation polynomial According to a given array SPLINE THEORY solution examples (Fig. 3) the Lagrange interpolation polynomial is determined by the formula It is advisable to consider the properties of the Lagrange interpolation polynomial from two opposite positions, discussing the main advantages separately from the disadvantages. The main advantages of the 1st approach: 1) the graph of the Lagrange interpolation polynomial passes through each point of the array, 2) the constructed function is easily described (the number of coefficients of the Lagrange interpolation polynomial on the grid u to be determined is equal to m + 1), 3) the constructed function has continuous derivatives of any order, 4) given an array, the interpolation polynomial is uniquely determined. The main disadvantages of the 1st approach: 1) the degree of the Lagrange interpolation polynomial depends on the number of grid nodes, and the larger this number, the higher the degree of the interpolation polynomial and, therefore, the more calculations are required, 2) changing at least one point in the array requires complete recalculation of the coefficients of the Lagrange interpolation polynomial, 3) adding a new point to the array increases the degree of the Lagrange interpolation polynomial by one and even leads to a complete recalculation of its coefficients, 4) with unlimited mesh refinement, the degree of the Lagrange interpolation polynomial increases indefinitely. The behavior of the Lagrange interpolation polynomial under unlimited mesh refinement generally requires special attention. Comments A. Approximation of a continuous function by a polynomial. It is known (Weierstrass, 1885) that any continuous (and even more so smooth) function on an interval can be approximated as well as desired on this interval by a polynomial. Let us describe this fact in the language of formulas. Let f(x) be a function continuous on the segment [a, 6]. Then for any e > 0 there is a polynomial Рn(x) such that for any x from the interval [a, 6] the inequality will be satisfied (Fig. 4) , there are infinitely many. On the segment [a, 6] we construct a grid w. It is clear that its nodes, generally speaking, do not coincide with the intersection points of the graphs of the polynomial Pn(x) and the function f(x) (Fig. 5). Therefore, for the grid taken, the polynomial Pn(x) is not an interpolation polynomial. When a continuous function is approximated by an Jla-grajj interpolation polynomial, its graph not only does not have to be close to the graph of the function f(x) at every point of the interval [a, b), but can deviate from this function as much as desired. Let's give two examples. Example 1 (Rung, 1901). With an unlimited increase in the number of nodes for a function on the interval [-1, 1], the limit equality is fulfilled (Fig. 6) Example 2 (Berichtein, 1912). A sequence of Lagrange interpolation polynomials constructed on uniform grids nm for a continuous function /(x) = |x| on the segment with an increase in the number of nodes m does not tend to the function f(x) (Fig. 7). Approach 2nd. Piecewise linear interpolation If the smoothness of the interpolated function is abandoned, the ratio between the number of advantages and the number of disadvantages can be noticeably changed in the direction of the former. Let us construct a piecewise linear function by successively connecting points (xit y,) with straight line segments (Fig. 8). The main advantages of the 2nd approach are: 1) the graph of a piecewise linear function passes through each point of the array, 2) the constructed function is easily described (the number of coefficients of the corresponding linear functions to be determined for the grid (1) is 2m), 3) the constructed function is defined by a given array unambiguously, 4) the degree of polynomials used to describe the interpolation function does not depend on the number of grid nodes (equal to 1), 5) changing one point in the array requires the calculation of four numbers (the coefficients of two rectilinear links emanating from the new point), 6) adding additional point to the array requires the calculation of four coefficients. The piecewise linear function behaves quite well when refining the grid. i The main drawback of the 2nd approach is that the approximating piecewise linear function is not smooth: the first derivatives suffer a discontinuity at the grid nodes (interpolation ears). Approach 3rd. Spline interpolation The proposed approaches can be combined so that the number of listed advantages of both approaches is preserved while reducing the number of disadvantages. This can be done by constructing a smooth interpolating spline function of degree p. The main advantages of the 3rd approach: 1) the graph of the constructed function passes through each point of the array, 2) the constructed function is relatively easy to describe (the number of coefficients of the corresponding polynomials to be determined for the grid (1) is 3) the constructed function is uniquely determined by a given array, 4) the degree polynomials does not depend on the number of grid nodes and, therefore, does not change with its increase, 5) the constructed function has continuous derivatives up to the order p - 1 inclusive, 6) the constructed function has good approximation properties. Brief reference. The proposed name - spline - is not accidental - the smooth piecewise polynomial functions introduced by us and drawing splines are closely related. Consider a flexible, ideally thin ruler passing through the reference points of the array located on the (x, y) plane. According to the Bernoulli-Euler law, the linearized equation of a curved ruler has the form The function S(x), which describes the rulers, is a polynomial of the third degree between each and two neighboring points of the array (supports) and is twice continuously differentiable on the entire interval (a, 6). Comment. 06 Interpolation of a continuous function In contrast to the Lagrange interpolation polynomials, a sequence of interpolation cubic splines on a uniform grid always converges to an interpolated continuous function, and with the improvement of the differential properties of this function, the rate of convergence increases. Example. For a function, a cubic spline on a grid with the number of nodes m = 6 gives an approximation error of the same order as the interpolation polynomial Ls(z), and on a grid with the number of nodes m = 21 this error is so small that on the scale of an ordinary book drawing it simply cannot be can be shown (Fig. 10) (the interpolation polynomial 1>2o(r) gives in this case an error of about 10,000 W). Properties of an interpolated cubic spline A. Approximation properties of a cubic spline. The approximation properties of the interpolating spline depend on the smoothness of the function f(x) - the higher the smoothness of the interpolated function, the higher the order of approximation, and when the mesh is refined, the higher the convergence rate. If the interpolated function f(x) is continuous on the interval If the interpolated function f(x) has a continuous first derivative on the interval [a, 6], that is, an interpolation spline that satisfies the boundary conditions of the 1st or 3rd type, then for h we have In this case, not only the spline converges to the interpolated function, but also the derivative of the spline converges to the derivative of this function. If the spline S(x) approximates the function f(x) on the segment [a, b], and its first and second derivatives approximate the function B, respectively. Extremal property of a cubic spline. The interpolating cubic spline has another useful property. Consider the following example. example. Construct a function /(x) minimizing the functional on the class of functions from the space C2 whose graphs pass through the points of the array x) that satisfies the boundary conditions delivers an extremum (minimum) to the functional. Remark 2. It is interesting to note that the interpolating cubic spline has the extremal property described above on a very wide class of functions, namely, on the class |0, 5]. 1.2. Smoothing cubic splines On the formulation of the smoothing problem Let a grid and a set of numbers be given. In fact, this means that an interval is specified for each, and any number from this interval can be taken as the value of y, . It is convenient to interpret the values of y, for example, as the results of measurements of some function y(x) for given values of the variable x, containing a random error. When solving the problem of restoring a function from such "experimental" values, it is hardly advisable to use interpolation, since the interpolation function will obediently reproduce bizarre oscillations caused by a random component in the array (y,). A more natural approach is based on a smoothing procedure designed to somehow reduce the element of randomness as a result of measurements. Usually in such problems it is required to find a function whose values for x = x, * = 0, 1, .... m, would fall into the corresponding intervals and which, in addition, would have sufficiently good properties. For example, it would have continuous first and second derivatives, or its graph would not be too strongly curved, that is, it would not have strong oscillations. A problem of this kind also arises when, according to a given (exactly) array, it is required to construct a function that would pass through non-given points, but near them and, moreover, change quite smoothly. In other words, the desired function smoothed the given array, as it were, and did not interpolate it. Let a grid w and two sets of numbers be given. SPLINE THEORY examples of solutions Problem. Construct a smooth function on the segment [a, A] whose values at the nodes of the grid and differ from the numbers y by the given values. The formulated smoothing problem is recovery smooth function given in a table. It is clear that such a problem has many different solutions. By imposing additional conditions on the constructed function, we can achieve the necessary uniqueness. Definition of a smoothing cubic spline A smoothing cubic spline S(x) on a mesh w is a function that 1) on each of the segments is a polynomial of the third degree, 2) is twice continuously differentiable on the segment [a, 6], that is, belongs to the class C2 [a , b], 3) delivers a minimum to the functional where are given numbers, 4) satisfies the boundary conditions of one of the three types indicated below. Boundary (boundary) conditions Boundary conditions are specified as restrictions on the values of the spline and its derivatives at the boundary nodes of the mesh w. A. Boundary conditions of the 1st type. - at the end of the interval [a, b) the values of the first derivative of the desired function are given. Boundary conditions of the 2nd type. - the second derivatives of the desired function at the ends of the interval (a, b] are equal to zero. B. Boundary conditions of the 3rd type are called periodic. Theorem. Cubic spline S (x), minimizing the functional (4) and satisfying the boundary conditions of one of the three indicated types is uniquely defined. Definition. A cubic spline that minimizes the functional J(f) and satisfies the boundary conditions of the i-type is called a smoothing spline of the i-type. this segment by four coefficients.Total segments - m.So, in order to completely define the spline, you need to find 4m numbers. The condition means the continuity of the function 5(ar) and all derivatives at all internal nodes of the grid o. "The number of such nodes is m - 1 Thus, to find the coefficients of all polynomials, 3(m - 1) conditions (equations) are obtained. for which the number of quantities to be determined is 2m + 2. On each of the intervals, the smoothing spline function is sought in the following form Let us first describe how the quantities n* are found. For boundary conditions of the 1st and 2nd types, the system of linear equations for determining the values of Hi is written in the following form where are the known numbers). The coefficients depend on the choice of boundary conditions. Boundary conditions of the 1st type: Boundary conditions of the 2nd type: In the case of boundary conditions of the 3rd type, the system for determining numbers is written as follows: moreover, all coefficients are calculated by formulas (5) (the quantities with indices k and m + k are considered equal to : Important* note. The matrices of systems are not degenerate, and therefore each of these systems has a unique solution. If the numbers n, - are found, then the quantities are easily determined by the formulas If everything and the smoothing spline turns out to be an interpolation one. This means, in particular, that the more precisely the values are given, the smaller the prescale value of the corresponding weighting coefficients. If, on the other hand, it is necessary that the spline pass through the point (x^, yk), then the weight factor p\ corresponding to it must be set equal to zero. In practical calculations, the most important is the choice of values pi-Let D, - the measurement error of the value y,. Then it is natural to require that the smoothing spline satisfies the condition or, which is the same. In the simplest case, the weight coefficients pi can be given, for example, in the form - where c is some sufficiently small constant. However, such a choice of weights p, does not allow using the "corridor" due to errors in the values of y, -. A more rational, but also more time-consuming algorithm for determining the values of p, - may look as follows. If the values are found at the fc-th iteration, then it is assumed where e is a small number, which is chosen experimentally taking into account the bit grid of the computer, the values of D, and the accuracy of solving the system of linear algebraic equations. If at the fc-th iteration at the point i, condition (6) is violated, then the last formula will ensure a decrease in the corresponding weight coefficient p,. If then, at the next iteration, an increase in p, leads to a more complete use of the "corridor" (6) and, ultimately, a more smoothly changing spline. A bit of theory A. Substantiation of formulas for calculating the coefficients of the interpolation cubic spline. We introduce the notation where m, are unknown quantities. Their number is equal to m + 1. The spline, written in the form where it satisfies the interpolation conditions and is continuous on the entire interval [a, b\: putting in the formula, we obtain, respectively. In addition, it has a continuous first derivative on the interval [a, 6]: differentiating relation (7) and setting, we get the corresponding. actually. Let us show that the numbers m can be chosen so that the spline function (7) has a continuous second derivative on the interval [a, 6]. Calculate the second derivative of the spline on the interval: At the point x, - 0 (at t = 1) we have Calculate the second derivative of the spline on the interval At the point we have From the condition of continuity of the second derivative at the internal grid nodes a; we obtain the m - 1 relation where Adding to these m - 1 equations two more, arising from and from the boundary conditions, we obtain a system of m + 1 linear algebraic equations with m + I unknown miy i = 0, 1. ... , m. The system of equations for calculating the values of gw in the case of boundary conditions of the 1st and 2nd types has the form where (boundary conditions of the 1st type), (boundary conditions of the 2nd type). For periodic boundary conditions (boundary conditions of the 3rd type), the grid o; lengthen by one more node and assume Then the system for determining the values of r* will have the form continuity at the second and (th - !) th grid nodes. We have From the last two relations, we obtain the missing two equations that correspond to the boundary conditions of the 4th type: Excluding the unknown r0 from the equations, and the unknown pc from the equations, as a result we obtain a system of equations Note that the number of unknowns in this system is equal to r - I. 6. Substantiation of the formulas for calculating the efficiency of a smoothing subic spline. We introduce the notation where Zi and nj are still unknown quantities. Their number is equal to 2m + 2. The spline function written in the form is continuous on the entire interval (a, 6]: putting in this formula, we obtain, respectively. Let us show that the numbers z, and n, can be chosen so that the spline written in the form ( 8), had a continuous first derivative on the interval [a, 6] Calculate the first derivative of the spline S(x) on the interval : At a point, we have From the condition of continuity of the first derivative of the spline at the internal nodes of the grid and --> we obtain an m - 1 relation. It is convenient to write this relationship in matrix form. relation (8) and setting, we obtain, respectively, Yeshe olyu matrix relation is obtained from the condition of the minimum of the functional (4). We have The last two matrix equalities can be considered as a linear system of 2m + 2 linear algebraic equations in 2m + 2 unknowns. Replacing the column r in the first equality with its expression obtained from relation (9), we arrive at the matrix equation SPLINE THEORY examples of solutions for determining the column M. This equation has a unique solution due to the fact that the matrix A + 6HRH7 is always nondegenerate. Finding him, we easily identify Mr. Eamshine. The elements of the trianglemagolal matrices A and H determine n only by the grid parameters u (with steps hi) and do not depend on the values yj. Linear Space of Cubic Spline Functions The set of cubic splines constructed on the segment [a, 6) by the node wcra + l is a linear space of dimension m + 3: 1) the sum of two cubic splines constructed by the grid u> and the product of a cubic spline , built on the grid u>, an arbitrary number more secret are cubic splines built on this grid, 2) any cubic spline built on the grid and from the node is completely determined by m + 1 by the value of the values \u200b\u200bof y "at these nodes and two boundary conditions - just + 3 parameters. Choosing in this space a basis consisting of m + 3 linearly independent splines, we can write an arbitrary cubic spline a(x) as a linear combination of them in a unique way. Comment. Such a spline specification is widely used in computational practice. Particularly convenient is the basis, consisting of the so-called cubic B-splines (basic, or fundamental, splines). The use of D-splines can significantly reduce the requirements for computer memory. L-splines. B -spline of zero degree, built on a number line along the grid w, is the function of the fork B -spline of degree k ^ I, built on a number line along the grid u, is determined by the recursive formula second in\7\x) degrees are shown in Fig. 11 and 12 respectively.B-spline of arbitrary degree k can be different from zero only on a certain segment (defined by k + 2 nodes).It is more convenient to number cubic B-splines so that the spline B,-3* (n) was different from zero on the segment ir,-+2]. Let us give a formula for a cubic spline of the third degree for the case of a uniform grid (with a step A). We have in other cases. A typical plot of a cubic B-spline is presented in Fig. 13. The function a) is twice continuously differentiable on a segment, that is, it belongs to the class C2[a, "), c) is nonzero only on four successive segments extended grid w * mo It is necessary to construct a family of m + 3 cubic B-splines: This family forms a basis in the space of cubic splines on the segment (a, b]. Thus, an arbitrary cubic spline S(z) constructed on the segment |s, 6] of the grid o; from +1 nodes, can be represented on this segment as a linear combination. The coefficients ft of this expansion are uniquely determined by the conditions of the problem. ... In the case when the values \u200b\u200bof the function at the nodes of the grid and the values \u200b\u200bof the first derivative of the function at the ends of the grid "(problem of interpolations with boundary conditions of the first kind), these coefficients are calculated from the system of the following form i and &m+i, we obtain a linear system with unknowns 5q, ... , bm and a three-diagonal matrix.The condition provides diagonal dominance and, therefore, the possibility of using the sweep method to resolve it. interpolation problems Zmmchm* 2. In comparison with the algorithms described in section 1.1, the use of the R-spline in interpolation problems * reduces the amount of information stored, that is, significantly reduces the requirements for computer memory, although it leads to an increase in the number of operations. Construction of spline curves using spline functions Above, arrays were considered, the points of which were numbered so that their abscissas formed a strictly increasing sequence. For example, the case depicted in Fig. 14, when different points of the array have the same abscissa, was not allowed. This circumstance determined both the choice of the class of approximating curves (traffic of functions) and the method of constructing them. However, the method proposed above makes it possible to quite successfully construct an interpolation curve in a more general case, when the numbering of array points and their location on the plane, as a rule, are not related (Fig. 15). Moreover, when posing the problem of constructing an interpolation curve, we can consider the given array to be nonplanar, that is, it is clear that in order to solve this general problem, it is necessary to significantly expand the class of admissible curves, including both closed curves, and curves with self-intersection points, and spatial curves. It is convenient to describe such curves using parametric equations Let us require. additionally, so that the functions have sufficient smoothness, for example, they belong to the class C1 [a, /0] or to the class To find the parametric equations of a curve successively passing through all points of the array, proceed as follows. 1st step. On an arbitrarily taken segment, it is drawn along the three nearest points.

Cubic spline interpolation

In recent years, a new branch of modern computational mathematics has been intensively developing - the theory splines. Splines make it possible to effectively solve the problems of processing experimental dependencies between parameters that have a rather complex structure.

The methods of local interpolation discussed above are essentially the simplest spline of the first degree (for linear interpolation) and the second degree (for quadratic interpolation).

Cubic splines have found the widest practical application, due to their simplicity. The main ideas of the theory of cubic splines were formed as a result of attempts to mathematically describe flexible rails made of elastic material (mechanical splines), which have long been used by draftsmen in cases where it became necessary to draw a sufficiently smooth curve through given points. It is known that a rail made of an elastic material, fixed at certain points and being in an equilibrium position, takes a form in which its energy is minimal. This fundamental property makes it possible to effectively use splines in solving practical problems of processing experimental information.

In general, for a function y=f(x) it is required to find an approximation y= j(x) so that f(x i)= j(x i) at points x = x i , a at other points of the segment [ a, b] values

functions f(x) and j(x) were close to each other. With a small number of experimental points (for example, 6-8), one of the methods for constructing interpolation polynomials can be used to solve the interpolation problem. However, with a large number of nodes, interpolation polynomials become practically unusable. This is due to the fact that the degree of the interpolation polynomial is only one less than the number of experimental values of the functions. It is possible, of course, to divide the segment on which the function is defined into segments containing a small number of experimental points, and for each of them construct interpolation polynomials. However, in this case, the approximating function will have points where the derivative is not continuous, i.e., the graph of the function will contain “break” points.

Cubic splines do not have this shortcoming. Studies of beam theory have shown that a flexible thin beam between two nodes is described quite well by a cubic polynomial, and since it does not collapse, the approximating function must be at least continuously differentiable. This means that the functions j(x), j'(x), j"(x) must be continuous on the segment [ a, b].

Cubic interpolation spline , appropriate for this function f(x) and given nodes x i , called a function y(x), satisfying the following conditions:

1. on each segment [ x i — 1 , x i], i = 1, 2, ..., n function y(x) is a third degree polynomial,

Function y(x), and also its first and second derivatives are continuous on the interval [ a,b],

cubic spline is glued together from polynomials of the third degree, which for i th section are written as follows:

For the entire interval, respectively P cubic polynomials differing in coefficients ai, b i, c i, d i. Most often, the nodes during spline interpolation are evenly spaced, i.e. Xi +1 -Xi = const = h (although this is not required).

It is necessary to find four coefficients under the condition that each polynomial passes through two points (x i,y i) and (x i +1 ,y i +1 ) , which results in the following obvious equations:

The first condition corresponds to the passage of the polynomial through the starting point, the second - through the end point. It is impossible to find all the coefficients from these equations, since there are fewer conditions than the required parameters. Therefore, these conditions are supplemented by the conditions of smoothness of the function (ie, continuity of the first derivative) and smoothness of the first derivative (ie, continuity of the second derivative) at the interpolation nodes. Mathematically, these conditions are written as equalities, respectively, of the first and second derivatives at the end i th and at the beginning ( i+1 )-th plots.

Since and , then

(y’ (x i +1 ) in the end i-th section is equal to u'(Xi +1 ) at first ( i+1 )-th),

(at"(Xi +1 ) in the end i-th section is equal to y" (xi +1 ) at first ( i+1)-th).

The result is a system of linear equations (for all sections) containing 4n - 2 equations with 4n unknowns (unknowns a 1 , a 2 ,…, a n , b 1 ,…, d n - spline coefficients). To solve the system, two boundary conditions of one of the following types are added (1 is more often used):

Joint solution of 4n equations allows finding all 4n coefficients.

To restore the derivatives, one can differentiate the corresponding cubic polynomial on each section. If it is necessary to determine the derivatives at the nodes, there are special techniques that reduce the definition of derivatives to solving a simpler system of equations with respect to the desired second or first order derivatives. An important advantage of cubic spline interpolation is obtaining a function that has the minimum possible curvature. The disadvantages of spline interpolation include the need to obtain a relatively large number of parameters.

Let's solve the interpolation problem using the MathCAD program. To do this, we use the built-in function interp(VS,x,y,z) . Variables x and y set the coordinates of the nodal points, z is a function argument, VS defines the type

boundary conditions at the ends of the interval.

We define interpolation functions for three types of cubic spline

Here cspline (VX , VY) returns a vector VS second derivatives when approaching at reference points to a cubic polynomial;

pspline(VX, VY) returns a vector VS second derivatives when approaching the reference points to the parabolic curve;

lspline(VX, VY) returns a vector VS second derivatives when approaching the reference points of the straight line;

interp(VS, VX, VY, x) returns a value y(x) for given vectors VS, VX, VY and set value x.

We calculate the values of interpolation functions at given points and compare the results with the exact values

Note that the results of interpolation by different types of cubic splines are practically the same at the internal points of the interval and coincide with the exact values of the function. Near the edges of the interval, the difference becomes more noticeable, and when extrapolated outside the given interval, different types of splines give significantly different results. For greater clarity, we present the results on the graph (Fig. 3.5)

Rice. 3.5 Cubic spline interpolation

If the function is specified discretely, then data matrices are specified for interpolation.

In global interpolation, polynomial interpolation is most often used n th degree or Lagrange interpolation.

The classical approach is based on the requirement of strict matching of values f(X) and j(X) at points x i(i = 0, 1, 2, … n).

We will look for the interpolation function j(X) as a degree polynomial n.

This polynomial has n+ 1 coefficient. It is natural to assume that n+ 1 conditions

j(x 0) = y 0 , j(x 1) = y 1 , . . ., j(x n) = y n (3.4)

superimposed on the polynomial

make it possible to uniquely determine its coefficients. Indeed, requiring j(X) fulfillment of conditions (3.4) , we get the system n+ 1 equations with n+ 1 unknown:

(3.6)

Solving this system for the unknowns a 0 , a 1 , …, an we obtain an analytic expression for the polynomial (3.5). System (3.6) always has a unique solution , because its determinant

known in algebra as Vandermonde determinant, different from zero . this implies , what is the interpolation polynomial j(X) for the function f(X) given in a table exists and is unique.

The resulting curve equation passes exactly through the given points. Outside the interpolation nodes, the mathematical model can have a significant error

Lagrange interpolation formula

Let the values of some function be known f(X) in n+ 1 different arbitrary points y i = f(x i) , i = 0,…, P. To interpolate (restore) a function at some point X, belonging to the segment [ x 0 ,x p], it is necessary to construct an interpolation polynomial of the nth order, which in the Lagrange method is represented as follows:

And it is easy to see that Qj(x i) = 0, if i¹ j, and Qj(x i) =1, if i= j. If we expand the product of all brackets in the numerator (in the denominator all brackets are numbers), then we get a polynomial of the nth order from X, since the numerator contains n factors of the first order. Therefore, the Lagrange interpolation polynomial is nothing more than an ordinary polynomial of the nth order, despite the specific notation.

Estimate interpolation error at a point X from [ x 0 ,xn] (i.e. solve the second

interpolation problem) can be given by the formula

In the formula - the maximum value of the (n+1)-th derivative of the original function f(X) on the segment [ x 0 ,xn]. Therefore, in order to estimate the interpolation error, some additional information about the original function is needed (this should be clear, since an infinite number of different functions can pass through the given initial points, for which the error will be different). Such information is the derivative of n + 1 order, which is not so easy to find. Below will be shown how to get out of this situation. We also note that the application of the error formula is possible only if the function is differentiable n + 1 times.

For building Lagrange interpolation formula in MathCAD it is convenient to use the function if.

if (cond, x, y)

Returns the value of x if cond is not 0 (true). Returns y if cond is 0 (false) (figure 3.6).

Interpolation formulas of Lagrange, Newton and Stirling, etc. when using a large number of interpolation nodes on the entire segment [ a, b] often lead to a poor approximation due to the accumulation of errors in the calculation process. In addition, due to the divergence of the interpolation process, increasing the number of nodes does not necessarily lead to an increase in accuracy. To reduce errors, the entire segment [ a, b] is divided into partial segments and on each of them the function is replaced by an approximately polynomial of low degree. It is called piecewise polynomial interpolation.

One of the methods of interpolation on the entire segment [ a, b] is spline interpolation.

Spline is called a piecewise polynomial function defined on the interval [ a, b] and having a certain number of continuous derivatives on this interval. The advantages of spline interpolation over conventional interpolation methods are in the convergence and stability of the computational process.

Consider one of the most common cases in practice - interpolation of the function cubic spline.
Let on the segment [ a, b] is a continuous function. Let's introduce a partition of the segment:

and denote , .

A spline corresponding to a given function and interpolation nodes (6) is a function that satisfies the following conditions:

1) on each segment , the function is a cubic polynomial;

2) the function , as well as its first and second derivatives are continuous on the segment [ a, b] ;

The third condition is called interpolation condition. A spline defined by conditions 1) - 3) is called interpolating cubic spline.

Consider a method for constructing a cubic spline.

On each of the segments, we will look for a spline function in the form of a polynomial of the third degree:

(7)

where desired coefficients.

We differentiate (7) three times with respect to X:

whence it follows

From the interpolation condition 3) we obtain:

It follows from the continuity conditions of the function.

A person can recognize his abilities only by trying to apply them. (Seneca)

Interpolation by splines: an example of building a spline in a program STATISTICS

Data structure

Let us be given the values of an unknown function on some set of points. (In fact, the variable y are the function values y=sinx on points from the segment .)

We construct an interpolation curve from these data using the program STATISTICS.

Step 1 Let's choose 2M Graphs - Scatterplot on the menu Graphic arts.

Step 2 Let's open a tab Additionally, choose as variables x and y, as a fit - Splines.

Press the OK button, and the constructed scatterplot will appear on the screen, on which the blue marks indicate the initial values between which the interpolating curve is drawn.

Let's change the number of points.

Now we have a set of twenty points as initial data.

By repeating the steps above, we get:

Let's also try to build a spline on a set of fifty points.

Fragment of the table of initial data:

Result:

And finally, let's try to build a spline using points randomly thrown onto a segment.

Initial data (table fragment):

A graph constructed in a similar way:

Now let's compare the obtained results with the original function y=sinx, whose graph looks like this:

As you can see, splines quite accurately interpolate the original function.

It can be noted that if the original function oscillates strongly, then the number of points should be large - of the order of the number of periods, but in practice such cases are rare.

Real Example: Clinical Drug Trial

Let's return to the real-life example of using splines in clinical drug trials, which was already mentioned at the very beginning.

A very important characteristic of a medicinal product is the so-called. AUC (Area under the plasma drug concentration-time curve) - area under the concentration-time curve.

This curve reflects the actual effect of the drug on the human body after the introduction of a certain dose. The AUC value is measured in mg h / l. The area under the curve depends on the rate of elimination of the drug from the body and the dose administered. The total amount of the drug removed from the body can be calculated by summing or integrating the amount of the drug removed at each individual point in time.

The AUC value is directly proportional to the dose of the administered drug for drugs with linear pharmacokinetics and inversely proportional to the so-called. clearance of the drug. The greater the clearance, the less time the drug is in the circulatory system and the faster its plasma concentration decreases. In this case, the effect of the drug on the body and the area under the concentration-time curve is less.

In clinical studies, the dependence of drug concentration in the blood over time can be determined by measuring the concentration at individual points in time. Then, the concentration graph is plotted and the AUC is estimated.

Often, the trapezoidal method is used to estimate AUC: the area under the concentration-time graph is divided into trapezoids and AUC is calculated by summing the areas of these trapezoids (which is essentially equivalent to interpolation by linear functions).

AUC=AUC0-2+AUC2-4+AUC4-6+AUC6-8+AUC8-10+AUC10-12+AUClast-infinity

In the same article, we will give an example of a more accurate AUC estimate obtained when the concentration function is interpolated by cubic splines.

Let there be concentration data obtained during the study:

Let's build a scatterplot on them and interpolate the values using a spline in the program STATISTICS.

As can be seen from the graph, the maximum concentration value C pmax = 29.78 mg/l corresponds to the time t max = 8 hours. Let's use the graph data editor and get the fitting values: