1 /*
2 * - - - - - - - - - - - - -
3 * g a l _ r k f c k s 4 5
4 * - - - - - - - - - - - - -
5 *
6 * This routine is part of the General Astrodynamics Library
7 *
8 * Description:
9 *
10 * This routine takes a Runge-Kutte-Fehlberg-Cash-Karp 4(5) step
11 *
12 * Status:
13 *
14 * Math special function.
15 *
16 * Given:
17 *
18 * y d[] dependent variable vector
19 * dydx d[] derivative of dependent variable vector
20 * n i Number of equations to integrate
21 * x d Independent variable value
22 * h d Stepsize
23 * *derivs func User defined function for calculating the right hand side derivatives
24 * *derivsp i Pointer to parameters structure for derivs routine
25 *
26 * Returned:
27 *
28 * yout d[] Ending y values
29 * yerr d[] Errors
30 * gal_rkfs45 i Status
31 * 0 = Success
32 * 1 = Failed to allocate memory
33 *
34 * Notes:
35 *
36 * 1) The parameters (Cash-Karp version) from Numerical Recipes for RKF45
37 * These values are taken from the c code and not from the table on page 717
38 * which has different values (which don't work, but look like they almost do). The
39 * Cash-Karp values seem to make the routine a bit faster compared to the Fehlberg values.
40 *
41 * References:
42 *
43 * Numerical Receipes in C
44 * The Art of Scientific Computing
45 * Second Edition
46 * by William H. Press, Saul A. Teukolsky, William T. Vettering & Brian P. Flannery
47 * Pages 710 - 722
48 *
49 * This revision:
50 *
51 * 2008 May 18
52 *
53 * Copyright (C) 2008 Paul C. L. Willmott. See notes at end.
54 *
55 *-----------------------------------------------------------------------
56 */
57
58 #include <stdlib.h>
59 #include "gal_const.h"
60 #include "gal_rkfcks45.h"
61
gal_rkfcks45(double y[],double dydx[],int n,double x,double h,double yout[],double yerr[],void (* derivs)(double,double[],double[],int *),int * derivsp)62 int gal_rkfcks45
63 (
64 double y[],
65 double dydx[],
66 int n,
67 double x,
68 double h,
69 double yout[],
70 double yerr[],
71 void ( *derivs ) ( double, double [], double [], int * ),
72 int *derivsp
73 )
74
75 {
76
77 struct {
78 double a[6] ;
79 double b[6][5] ;
80 double c[6] ;
81 double chat[6] ;
82 } param = {
83
84 {
85 0.0,
86 0.2,
87 0.3,
88 0.6,
89 1.0,
90 0.875,
91 } ,
92
93 {
94 { 0.0, 0.0, 0.0, 0.0, 0.0, } ,
95 { 0.2, 0.0, 0.0, 0.0, 0.0, } ,
96 { 3.0/40.0, 9.0/40.0, 0.0, 0.0, 0.0, } ,
97 { 0.3, -0.9, 1.2, 0.0, 0.0, } ,
98 { -11.0/54.0, 2.5, -70.0/27.0, 35.0/27.0, 0.0, } ,
99 { 1631.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0, } ,
100 } ,
101
102 {
103 37.0 / 378.0,
104 0.0,
105 250.0 / 621.0,
106 125.0 / 594.0,
107 0.0,
108 512.0 / 1771.0,
109 } ,
110
111 {
112 2825.0 / 27648.0,
113 0.0,
114 18575.0 / 48384.0,
115 13525.0 / 55296.0,
116 277.0 / 14336.0,
117 1.0 / 4.0,
118 } ,
119
120 } ;
121
122 #define A(i) param.a[i]
123 #define B(i,j) param.b[i][j]
124 #define C(i) param.c[i]
125 #define CHAT param.chat[i]
126 #define DC(i) (param.c[i] - param.chat[i])
127
128 int i ;
129 double *ak2, *ak3, *ak4, *ak5, *ak6, *ytemp ;
130
131 /*
132 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133 */
134
135 /*
136 * Allocate Workspace
137 */
138
139 ak2 = malloc ( n * sizeof ( double ) ) ;
140 if ( ak2 == NULL ) {
141 return 1 ;
142 }
143 ak3 = malloc ( n * sizeof ( double ) ) ;
144 if ( ak3 == NULL ) {
145 free ( ak2 ) ;
146 return 1 ;
147 }
148 ak4 = malloc ( n * sizeof ( double ) ) ;
149 if ( ak4 == NULL ) {
150 free ( ak3 ) ;
151 free ( ak2 ) ;
152 return 1 ;
153 }
154 ak5 = malloc ( n * sizeof ( double ) ) ;
155 if ( ak5 == NULL ) {
156 free ( ak4 ) ;
157 free ( ak3 ) ;
158 free ( ak2 ) ;
159 return 1 ;
160 }
161 ak6 = malloc ( n * sizeof ( double ) ) ;
162 if ( ak6 == NULL ) {
163 free ( ak5 ) ;
164 free ( ak4 ) ;
165 free ( ak3 ) ;
166 free ( ak2 ) ;
167 return 1 ;
168 }
169
170 ytemp = malloc ( n * sizeof ( double ) ) ;
171 if ( ytemp == NULL ) {
172 free ( ak6 ) ;
173 free ( ak5 ) ;
174 free ( ak4 ) ;
175 free ( ak3 ) ;
176 free ( ak2 ) ;
177 return 1 ;
178 }
179
180 /*
181 * Do it!
182 */
183
184 for ( i = 0; i < n; i++ ) {
185 ytemp[i] = y[i] + B(1,0) * h * dydx[i] ;
186 }
187 ( *derivs ) ( x + A(1) * h, ytemp, ak2, derivsp ) ;
188
189 for ( i = 0; i < n; i++ ) {
190 ytemp[i] = y[i] + h * (B(2,0) * dydx[i] + B(2,1) * ak2[i]) ;
191 }
192 ( *derivs ) ( x + A(2) * h, ytemp, ak3, derivsp ) ;
193
194 for ( i = 0; i < n; i++ ) {
195 ytemp[i] = y[i] + h * (B(3,0) * dydx[i] + B(3,1) * ak2[i] + B(3,2) * ak3[i]) ;
196 }
197 ( *derivs ) ( x + A(3) * h, ytemp, ak4, derivsp ) ;
198
199 for ( i = 0; i < n; i++ ) {
200 ytemp[i] = y[i] + h * (B(4,0) * dydx[i] + B(4,1) * ak2[i] + B(4,2) * ak3[i] + B(4,3) * ak4[i]) ;
201 }
202 ( *derivs ) ( x + A(4) * h, ytemp, ak5, derivsp ) ;
203
204 for ( i = 0; i < n; i++ ) {
205 ytemp[i] = y[i] + h * (B(5,0) * dydx[i] + B(5,1) * ak2[i] + B(5,2) * ak3[i] + B(5,3) * ak4[i] + B(5,4) * ak5[i]) ;
206 }
207 ( *derivs )(x + A(5) * h, ytemp, ak6, derivsp ) ;
208
209 for ( i = 0; i < n; i++ ) {
210 yout[i] = y[i] + h * (C(0) * dydx[i] + C(1) * ak2[i] + C(2) * ak3[i] + C(3) * ak4[i] + C(4) * ak5[i] + C(5) * ak6[i]) ;
211 }
212
213 for ( i = 0; i < n; i++ ) {
214 yerr[i] = h * (DC(0) * dydx[i] + DC(1) * ak2[i] + DC(2) * ak3[i] + DC(3) * ak4[i] + DC(4) * ak5[i] + DC(5) * ak6[i]) ;
215 }
216
217 /*
218 * De-allocate Workspace Memory
219 */
220
221 free ( ytemp ) ;
222 free ( ak6 ) ;
223 free ( ak5 ) ;
224 free ( ak4 ) ;
225 free ( ak3 ) ;
226 free ( ak2 ) ;
227
228 return GAL_SUCCESS ;
229
230 /*
231 * Finished.
232 */
233
234 }
235
236 /*
237 * gal - General Astrodynamics Library
238 * Copyright (C) 2008 Paul C. L. Willmott
239 *
240 * This program is free software; you can redistribute it and/or modify
241 * it under the terms of the GNU General Public License as published by
242 * the Free Software Foundation; either version 2 of the License, or
243 * (at your option) any later version.
244 *
245 * This program is distributed in the hope that it will be useful,
246 * but WITHOUT ANY WARRANTY; without even the implied warranty of
247 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
248 * GNU General Public License for more details.
249 *
250 * You should have received a copy of the GNU General Public License along
251 * with this program; if not, write to the Free Software Foundation, Inc.,
252 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
253 *
254 * Contact:
255 *
256 * Paul Willmott
257 * vp9mu@amsat.org
258 */
259
260