git@sv / gsl.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
history | raw | HEAD
preparing 2.8 release
[gsl.git] / contrib / euler.c
1 /* euler.c
2  * 
3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
4  * Copyright (C) 2004 Alfonso Acosta & Daniel Rodríguez
5  * 
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 3 of the License, or (at
9  * your option) any later version.
10  * 
11  * This program is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14  * General Public License for more details.
15  * 
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19  */
20
21 /* Euler Method*/
22
23 /* Authors:  Alfonso Acosta & Daniel Rodríguez
24  */
25
26 #include <gsl/gsl_errno.h>
27 #include <gsl/gsl_odeiv.h>
28
29 #include "odeiv_util.h"
30
31 typedef struct
32 {
33   double *k;
34 }
35 euler_state_t;
36
37 static void *
38 euler_alloc (size_t dim)
39 {
40   euler_state_t *state = (euler_state_t *) malloc (sizeof (euler_state_t));
41
42   if (state == 0)
43     {
44       GSL_ERROR_NULL ("failed to allocate space for euler_state", GSL_ENOMEM);
45     }
46
47   state->k = (double *) malloc (dim * sizeof (double));
48
49   if (state->k == 0)
50     {
51       free (state);
52       GSL_ERROR_NULL ("failed to allocate space for k", GSL_ENOMEM);
53     }
54
55   return state;
56 }
57
58
59 static int
60 euler_apply (void *vstate,
61            size_t dim,
62            double t,
63            double h,
64            double y[],
65            double yerr[],
66            const double dydt_in[],
67            double dydt_out[], 
68            const gsl_odeiv_system * sys)
69 {
70   euler_state_t *state = (euler_state_t *) vstate;
71   
72   double temp;
73   size_t i;
74   int status = 0;
75
76   double *const k = state->k;
77
78
79   if (dydt_in != NULL)
80     {
81       DBL_MEMCPY (k, dydt_in, dim);
82     }
83   else
84     {
85       int s = GSL_ODEIV_FN_EVAL (sys, t, y, k);
86       GSL_STATUS_UPDATE (&status, s);
87     }
88   
89
90   for (i = 0; i < dim; i++)
91     {
92       temp = y[i]; /* save y[i] */ 
93       y[i] = h * k[i];   
94       yerr[i] = h * y[i];
95       y[i] += temp;
96       if (dydt_out != NULL)
97           dydt_out[i] = k[i];
98     }
99
100   return status;
101 }
102
103 static int
104 euler_reset (void *vstate, size_t dim)
105 {
106   euler_state_t *state = (euler_state_t *) vstate;
107
108   DBL_ZERO_MEMSET (state->k, dim);
109
110   return GSL_SUCCESS;
111 }
112
113 static unsigned int
114 euler_order (void *vstate)
115 {
116   euler_state_t *state = (euler_state_t *) vstate;
117   state = 0; /* prevent warnings about unused parameters */
118   return 1;
119 }
120
121 static void
122 euler_free (void *vstate)
123 {
124   euler_state_t *state = (euler_state_t *) vstate;
125   free (state->k);
126   free (state);
127 }
128
129 static const gsl_odeiv_step_type euler_type = { "euler",    /* name */
130   1,                            /* can use dydt_in */
131   0,                            /* gives exact dydt_out */
132   &euler_alloc,
133   &euler_apply,
134   &euler_reset,
135   &euler_order,
136   &euler_free
137 };
138
139 const gsl_odeiv_step_type *gsl_odeiv_step_euler = &euler_type;
GNU Scientific Library