c3-3 (779473)
Текст из файла
3.3 Cubic Spline Interpolation113}CITED REFERENCES AND FURTHER READING:Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§2.2. [1]Gear, C.W. 1971, Numerical Initial Value Problems in Ordinary Differential Equations (EnglewoodCliffs, NJ: Prentice-Hall), §6.2.Cuyt, A., and Wuytack, L. 1987, Nonlinear Methods in Numerical Analysis (Amsterdam: NorthHolland), Chapter 3.3.3 Cubic Spline InterpolationGiven a tabulated function yi = y(xi ), i = 1...N , focus attention on oneparticular interval, between xj and xj+1 . Linear interpolation in that interval givesthe interpolation formulay = Ayj + Byj+1(3.3.1)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).c=vector(1,n);d=vector(1,n);hh=fabs(x-xa[1]);for (i=1;i<=n;i++) {h=fabs(x-xa[i]);if (h == 0.0) {*y=ya[i];*dy=0.0;FREERETURN} else if (h < hh) {ns=i;hh=h;}c[i]=ya[i];d[i]=ya[i]+TINY;The TINY part is needed to prevent a rare zero-over-zero}condition.*y=ya[ns--];for (m=1;m<n;m++) {for (i=1;i<=n-m;i++) {w=c[i+1]-d[i];h=xa[i+m]-x;h will never be zero, since this was tested in the initialt=(xa[i]-x)*d[i]/h;izing loop.dd=t-c[i+1];if (dd == 0.0) nrerror("Error in routine ratint");This error condition indicates that the interpolating function has a pole at therequested value of x.dd=w/dd;d[i]=c[i+1]*dd;c[i]=t*dd;}*y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));}FREERETURN114Chapter 3.Interpolation and ExtrapolationwhereA≡xj+1 − xxj+1 − xjB ≡ 1−A =x − xjxj+1 − xj(3.3.2)00y = Ayj + Byj+1 + Cyj00 + Dyj+1(3.3.3)where A and B are defined in (3.3.2) andC≡1 3(A − A)(xj+1 − xj )26D≡1 3(B − B)(xj+1 − xj )26(3.3.4)Notice that the dependence on the independent variable x in equations (3.3.3) and(3.3.4) is entirely through the linear x-dependence of A and B, and (through A andB) the cubic x-dependence of C and D.We can readily check that y00 is in fact the second derivative of the newinterpolating polynomial.
We take derivatives of equation (3.3.3) with respectto x, using the definitions of A, B, C, D to compute dA/dx, dB/dx, dC/dx, anddD/dx. The result isyj+1 − yj3A2 − 13B 2 − 1dy00=(xj+1 − xj )yj00 +(xj+1 − xj )yj+1−(3.3.5)dxxj+1 − xj66for the first derivative, andd2 y00= Ayj00 + Byj+1dx2(3.3.6)for the second derivative. Since A = 1 at xj , A = 0 at xj+1 , while B is just theother way around, (3.3.6) shows that y00 is just the tabulated second derivative, andalso that the second derivative will be continuous across (e.g.) the boundary betweenthe two intervals (xj−1 , xj ) and (xj , xj+1 ).Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Equations (3.3.1) and (3.3.2) are a special case of the general Lagrange interpolationformula (3.1.1).Since it is (piecewise) linear, equation (3.3.1) has zero second derivative inthe interior of each interval, and an undefined, or infinite, second derivative at theabscissas xj .
The goal of cubic spline interpolation is to get an interpolation formulathat is smooth in the first derivative, and continuous in the second derivative, bothwithin an interval and at its boundaries.Suppose, contrary to fact, that in addition to the tabulated values of yi , wealso have tabulated values for the function’s second derivatives, y00 , that is, a setof numbers yi00 . Then, within each interval, we can add to the right-hand side ofequation (3.3.1) a cubic polynomial whose second derivative varies linearly from a00on the right. Doing so, we will have the desiredvalue yj00 on the left to a value yj+1continuous second derivative.
If we also construct the cubic polynomial to havezero values at xj and xj+1 , then adding it in will not spoil the agreement with thetabulated functional values yj and yj+1 at the endpoints xj and xj+1 .A little side calculation shows that there is only one way to arrange thisconstruction, namely replacing (3.3.1) by3.3 Cubic Spline Interpolation115xj+1 − xj−1 00 xj+1 − xj 00yj+1 − yjyj − yj−1xj − xj−1 00yj−1 +yj +yj+1 =−636xj+1 − xjxj − xj−1(3.3.7)These are N − 2 linear equations in the N unknowns yi00 , i = 1, . . .
, N . Thereforethere is a two-parameter family of possible solutions.For a unique solution, we need to specify two further conditions, typically takenas boundary conditions at x1 and xN . The most common ways of doing this are either00equal to zero, giving the so-called natural• set one or both of y100 and yNcubic spline, which has zero second derivative on one or both of itsboundaries, or00• set either of y100 and yNto values calculated from equation (3.3.5) so asto make the first derivative of the interpolating function have a specifiedvalue on either or both boundaries.One reason that cubic splines are especially practical is that the set of equations(3.3.7), along with the two additional boundary conditions, are not only linear, butalso tridiagonal.
Each yj00 is coupled only to its nearest neighbors at j ± 1. Therefore,the equations can be solved in O(N ) operations by the tridiagonal algorithm (§2.4).That algorithm is concise enough to build right into the spline calculational routine.This makes the routine not completely transparent as an implementation of (3.3.7),so we encourage you to study it carefully, comparing with tridag (§2.4). Arraysare assumed to be unit-offset. If you have zero-offset arrays, see §1.2.#include "nrutil.h"void spline(float x[], float y[], int n, float yp1, float ypn, float y2[])Given arrays x[1..n] and y[1..n] containing a tabulated function, i.e., yi = f (xi ), withx1 < x2 < .
. . < xN , and given values yp1 and ypn for the first derivative of the interpolatingfunction at points 1 and n, respectively, this routine returns an array y2[1..n] that containsthe second derivatives of the interpolating function at the tabulated points xi . If yp1 and/orypn are equal to 1 × 1030 or larger, the routine is signaled to set the corresponding boundarycondition for a natural spline, with zero second derivative on that boundary.{int i,k;float p,qn,sig,un,*u;u=vector(1,n-1);if (yp1 > 0.99e30)The lower boundary condition is set either to be “naty2[1]=u[1]=0.0;ural”else {or else to have a specified first derivative.y2[1] = -0.5;u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);}Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).The only problem now is that we supposed the yi00 ’s to be known, when, actually,they are not. However, we have not yet required that the first derivative, computedfrom equation (3.3.5), be continuous across the boundary between two intervals.
Thekey idea of a cubic spline is to require this continuity and to use it to get equationsfor the second derivatives yi00 .The required equations are obtained by setting equation (3.3.5) evaluated forx = xj in the interval (xj−1 , xj ) equal to the same equation evaluated for x = xj butin the interval (xj , xj+1). With some rearrangement, this gives (for j = 2, . . . , N −1)116Chapter 3.Interpolation and Extrapolation}It is important to understand that the program spline is called only once toprocess an entire tabulated function in arrays xi and yi . Once this has been done,values of the interpolated function for any value of x are obtained by calls (as manyas desired) to a separate routine splint (for “spline interpolation”):void splint(float xa[], float ya[], float y2a[], int n, float x, float *y)Given the arrays xa[1..n] and ya[1..n], which tabulate a function (with the xai ’s in order),and given the array y2a[1..n], which is the output from spline above, and given a value ofx, this routine returns a cubic-spline interpolated value y.{void nrerror(char error_text[]);int klo,khi,k;float h,b,a;klo=1;We will find the right place in the table by means ofkhi=n;bisection.
This is optimal if sequential calls to thiswhile (khi-klo > 1) {routine are at random values of x. If sequential callsk=(khi+klo) >> 1;are in order, and closely spaced, one would do betterif (xa[k] > x) khi=k;to store previous values of klo and khi and test ifelse klo=k;they remain appropriate on the next call.}klo and khi now bracket the input value of x.h=xa[khi]-xa[klo];if (h == 0.0) nrerror("Bad xa input to routine splint"); The xa’s must be disa=(xa[khi]-x)/h;tinct.b=(x-xa[klo])/h;Cubic spline polynomial is now evaluated.*y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;}CITED REFERENCES AND FURTHER READING:De Boor, C. 1978, A Practical Guide to Splines (New York: Springer-Verlag).Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for MathematicalComputations (Englewood Cliffs, NJ: Prentice-Hall), §§4.4–4.5.Stoer, J., and Bulirsch, R.
1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§2.4.Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York:McGraw-Hill), §3.8.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
