Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInspiral 5.0.3.1-6c6b863
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALRungeKutta4.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 Duncan Brown, Gareth Jones, Jolien Creighton, B.S. Sathyaprakash, Craig Robinson
3*
4* This program is free software; you can redistribute it and/or modify
5* it under the terms of the GNU General Public License as published by
6* the Free Software Foundation; either version 2 of the License, or
7* (at your option) any later version.
8*
9* This program is distributed in the hope that it will be useful,
10* but WITHOUT ANY WARRANTY; without even the implied warranty of
11* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12* GNU General Public License for more details.
13*
14* You should have received a copy of the GNU General Public License
15* along with with program; see the file COPYING. If not, write to the
16* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17* MA 02110-1301 USA
18*/
19
20/**
21 * \author Robinson, C. A.
22 * \file
23 * \ingroup LALInspiral_h
24 *
25 * \brief NONE
26 *
27 * ### Prototypes ###
28 *
29 * <tt>LALRungeKutta4()</tt>
30 * <ul>
31 * <li> \c n: The number of coupled equations being integrated.</li>
32 * <li> \c yout: The output values for the system after the time-step.</li>
33 * <li> \c input: The input for the system</li>
34 * <li> \c integrator Required for the GSL integratior. Created using XLALRungeKutta4Init().</li>
35 * <li> \c params Parameters to be passed to the derivative function</li>
36 * </ul>
37 *
38 * ### Description ###
39 *
40 * The code \ref LALRungeKutta4.c solves a system of \f$n\f$ coupled first--order differential equations.
41 * Internally, it uses the gsl routines for performing adaptive step evolution of the system, but to the outside
42 * user, it returns results for a fixed step size.
43 *
44 * Prior to evolving a system using LALRungeKutta4(), it is necessary to create the GSL integrator using
45 * XLALRungeKutta4Init(). Once the evolution of the system has finished, this integrator should then
46 * be freed using XLALRungeKutta4Free().
47 *
48 * ### Algorithm ###
49 *
50 *
51 * ### Uses ###
52 *
53 * None.
54 *
55 * ### Notes ###
56 *
57 */
58
59#include <lal/XLALGSL.h>
60#include <lal/LALInspiral.h>
61
62#ifdef __GNUC__
63#define UNUSED __attribute__ ((unused))
64#else
65#define UNUSED
66#endif
67
70 void *params;
71};
72
74 REAL8 t,
75 const REAL8 y[],
76 REAL8 dydx[],
77 void *params);
78
79/* Function for allocating memory and setting up the GSL integrator */
80
82 rk4In *input
83 )
84{
85
86 rk4GSLIntegrator *integrator = NULL;
87
88 /* Check we have an input */
89 if (!input)
91
92 /* Allocate memory for the integrator structure */
93 if (!(integrator = (rk4GSLIntegrator *) LALCalloc(1, sizeof(rk4GSLIntegrator))))
94 {
96 }
97
98 integrator->input = input;
99
100 /* Set the algorithm to 4th-order Runge-Kutta */
101 integrator->type = gsl_odeiv_step_rkf45;
102
103 /* Allocate memory for data values */
104 if (!(integrator->y = (REAL8 *) LALMalloc(n * sizeof(REAL8))))
105 {
106 LALFree(integrator);
108 }
109
110 /* Initialise GSL integrator */
111 XLAL_CALLGSL( integrator->step = gsl_odeiv_step_alloc(integrator->type, n) );
112 XLAL_CALLGSL( integrator->control = gsl_odeiv_control_standard_new(1.0e-2, 1.0e-2, 1.0, 1.0) );
113 XLAL_CALLGSL( integrator->evolve = gsl_odeiv_evolve_alloc(n) );
114
115 /* Check the integrator is allocated correctly */
116 if (!(integrator->step) || !(integrator->control) || !(integrator->evolve))
117 {
118 XLALRungeKutta4Free( integrator );
120 }
121
122 return integrator;
123}
124
125
126void
129 REAL8Vector *yout,
130 rk4GSLIntegrator *integrator,
131 void *params
132 )
133{
134
136
137 XLAL_PRINT_DEPRECATION_WARNING("XLALRungeKutta4");
138
139 if ( XLALRungeKutta4( yout, integrator, params ) == XLAL_FAILURE )
140 ABORTXLAL( status );
141
142 RETURN( status );
143}
144
145
146int
148 REAL8Vector *yout,
149 rk4GSLIntegrator *integrator,
150 void *params
151 )
152{
153
154 int gslStatus;
155
156 INT4 i;
157 REAL8 t = 0.0;
158 struct RungeGSLParams gslParams;
159 rk4In *input = NULL;
160 REAL8 h;
161 gsl_odeiv_system sys;
162
163 if ( !yout )
165
166 if ( !yout->data )
168
169 if ( !integrator )
171
172 if ( !integrator->input )
174
175 if ( !params )
177
178 /* Initialise GSL integrator */
179
180 input = integrator->input;
181 h = input->h;
182
183 gslParams.input = input;
184 gslParams.params = params;
185
186 sys.function = derivativeGSLWrapper;
187 sys.jacobian = NULL;
188 sys.dimension = input->n;
189 sys.params = &gslParams;
190
191
192 memcpy( integrator->y, input->y->data, input->n * sizeof(REAL8));
193
194 /* Evolve the system */
195 while (t < input->h)
196 {
197 REAL8 tOld = t;
198 XLAL_CALLGSL( gslStatus = gsl_odeiv_evolve_apply(integrator->evolve,
199 integrator->control, integrator->step, &sys,
200 &t, input->h, &h, integrator->y) );
201
202 /*printf("h = %e, t = %e\n", h, t);*/
203 if ( gslStatus != GSL_SUCCESS )
204 {
205 XLALPrintError( "Failure in gsl_odeiv_evolve_apply\n" );
207 }
208
209 /* In case integration becomes degenerate */
210 if (t == tOld)
211 {
212 for (i=0; i<input->n; i++)
213 yout->data[i] = 0.0;
214
215 XLALPrintError( "Time step grown too small!\n" );
217 }
218 }
219
220 memcpy( yout->data, integrator->y, input->n * sizeof(REAL8));
221
222 return XLAL_SUCCESS;
223}
224
225
226/* Function for freeing up memory for the GSL integrator */
227
229{
230
231
232 if (!integrator) XLAL_ERROR_VOID(XLAL_EFAULT);
233
234 /* Free the GSL integrator controls etc */
235 if (integrator->evolve) XLAL_CALLGSL( gsl_odeiv_evolve_free(integrator->evolve) );
236 if (integrator->control) XLAL_CALLGSL( gsl_odeiv_control_free(integrator->control) );
237 if (integrator->step) XLAL_CALLGSL( gsl_odeiv_step_free(integrator->step) );
238 LALFree( integrator->y );
239 LALFree( integrator );
240 return;
241}
242
243
244/* A simple wrapper function to allow GSL to use the LAL
245 derivative functions */
247 REAL8 UNUSED t,
248 const REAL8 y[],
249 REAL8 dydx[],
250 void *params)
251{
252 struct RungeGSLParams *in = (struct RungeGSLParams *)params;
253 REAL8Vector dyVect;
254 REAL8Vector *yVect = in->input->yt;
255
256 memcpy(yVect->data, y, in->input->n * sizeof(REAL8));
257
258 dyVect.length = in->input->n;
259 dyVect.data = dydx;
260 in->input->function(yVect, &dyVect, in->params);
261 return GSL_SUCCESS;
262}
#define LALCalloc(m, n)
#define LALMalloc(n)
#define LALFree(p)
void XLALRungeKutta4Free(rk4GSLIntegrator *integrator)
static int derivativeGSLWrapper(REAL8 t, const REAL8 y[], REAL8 dydx[], void *params)
void LALRungeKutta4(LALStatus *status, REAL8Vector *yout, rk4GSLIntegrator *integrator, void *params)
rk4GSLIntegrator * XLALRungeKutta4Init(INT4 n, rk4In *input)
int XLALRungeKutta4(REAL8Vector *yout, rk4GSLIntegrator *integrator, void *params)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define ABORTXLAL(sp)
#define XLAL_CALLGSL(statement)
double i
double REAL8
int32_t INT4
#define XLAL_ERROR_VOID(...)
#define XLAL_ERROR_NULL(...)
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_PRINT_DEPRECATION_WARNING(replacement)
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFUNC
XLAL_EFAILED
XLAL_FAILURE
list y
int n
REAL8 * data
Structure containing steps and controls for the GSL Runge-Kutta solver.
Definition: LALInspiral.h:637
gsl_odeiv_control * control
Definition: LALInspiral.h:640
gsl_odeiv_step * step
Definition: LALInspiral.h:639
const gsl_odeiv_step_type * type
Definition: LALInspiral.h:638
gsl_odeiv_evolve * evolve
Definition: LALInspiral.h:641
Structure used as an input to Runge-Kutta solver.
Definition: LALInspiral.h:620
REAL8Vector * yt
Definition: LALInspiral.h:625
REAL8 h
Definition: LALInspiral.h:628
INT4 n
Definition: LALInspiral.h:629
REAL8Vector * y
Definition: LALInspiral.h:623
TestFunction * function
Definition: LALInspiral.h:621