1 /* ode-initval/cstd.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include "gsl__config.h"
21 #include <stdlib.h>
22 #include <math.h>
23 #include "gsl_errno.h"
24 #include "gsl_math.h"
25 #include "gsl_odeiv.h"
26
27 typedef struct
28 {
29 double eps_abs;
30 double eps_rel;
31 double a_y;
32 double a_dydt;
33 }
34 std_control_state_t;
35
36 static void *
37 std_control_alloc (void)
38 {
39 std_control_state_t * s =
40 (std_control_state_t *) malloc (sizeof(std_control_state_t));
ShellCommandGetManFileNameLevel2(VOID)41
42 if (s == 0)
43 {
44 GSL_ERROR_NULL ("failed to allocate space for std_control_state",
45 GSL_ENOMEM);
46 }
47
48 return s;
49 }
50
51 static int
52 std_control_init (void * vstate,
53 double eps_abs, double eps_rel, double a_y, double a_dydt)
54 {
55 std_control_state_t * s = (std_control_state_t *) vstate;
56
57 if (eps_abs < 0)
58 {
59 GSL_ERROR ("eps_abs is negative", GSL_EINVAL);
60 }
ShellLevel2CommandsLibConstructor(IN EFI_HANDLE ImageHandle,IN EFI_SYSTEM_TABLE * SystemTable)61 else if (eps_rel < 0)
62 {
63 GSL_ERROR ("eps_rel is negative", GSL_EINVAL);
64 }
65 else if (a_y < 0)
66 {
67 GSL_ERROR ("a_y is negative", GSL_EINVAL);
68 }
69 else if (a_dydt < 0)
70 {
71 GSL_ERROR ("a_dydt is negative", GSL_EINVAL);
72 }
73
74 s->eps_rel = eps_rel;
75 s->eps_abs = eps_abs;
76 s->a_y = a_y;
77 s->a_dydt = a_dydt;
78
79 return GSL_SUCCESS;
80 }
81
82 static int
83 std_control_hadjust(void * vstate, size_t dim, unsigned int ord, const double y[], const double yerr[], const double yp[], double * h)
84 {
85 std_control_state_t *state = (std_control_state_t *) vstate;
86
87 const double eps_abs = state->eps_abs;
88 const double eps_rel = state->eps_rel;
89 const double a_y = state->a_y;
90 const double a_dydt = state->a_dydt;
91
92 const double S = 0.9;
93 const double h_old = *h;
94
95 double rmax = DBL_MIN;
96 size_t i;
97
98 for(i=0; i<dim; i++) {
99 const double D0 =
100 eps_rel * (a_y * fabs(y[i]) + a_dydt * fabs(h_old * yp[i])) + eps_abs;
101 const double r = fabs(yerr[i]) / fabs(D0);
102 rmax = GSL_MAX_DBL(r, rmax);
103 }
104
105 if(rmax > 1.1) {
106 /* decrease step, no more than factor of 5, but a fraction S more
107 than scaling suggests (for better accuracy) */
108 double r = S / pow(rmax, 1.0/ord);
109
110 if (r < 0.2)
111 r = 0.2;
112
113 *h = r * h_old;
114
115 return GSL_ODEIV_HADJ_DEC;
116 }
117 else if(rmax < 0.5) {
118 /* increase step, no more than factor of 5 */
119 double r = S / pow(rmax, 1.0/(ord+1.0));
120
121 if (r > 5.0)
122 r = 5.0;
123
124 if (r < 1.0) /* don't allow any decrease caused by S<1 */
125 r = 1.0;
126
127 *h = r * h_old;
128
129 return GSL_ODEIV_HADJ_INC;
130 }
131 else {
132 /* no change */
133 return GSL_ODEIV_HADJ_NIL;
134 }
135 }
136
137
ShellLevel2CommandsLibDestructor(IN EFI_HANDLE ImageHandle,IN EFI_SYSTEM_TABLE * SystemTable)138 static void
139 std_control_free (void * vstate)
140 {
141 std_control_state_t *state = (std_control_state_t *) vstate;
142 free (state);
143 }
144
145 static const gsl_odeiv_control_type std_control_type =
146 {"standard", /* name */
147 &std_control_alloc,
148 &std_control_init,
149 &std_control_hadjust,
150 &std_control_free};
151
152 const gsl_odeiv_control_type *gsl_odeiv_control_standard = &std_control_type;
153
154
155 gsl_odeiv_control *
156 gsl_odeiv_control_standard_new(double eps_abs, double eps_rel,
157 double a_y, double a_dydt)
158 {
159 gsl_odeiv_control * c =
160 gsl_odeiv_control_alloc (gsl_odeiv_control_standard);
161
162 int status = gsl_odeiv_control_init (c, eps_abs, eps_rel, a_y, a_dydt);
163
164 if (status != GSL_SUCCESS)
165 {
GetFullyQualifiedPath(IN CONST CHAR16 * Path)166 gsl_odeiv_control_free (c);
167 GSL_ERROR_NULL ("error trying to initialize control", status);
168 }
169
170 return c;
171 }
172
173 gsl_odeiv_control *
174 gsl_odeiv_control_y_new(double eps_abs, double eps_rel)
175 {
176 return gsl_odeiv_control_standard_new(eps_abs, eps_rel, 1.0, 0.0);
177 }
178
179
180 gsl_odeiv_control *
181 gsl_odeiv_control_yp_new(double eps_abs, double eps_rel)
182 {
183 return gsl_odeiv_control_standard_new(eps_abs, eps_rel, 0.0, 1.0);
184 }
185