Line data Source code
1 : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : /*
3 : * This file is part of the LibreOffice project.
4 : *
5 : * This Source Code Form is subject to the terms of the Mozilla Public
6 : * License, v. 2.0. If a copy of the MPL was not distributed with this
7 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 : *
9 : * This file incorporates work covered by the following license notice:
10 : *
11 : * Licensed to the Apache Software Foundation (ASF) under one or more
12 : * contributor license agreements. See the NOTICE file distributed
13 : * with this work for additional information regarding copyright
14 : * ownership. The ASF licenses this file to you under the Apache
15 : * License, Version 2.0 (the "License"); you may not use this file
16 : * except in compliance with the License. You may obtain a copy of
17 : * the License at http://www.apache.org/licenses/LICENSE-2.0 .
18 : */
19 :
20 : // Natural, Clamped, or Periodic Cubic Splines
21 : //
22 : // Input: A list of N+1 points (x_i,a_i), 0 <= i <= N, which are sampled
23 : // from a function, a_i = f(x_i). The function f is unknown. Boundary
24 : // conditions are
25 : // (1) Natural splines: f"(x_0) = f"(x_N) = 0
26 : // (2) Clamped splines: f'(x_0) and f'(x_N) are user-specified.
27 : // (3) Periodic splines: f(x_0) = f(x_N) [in which case a_N = a_0 is
28 : // required in the input], f'(x_0) = f'(x_N), and f"(x_0) = f"(x_N).
29 : //
30 : // Output: b_i, c_i, d_i, 0 <= i <= N-1, which are coefficients for the cubic
31 : // spline S_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3 for
32 : // x_i <= x < x_{i+1}.
33 : //
34 : // The natural and clamped algorithms were implemented from
35 : //
36 : // Numerical Analysis, 3rd edition
37 : // Richard L. Burden and J. Douglas Faires
38 : // Prindle, Weber & Schmidt
39 : // Boston, 1985, pp. 122-124.
40 : //
41 : // The algorithm sets up a tridiagonal linear system of equations in the
42 : // c_i values. This can be solved in O(N) time.
43 : //
44 : // The periodic spline algorithm was implemented from my own derivation. The
45 : // linear system of equations is not tridiagonal. For now I use a standard
46 : // linear solver that does not take advantage of the sparseness of the
47 : // matrix. Therefore for very large N, you may have to worry about memory
48 : // usage.
49 :
50 : #include <sal/config.h>
51 :
52 : #include "cspline.h"
53 : #include "solver.h"
54 :
55 0 : void NaturalSpline (int N, double* x, double* a, double*& b, double*& c,
56 : double*& d)
57 : {
58 0 : const double oneThird = 1.0/3.0;
59 :
60 : int i;
61 0 : double* h = new double[N];
62 0 : double* hdiff = new double[N];
63 0 : double* alpha = new double[N];
64 :
65 0 : for (i = 0; i < N; i++){
66 0 : h[i] = x[i+1]-x[i];
67 : }
68 :
69 0 : for (i = 1; i < N; i++)
70 0 : hdiff[i] = x[i+1]-x[i-1];
71 :
72 0 : for (i = 1; i < N; i++)
73 : {
74 0 : double numer = 3.0*(a[i+1]*h[i-1]-a[i]*hdiff[i]+a[i-1]*h[i]);
75 0 : double denom = h[i-1]*h[i];
76 0 : alpha[i] = numer/denom;
77 : }
78 :
79 0 : double* ell = new double[N+1];
80 0 : double* mu = new double[N];
81 0 : double* z = new double[N+1];
82 : double recip;
83 :
84 0 : ell[0] = 1.0;
85 0 : mu[0] = 0.0;
86 0 : z[0] = 0.0;
87 :
88 0 : for (i = 1; i < N; i++)
89 : {
90 0 : ell[i] = 2.0*hdiff[i]-h[i-1]*mu[i-1];
91 0 : recip = 1.0/ell[i];
92 0 : mu[i] = recip*h[i];
93 0 : z[i] = recip*(alpha[i]-h[i-1]*z[i-1]);
94 : }
95 0 : ell[N] = 1.0;
96 0 : z[N] = 0.0;
97 :
98 0 : b = new double[N];
99 0 : c = new double[N+1];
100 0 : d = new double[N];
101 :
102 0 : c[N] = 0.0;
103 :
104 0 : for (i = N-1; i >= 0; i--)
105 : {
106 0 : c[i] = z[i]-mu[i]*c[i+1];
107 0 : recip = 1.0/h[i];
108 0 : b[i] = recip*(a[i+1]-a[i])-h[i]*(c[i+1]+2.0*c[i])*oneThird;
109 0 : d[i] = oneThird*recip*(c[i+1]-c[i]);
110 : }
111 :
112 0 : delete[] h;
113 0 : delete[] hdiff;
114 0 : delete[] alpha;
115 0 : delete[] ell;
116 0 : delete[] mu;
117 0 : delete[] z;
118 0 : }
119 :
120 0 : void PeriodicSpline (int N, double* x, double* a, double*& b, double*& c,
121 : double*& d)
122 : {
123 0 : double* h = new double[N];
124 : int i;
125 0 : for (i = 0; i < N; i++)
126 0 : h[i] = x[i+1]-x[i];
127 :
128 0 : mgcLinearSystemD sys;
129 0 : double** mat = mgcLinearSystemD::NewMatrix(N+1); // guaranteed to be zeroed memory
130 0 : c = mgcLinearSystemD::NewVector(N+1); // guaranteed to be zeroed memory
131 :
132 : // c[0] - c[N] = 0
133 0 : mat[0][0] = +1.0f;
134 0 : mat[0][N] = -1.0f;
135 :
136 : // h[i-1]*c[i-1]+2*(h[i-1]+h[i])*c[i]+h[i]*c[i+1] =
137 : // 3*{(a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]}
138 0 : for (i = 1; i <= N-1; i++)
139 : {
140 0 : mat[i][i-1] = h[i-1];
141 0 : mat[i][i ] = 2.0f*(h[i-1]+h[i]);
142 0 : mat[i][i+1] = h[i];
143 0 : c[i] = 3.0f*((a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]);
144 : }
145 :
146 : // "wrap around equation" for periodicity
147 : // h[N-1]*c[N-1]+2*(h[N-1]+h[0])*c[0]+h[0]*c[1] =
148 : // 3*{(a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]}
149 0 : mat[N][N-1] = h[N-1];
150 0 : mat[N][0 ] = 2.0f*(h[N-1]+h[0]);
151 0 : mat[N][1 ] = h[0];
152 0 : c[N] = 3.0f*((a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]);
153 :
154 : // solve for c[0] through c[N]
155 0 : mgcLinearSystemD::Solve(N+1,mat,c);
156 :
157 0 : const double oneThird = 1.0/3.0;
158 0 : b = new double[N];
159 0 : d = new double[N];
160 0 : for (i = 0; i < N; i++)
161 : {
162 0 : b[i] = (a[i+1]-a[i])/h[i] - oneThird*(c[i+1]+2.0f*c[i])*h[i];
163 0 : d[i] = oneThird*(c[i+1]-c[i])/h[i];
164 : }
165 :
166 0 : delete[] h;
167 0 : mgcLinearSystemD::DeleteMatrix(N+1,mat);
168 0 : }
169 :
170 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
