1 /*************************************************************************
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
5 * Copyright 2008 by Sun Microsystems, Inc.
7 * OpenOffice.org - a multi-platform office productivity suite
9 * $RCSfile: cspline.cpp,v $
12 * This file is part of OpenOffice.org.
14 * OpenOffice.org is free software: you can redistribute it and/or modify
15 * it under the terms of the GNU Lesser General Public License version 3
16 * only, as published by the Free Software Foundation.
18 * OpenOffice.org is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU Lesser General Public License version 3 for more details
22 * (a copy is included in the LICENSE file that accompanied this code).
24 * You should have received a copy of the GNU Lesser General Public License
25 * version 3 along with OpenOffice.org. If not, see
26 * <http://www.openoffice.org/license.html>
27 * for a copy of the LGPLv3 License.
29 ************************************************************************/
31 // Natural, Clamped, or Periodic Cubic Splines
33 // Input: A list of N+1 points (x_i,a_i), 0 <= i <= N, which are sampled
34 // from a function, a_i = f(x_i). The function f is unknown. Boundary
36 // (1) Natural splines: f"(x_0) = f"(x_N) = 0
37 // (2) Clamped splines: f'(x_0) and f'(x_N) are user-specified.
38 // (3) Periodic splines: f(x_0) = f(x_N) [in which case a_N = a_0 is
39 // required in the input], f'(x_0) = f'(x_N), and f"(x_0) = f"(x_N).
41 // Output: b_i, c_i, d_i, 0 <= i <= N-1, which are coefficients for the cubic
42 // spline S_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3 for
43 // x_i <= x < x_{i+1}.
45 // The natural and clamped algorithms were implemented from
47 // Numerical Analysis, 3rd edition
48 // Richard L. Burden and J. Douglas Faires
49 // Prindle, Weber & Schmidt
50 // Boston, 1985, pp. 122-124.
52 // The algorithm sets up a tridiagonal linear system of equations in the
53 // c_i values. This can be solved in O(N) time.
55 // The periodic spline algorithm was implemented from my own derivation. The
56 // linear system of equations is not tridiagonal. For now I use a standard
57 // linear solver that does not take advantage of the sparseness of the
58 // matrix. Therefore for very large N, you may have to worry about memory
62 //-----------------------------------------------------------------------------
63 void NaturalSpline (int N
, double* x
, double* a
, double*& b
, double*& c
,
66 const double oneThird
= 1.0/3.0;
69 double* h
= new double[N
];
70 double* hdiff
= new double[N
];
71 double* alpha
= new double[N
];
73 for (i
= 0; i
< N
; i
++){
77 for (i
= 1; i
< N
; i
++)
78 hdiff
[i
] = x
[i
+1]-x
[i
-1];
80 for (i
= 1; i
< N
; i
++)
82 double numer
= 3.0*(a
[i
+1]*h
[i
-1]-a
[i
]*hdiff
[i
]+a
[i
-1]*h
[i
]);
83 double denom
= h
[i
-1]*h
[i
];
84 alpha
[i
] = numer
/denom
;
87 double* ell
= new double[N
+1];
88 double* mu
= new double[N
];
89 double* z
= new double[N
+1];
96 for (i
= 1; i
< N
; i
++)
98 ell
[i
] = 2.0*hdiff
[i
]-h
[i
-1]*mu
[i
-1];
101 z
[i
] = recip
*(alpha
[i
]-h
[i
-1]*z
[i
-1]);
112 for (i
= N
-1; i
>= 0; i
--)
114 c
[i
] = z
[i
]-mu
[i
]*c
[i
+1];
116 b
[i
] = recip
*(a
[i
+1]-a
[i
])-h
[i
]*(c
[i
+1]+2.0*c
[i
])*oneThird
;
117 d
[i
] = oneThird
*recip
*(c
[i
+1]-c
[i
]);
128 void PeriodicSpline (int N
, double* x
, double* a
, double*& b
, double*& c
,
131 double* h
= new double[N
];
133 for (i
= 0; i
< N
; i
++)
136 mgcLinearSystemD sys
;
137 double** mat
= sys
.NewMatrix(N
+1); // guaranteed to be zeroed memory
138 c
= sys
.NewVector(N
+1); // guaranteed to be zeroed memory
144 // h[i-1]*c[i-1]+2*(h[i-1]+h[i])*c[i]+h[i]*c[i+1] =
145 // 3*{(a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]}
146 for (i
= 1; i
<= N
-1; i
++)
148 mat
[i
][i
-1] = h
[i
-1];
149 mat
[i
][i
] = 2.0f
*(h
[i
-1]+h
[i
]);
151 c
[i
] = 3.0f
*((a
[i
+1]-a
[i
])/h
[i
] - (a
[i
]-a
[i
-1])/h
[i
-1]);
154 // "wrap around equation" for periodicity
155 // h[N-1]*c[N-1]+2*(h[N-1]+h[0])*c[0]+h[0]*c[1] =
156 // 3*{(a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]}
157 mat
[N
][N
-1] = h
[N
-1];
158 mat
[N
][0 ] = 2.0f
*(h
[N
-1]+h
[0]);
160 c
[N
] = 3.0f
*((a
[1]-a
[0])/h
[0] - (a
[0]-a
[N
-1])/h
[N
-1]);
162 // solve for c[0] through c[N]
163 sys
.Solve(N
+1,mat
,c
);
165 const double oneThird
= 1.0/3.0;
168 for (i
= 0; i
< N
; i
++)
170 b
[i
] = (a
[i
+1]-a
[i
])/h
[i
] - oneThird
*(c
[i
+1]+2.0f
*c
[i
])*h
[i
];
171 d
[i
] = oneThird
*(c
[i
+1]-c
[i
])/h
[i
];
174 sys
.DeleteMatrix(N
+1,mat
);