| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 /// \file ODESolver.cc
28 /// \brief Implementation of the ODESolver cla
29
30 #include "ODESolver.hh"
31 #include <iostream>
32 #include <cmath>
33
34 //....oooOO0OOooo........oooOO0OOooo........oo
35
36 std::vector<double> operator*(const std::vecto
37 {
38 std::vector<double> vout;
39 for (auto const val : v) vout.push_back(val*
40 return vout;
41 }
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44
45 std::vector<double> operator+(const std::vecto
46 {
47 std::vector<double> vout;
48 for (auto const val : v) vout.push_back(val
49 return vout;
50 }
51
52 //....oooOO0OOooo........oooOO0OOooo........oo
53
54 std::vector<double> operator+(const std::vecto
55 {
56 std::vector<double> vout;
57 for (size_t i=0;i<v1.size();i++) vout.push_b
58 return vout;
59 }
60
61 //....oooOO0OOooo........oooOO0OOooo........oo
62
63 ODESolver::ODESolver(): fNstepsForObserver(1)
64 {}
65
66 //....oooOO0OOooo........oooOO0OOooo........oo
67
68 double ODESolver::RungeKutta_Fehlberg( std::fu
69 func,std::vector<double> &y, double t, double
70 {
71 //based on https://en.wikipedia.org/wiki/Run
72 const int nk=6;
73 double h = stepsize;
74 double CH[nk]={47./450.,0,12./25.,32./255.,1
75 double CT[nk]={-1./150.,0.,3./100.,-16./75.,
76 double A[nk]={0.,2./9.,1./3.,3./4.,1.,5./6.}
77 double B21=2./9., B31=1./12., B41=69./128.,
78 double B32=1./4., B42=-243./128., B52=27./5.
79 double B43=135./64., B53=-27./5., B63=13./16
80 double B54=16./15., B64=4./27.;
81 double B65=5./144.;
82 double maxError = 1.;
83 std::vector<double> k1 = func(t+A[0]*h,y)*h;
84 std::vector<double> k2 = func(t+A[1]*h,y + k
85 std::vector<double> k3 = func(t+A[2]*h,y + k
86 std::vector<double> k4 = func(t+A[3]*h,y + k
87 std::vector<double> k5 = func(t+A[4]*h,y + k
88 std::vector<double> k6 = func(t+A[5]*h,y + k
89 y = y + k1*CH[0] + k2*CH[1] + k3*CH[2] + k4*
90 auto TE = k1*CT[0] + k2*CT[1] + k3*CT[2] + k
91 absValuesVector(TE);
92 maxError = *std::max_element(TE.begin(),TE.e
93
94 k1.clear(); k1.shrink_to_fit();
95 k2.clear(); k2.shrink_to_fit();
96 k3.clear(); k3.shrink_to_fit();
97 k4.clear(); k4.shrink_to_fit();
98 k5.clear(); k5.shrink_to_fit();
99 k6.clear(); k6.shrink_to_fit();
100 TE.clear(); TE.shrink_to_fit();
101 return maxError;
102 }
103
104 //....oooOO0OOooo........oooOO0OOooo........oo
105
106 void ODESolver::Embedded_RungeKutta_Fehlberg(
107 std::function<std::vector<double>(double,std
108 double start,double end,double stepsize,doub
109 std::vector<double> *time_observer,std::vect
110 {
111 double t = start;
112 double h = stepsize;
113 int nsteps = 0;
114 if (h < 0) h = (end - start)/(10000.);
115 if (time_observer) time_observer->push_back(
116 if (state_observer) state_observer->push_bac
117 auto ytemp = y;
118 while (t < end)
119 {
120 ytemp = y;
121 double maxerror = RungeKutta_Fehlberg(func
122 double scale = 0.9*std::pow(epsilon/maxerr
123
124 double hnew = h*scale;
125 while (maxerror > epsilon)
126 {
127 ytemp = y;
128 maxerror = RungeKutta_Fehlberg(func,ytem
129 scale = 0.9*std::pow(epsilon/maxerror,1.
130 hnew = hnew*scale;
131 }
132 h = hnew;
133 y = ytemp;
134 t += h;
135 if (t > end) break;
136 if ( time_observer || state_observer) {
137 nsteps++;
138 if (nsteps%fNstepsForObserver == 0) {
139 if (time_observer) time_observer->push
140 if (state_observer) state_observer->pu
141 nsteps = 0;
142 }
143 }
144 }
145
146 ytemp.clear(); ytemp.shrink_to_fit();
147 }
148
149 //....oooOO0OOooo........oooOO0OOooo........oo
150
151 void ODESolver::RungeKutta4(
152 std::function<std::vector<double>(double,std
153 double start,double end,double stepsize,
154 std::vector<double> *time_observer,std::ve
155 {
156 double t = start;
157 double h = stepsize;
158 int nsteps = 0;
159 if (h < 0) h = (end - start)/(10000.);
160 if (time_observer) time_observer->push_back(
161 if (state_observer) state_observer->push_bac
162 while (t < end)
163 {
164 std::vector<double> k1 = func(t,y)*h;
165 std::vector<double> k2 = func(t+0.5*h,y +
166 std::vector<double> k3 = func(t+0.5*h,y +
167 std::vector<double> k4 = func(t+h,y + k3)*
168 t += h;
169 if (t > end) break;
170 y = y +(k1 +k2*2+k3*2+k4)*(1./6.0);
171 if ( time_observer || state_observer) {
172 nsteps++;
173 if (nsteps%fNstepsForObserver == 0) {
174 if (time_observer) time_observer->push
175 if (state_observer) state_observer->pu
176 nsteps = 0;
177 }
178 }
179 k1.clear(); k1.shrink_to_fit();
180 k2.clear(); k2.shrink_to_fit();
181 k3.clear(); k3.shrink_to_fit();
182 k4.clear(); k4.shrink_to_fit();
183 }
184 }
185
186 //....oooOO0OOooo........oooOO0OOooo........oo