PLUTO
Functions
math_interp.c File Reference

Miscellaneous functions for handling interpolation. More...

#include "pluto.h"
Include dependency graph for math_interp.c:

Go to the source code of this file.

Functions

void MonotoneSplineCoeffs (double *x, double *y, double *dydx, int n, double *a, double *b, double *c, double *d)
 
void SplineCoeffs (double *x, double *f, double dfL, double dfR, int n, double *a, double *b, double *c, double *d)
 

Detailed Description

Miscellaneous functions for handling interpolation.

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Jan 2, 2014

Definition in file math_interp.c.

Function Documentation

void MonotoneSplineCoeffs ( double *  x,
double *  y,
double *  dydx,
int  n,
double *  a,
double *  b,
double *  c,
double *  d 
)

Compute the cubic spline coefficients for interpolating a monotonic dataset known at node values f[k] = f(x[k]) and derivative dfdx(x[k]). The resulting spline will be the smoothest curve that passes through the control points while preserving monotonocity of the dataset but not necessarily continuity of the second derivative. Coefficients will be computed in each interval x[k] < x < x[k+1].

The x[k] may not be equally spaced.

Parameters
[in]x1D array of abscissas
[in]f1D array of function values
[in]dfdx1D array of function derivative
[in]nthe number of function values.
[out]athe x^3 cubic coefficient
[out]bthe x^2 cubic coefficient
[out]cthe x cubic coefficient
[out]dthe (1) cubic coefficient

Reference:

  • "Monotonic cubic spline interpolation", G. Wolberg, I. Alfy, Proceedings of Computer Graphics International (1999)
  • "An energy-minimization framework for monotonic cubic spline interpolation", Wolberg & Alfy, J. of comput. and applied math. (2002) 143, 145

Definition at line 13 of file math_interp.c.

44 {
45  int k, nfail = 0;
46  char cm, c1, c2, c2a, c2b, c2c;
47  double m, alpha, beta, dx;
48 
49  for (k = 0; k < n-1; k++){
50  m = (y[k+1] - y[k])/(x[k+1] - x[k]);
51  alpha = dydx[k]/m;
52  beta = dydx[k+1]/m;
53 
54  /* -----------------------------------------------------
55  Check monotonicity conditions. Quit if not verified
56  ----------------------------------------------------- */
57 
58  cm = (DSIGN(dydx[k]) == DSIGN(dydx[k+1])) && (DSIGN(dydx[k]) == DSIGN(m));
59  c1 = (alpha + beta - 2.0) <= 0.0;
60  c2 = (alpha + 2.0*beta - 2.0) > 0.0;
61  c2a = (2.0*alpha + beta - 3.0) <= 0.0;
62  c2b = (alpha + 2.0*beta - 3.0) <= 0.0;
63  c2c = (alpha*alpha + alpha*(beta-6.0) + (beta-3.0)*(beta-3.0)) < 0.0;
64  if ( (cm && c1) || (c2 && cm && (c2a || c2b || c2c))){
65  dx = x[k+1] - x[k];
66  a[k] = dx*(-2.0*m + dydx[k] + dydx[k+1]);
67  b[k] = dx*( 3.0*m - 2.0*dydx[k] - dydx[k+1]);
68  c[k] = dx*dydx[k];
69  d[k] = y[k];
70  }else{
71  print1 ("! MonotoneSplineCoeffs(): monotonicity condition not ");
72  print1 ("satisifed in SPLINE1 \n");
73  print1 ("cm = %d, c1 = %d, c2 = %d, c2a = %d\n",cm,c1,c2,c2a);
74  QUIT_PLUTO(1);
75  }
76 
77  /* ---------------------------------------------------------
78  Reset a = b = 0 if the interpolant does not have
79  enough curvature.
80  This is done by checking the size of |y''/y| at x = 0.5
81  (midpoint).
82  --------------------------------------------------------- */
83 
84  alpha = 6.0*a[k]*0.5 + 2.0*b[k];
85  beta = d[k] + 0.5*(c[k] + 0.5*(b[k] + 0.5*a[k]));
86  m = fabs(alpha/beta);
87  if (m < 1.e-16){
88  nfail++;
89  a[k] = b[k] = 0.0;
90  }
91  }
92 
93 }
alpha
static double alpha
Definition: init.c:3
a
static double a
Definition: init.c:135
n
static int n
Definition: analysis.c:3
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
DSIGN
#define DSIGN(x)
Definition: macros.h:110
k
int k
Definition: analysis.c:2
menu.c
tuple c
Definition: menu.py:375
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

Here is the caller graph for this function:

