8
double *d3_np_fs ( int n, double a[], double b[] );
9
double *spline_cubic_set ( int n, double t[], double y[], int ibcbeg, double ybcbeg, int ibcend, double ybcend );
10
double spline_cubic_val ( int n, double t[], double tval, double y[], double ypp[], double *ypval, double *yppval );
17
CubicSpline(const CubicSpline& Other);
18
CubicSpline& operator = (const CubicSpline& Other);
20
CubicSpline(int numpoint, std::vector<double> x_array, std::vector<double> y_array, int ibcbeg = 0, double ybcbeg = 0, int ibcend = 0, double ybcend = 0);
21
void Set(int numpoint, std::vector<double> x_array, std::vector<double> y_array, int ibcbeg = 0, double ybcbeg = 0, int ibcend = 0, double ybcend = 0);
25
//**********************************************************************
26
// double* CubicSpline::Compute ( int n, double t[], double y[], int ibcbeg, double ybcbeg, int ibcend, double ybcend )
29
// SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic spline.
33
// For data interpolation, the user must call SPLINE_SET to determine
34
// the second derivative data, passing in the data to be interpolated,
35
// and the desired boundary conditions.
37
// The data to be interpolated, plus the SPLINE_SET output, defines
38
// the spline. The user may then call SPLINE_VAL to evaluate the
39
// spline at any point.
41
// The cubic spline is a piecewise cubic polynomial. The intervals
42
// are determined by the "knots" or abscissas of the data to be
43
// interpolated. The cubic spline has continous first and second
44
// derivatives over the entire interval of interpolation.
46
// For any point T in the interval T(IVAL), T(IVAL+1), the form of
50
// + B(IVAL) * ( T - T(IVAL) )
51
// + C(IVAL) * ( T - T(IVAL) )**2
52
// + D(IVAL) * ( T - T(IVAL) )**3
54
// If we assume that we know the values Y(*) and DDY(*), which represent
55
// the values and second derivatives of the spline at each knot, then
56
// the coefficients can be computed as:
59
// B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
60
// - ( DDY(IVAL+1) + 2 * DDY(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
61
// C(IVAL) = DDY(IVAL) / 2
62
// D(IVAL) = ( DDY(IVAL+1) - DDY(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
64
// Since the first derivative of the spline is
67
// + 2 * C(IVAL) * ( T - T(IVAL) )
68
// + 3 * D(IVAL) * ( T - T(IVAL) )**2,
70
// the requirement that the first derivative be continuous at interior
71
// knot I results in a total of N-2 equations, of the form:
73
// B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1)) + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))^2 = B(IVAL)
75
// or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
77
// ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) - ( DDY(IVAL) + 2 * DDY(IVAL-1) ) * H(IVAL-1) / 6 + DDY(IVAL-1) * H(IVAL-1)
78
// + ( DDY(IVAL) - DDY(IVAL-1) ) * H(IVAL-1) / 2
79
// = ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) - ( DDY(IVAL+1) + 2 * DDY(IVAL) ) * H(IVAL) / 6
83
// DDY(IVAL-1) * H(IVAL-1) + 2 * DDY(IVAL) * ( H(IVAL-1) + H(IVAL) ) + DDY(IVAL) * H(IVAL)
84
// = 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
86
// Boundary conditions must be applied at the first and last knots.
87
// The resulting tridiagonal system can be solved for the DDY values.
92
// Input, int N, the number of data points. N must be at least 2.
93
// In the special case where N = 2 and IBCBEG = IBCEND = 0, the
94
// spline will actually be linear.
96
// Input, double T[N], the knot values, that is, the points were data is
97
// specified. The knot values should be distinct, and increasing.
99
// Input, double Y[N], the data values to be interpolated.
101
// Input, int IBCBEG, left boundary condition flag:
102
// 0: the cubic spline should be a quadratic over the first interval;
103
// 1: the first derivative at the left endpoint should be YBCBEG;
104
// 2: the second derivative at the left endpoint should be YBCBEG.
106
// Input, double YBCBEG, the values to be used in the boundary
107
// conditions if IBCBEG is equal to 1 or 2.
109
// Input, int IBCEND, right boundary condition flag:
110
// 0: the cubic spline should be a quadratic over the last interval;
111
// 1: the first derivative at the right endpoint should be YBCEND;
112
// 2: the second derivative at the right endpoint should be YBCEND.
114
// Input, double YBCEND, the values to be used in the boundary
115
// conditions if IBCEND is equal to 1 or 2.
117
// Output, double SPLINE_CUBIC_SET[N], the second derivatives of the cubic spline.
119
double* Compute(int ibcbeg, double ybcbeg, int ibcend, double ybcend);
121
//**********************************************************************
122
// double Eval ( int n, double t[], double tval, double y[], double ddy[], double *dyval, double *ddyval )
125
// SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point.
129
// SPLINE_CUBIC_SET must have already been called to define the values of YPP.
131
// For any point T in the interval T(IVAL), T(IVAL+1), the form of
135
// + B * ( T - T(IVAL) )
136
// + C * ( T - T(IVAL) )**2
137
// + D * ( T - T(IVAL) )**3
141
// B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
142
// - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
144
// D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
156
// Input, int N, the number of knots.
158
// Input, double Y[N], the data values at the knots.
160
// Input, double T[N], the knot values.
162
// Input, double TVAL, a point, typically between T[0] and T[N-1], at
163
// which the spline is to be evalulated. If TVAL lies outside
164
// this range, extrapolation is used.
166
// Input, double Y[N], the data values at the knots.
168
// Input, double YPP[N], the second derivatives of the spline at
171
// Output, double *YPVAL, the derivative of the spline at TVAL.
173
// Output, double *YPPVAL, the second derivative of the spline at TVAL.
175
// Output, double SPLINE_VAL, the value of the spline at TVAL.
177
double Eval (double tval);
179
//**********************************************************************
180
// double* SolveTridiag ( int n, double a[], double b[] );
183
// D3_NP_FS factors and solves a D3 system.
187
// The D3 storage format is used for a tridiagonal matrix.
188
// The superdiagonal is stored in entries (1,2:N), the diagonal in
189
// entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
190
// original matrix is "collapsed" vertically into the array.
192
// This algorithm requires that each diagonal entry be nonzero.
193
// It does not use pivoting, and so can fail on systems that
194
// are actually nonsingular.
198
// Here is how a D3 matrix of order 5 would be stored:
201
// A11 A22 A33 A44 A55
206
// Input, int N, the order of the linear system.
208
// Input/output, double A[3*N].
209
// On input, the nonzero diagonals of the linear system.
210
// On output, the data in these vectors has been overwritten
211
// by factorization information.
213
// Input, double B[N], the right hand side.
215
// Output, double D3_NP_FS[N], the solution of the linear system.
216
// This is NULL if there was an error because one of the diagonal
219
double* SolveTridiag ( int n, double a[], double b[] );
231
int np; // number of points
added
removed
NuxCore/Math/Spline.h
