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 "solver.h"
51 :
52 0 : void NaturalSpline (int N, double* x, double* a, double*& b, double*& c,
53 : double*& d)
54 : {
55 0 : const double oneThird = 1.0/3.0;
56 :
57 : int i;
58 0 : double* h = new double[N];
59 0 : double* hdiff = new double[N];
60 0 : double* alpha = new double[N];
61 :
62 0 : for (i = 0; i < N; i++){
63 0 : h[i] = x[i+1]-x[i];
64 : }
65 :
66 0 : for (i = 1; i < N; i++)
67 0 : hdiff[i] = x[i+1]-x[i-1];
68 :
69 0 : for (i = 1; i < N; i++)
70 : {
71 0 : double numer = 3.0*(a[i+1]*h[i-1]-a[i]*hdiff[i]+a[i-1]*h[i]);
72 0 : double denom = h[i-1]*h[i];
73 0 : alpha[i] = numer/denom;
74 : }
75 :
76 0 : double* ell = new double[N+1];
77 0 : double* mu = new double[N];
78 0 : double* z = new double[N+1];
79 : double recip;
80 :
81 0 : ell[0] = 1.0;
82 0 : mu[0] = 0.0;
83 0 : z[0] = 0.0;
84 :
85 0 : for (i = 1; i < N; i++)
86 : {
87 0 : ell[i] = 2.0*hdiff[i]-h[i-1]*mu[i-1];
88 0 : recip = 1.0/ell[i];
89 0 : mu[i] = recip*h[i];
90 0 : z[i] = recip*(alpha[i]-h[i-1]*z[i-1]);
91 : }
92 0 : ell[N] = 1.0;
93 0 : z[N] = 0.0;
94 :
95 0 : b = new double[N];
96 0 : c = new double[N+1];
97 0 : d = new double[N];
98 :
99 0 : c[N] = 0.0;
100 :
101 0 : for (i = N-1; i >= 0; i--)
102 : {
103 0 : c[i] = z[i]-mu[i]*c[i+1];
104 0 : recip = 1.0/h[i];
105 0 : b[i] = recip*(a[i+1]-a[i])-h[i]*(c[i+1]+2.0*c[i])*oneThird;
106 0 : d[i] = oneThird*recip*(c[i+1]-c[i]);
107 : }
108 :
109 0 : delete[] h;
110 0 : delete[] hdiff;
111 0 : delete[] alpha;
112 0 : delete[] ell;
113 0 : delete[] mu;
114 0 : delete[] z;
115 0 : }
116 :
117 0 : void PeriodicSpline (int N, double* x, double* a, double*& b, double*& c,
118 : double*& d)
119 : {
120 0 : double* h = new double[N];
121 : int i;
122 0 : for (i = 0; i < N; i++)
123 0 : h[i] = x[i+1]-x[i];
124 :
125 0 : mgcLinearSystemD sys;
126 0 : double** mat = sys.NewMatrix(N+1); // guaranteed to be zeroed memory
127 0 : c = sys.NewVector(N+1); // guaranteed to be zeroed memory
128 :
129 : // c[0] - c[N] = 0
130 0 : mat[0][0] = +1.0f;
131 0 : mat[0][N] = -1.0f;
132 :
133 : // h[i-1]*c[i-1]+2*(h[i-1]+h[i])*c[i]+h[i]*c[i+1] =
134 : // 3*{(a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]}
135 0 : for (i = 1; i <= N-1; i++)
136 : {
137 0 : mat[i][i-1] = h[i-1];
138 0 : mat[i][i ] = 2.0f*(h[i-1]+h[i]);
139 0 : mat[i][i+1] = h[i];
140 0 : c[i] = 3.0f*((a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]);
141 : }
142 :
143 : // "wrap around equation" for periodicity
144 : // h[N-1]*c[N-1]+2*(h[N-1]+h[0])*c[0]+h[0]*c[1] =
145 : // 3*{(a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]}
146 0 : mat[N][N-1] = h[N-1];
147 0 : mat[N][0 ] = 2.0f*(h[N-1]+h[0]);
148 0 : mat[N][1 ] = h[0];
149 0 : c[N] = 3.0f*((a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]);
150 :
151 : // solve for c[0] through c[N]
152 0 : sys.Solve(N+1,mat,c);
153 :
154 0 : const double oneThird = 1.0/3.0;
155 0 : b = new double[N];
156 0 : d = new double[N];
157 0 : for (i = 0; i < N; i++)
158 : {
159 0 : b[i] = (a[i+1]-a[i])/h[i] - oneThird*(c[i+1]+2.0f*c[i])*h[i];
160 0 : d[i] = oneThird*(c[i+1]-c[i])/h[i];
161 : }
162 :
163 0 : delete[] h;
164 0 : sys.DeleteMatrix(N+1,mat);
165 0 : }
166 :
167 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