void SplineCoeffs ( double *  x,
double *  f,
double  dfL,
double  dfR,
int  n,
double *  a,
double *  b,
double *  c,
double *  d 
)

Compute the cubic spline coefficients for interpolating a dataset (not necessarily monotonic) known at node values f[k] = f(x[k]). The resulting spline will be continuous up to the second derivative. Coefficients will be computed in each interval x[k] < x < x[k+1].

The x[k] may not be equally spaced.

Parameters
[in]x1D array of abscissas
[in]f1D array of function values
[in]dfLderivative at the leftmost node
[in]dfRderivative at the rightmost node
[in]nthe number of function values.
[out]athe x^3 cubic coefficient
[out]bthe x^2 cubic coefficient
[out]cthe x cubic coefficient
[out]dthe (1) cubic coefficient

Reference:

  • "Numerical Recipe in C, second edition", Section 3.3.

Definition at line 95 of file math_interp.c.

120 {
121  int i, nfail=0;
122  double p, sig, qn, un, dx;
123  double *u, *d2f;
124  double xm, xp, det, aa,bb,cc; /* coefficients of the derivative (parabola) */
125 
126  if (u == NULL) {
127  u = ARRAY_1D(n, double);
128  d2f = ARRAY_1D(n, double);
129  }
130 
131  d2f[0] = -0.5;
132  u[0] = (3.0/(x[1] - x[0])*(f[1] - f[0])/(x[1]-x[0]) - dfL);
133 
134 /* ---------------------------------------------------------------
135  Start decomposition loop of the tridiagonal algorithm.
136  d2f[] and u[] are used for temporary storage of the
137  decomposed factors.
138  --------------------------------------------------------------- */
139 
140  for (i = 1; i <= n-2;i++) {
141  sig = (x[i] - x[i-1])/(x[i+1] - x[i-1]);
142  p = sig*d2f[i-1] + 2.0;
143  d2f[i] = (sig - 1.0)/p;
144  u[i] = (f[i+1] - f[i])/(x[i+1]-x[i]) - (f[i] - f[i-1])/(x[i]-x[i-1]);
145  u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
146  }
147  qn = 0.5;
148  un = (3.0/(x[n-1]-x[n-2]))*(dfR - (f[n-1] - f[n-2])/(x[n-1]-x[n-2]));
149 
150  d2f[n-1] = (un - qn*u[n-2])/(qn*d2f[n-2] + 1.0);
151  for (i = n-2; i >= 0; i--) d2f[i] = d2f[i]*d2f[i+1] + u[i];
152 
153 /* -----------------------------------------------------
154  Construct spline coefficients, Eq. [3.3.3], [3.3.4]
155  ----------------------------------------------------- */
156 
157  for (i = 0; i < n-1; i++){
158  dx = x[i+1] - x[i];
159  a[i] = dx*dx/6.0*(d2f[i+1] - d2f[i]);
160  b[i] = dx*dx*0.5*d2f[i];
161  c[i] = f[i+1] - f[i] - dx*dx/6.0*(2.0*d2f[i] + d2f[i+1]);
162  d[i] = f[i];
163 
164  /* -- check monotonicity (quit if not satisfied) -- */
165 
166  aa = 3.0*a[i];
167  bb = 2.0*b[i];
168  cc = 1.0*c[i];
169  det = bb*bb - 4.0*aa*cc;
170  if (det >= 0.0){
171  xm = (-bb + sqrt(det))/(2.0*aa);
172  xp = (-bb - sqrt(det))/(2.0*aa);
173 
174  if ((xm > 0.0 && xm < 1.0) || (xp > 0.0 && xp < 1.0)){
175  printf ("! SplineCoeffs(): cubic is not monotonic, in [%8.3e, %8.3e]\n",
176  x[i], x[i+1]);
177  printf ("! xm = %8.3e, xp = %8.3e\n",xm,xp);
178  QUIT_PLUTO(1);
179  }
180  }
181 
182  /* ---------------------------------------------------------
183  Reset a = b = 0 if the interpolant does not have
184  enough curvature.
185  This is done by checking the size of |y''/y| at x = 0.5
186  (midpoint).
187  --------------------------------------------------------- */
188 
189  xm = 6.0*a[i]*0.5 + 2.0*b[i];
190  xp = d[i] + 0.5*(c[i] + 0.5*(b[i] + 0.5*a[i]));
191  det = fabs(xm/xp);
192  if (det < 1.e-16){
193  nfail++;
194  a[i] = b[i] = 0.0;
195  }
196  }
197 }
a
static double a
Definition: init.c:135
n
static int n
Definition: analysis.c:3
menu.c
tuple c
Definition: menu.py:375
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
i
int i
Definition: analysis.c:2
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the caller graph for this function: