| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file ODESolver.cc
28 /// \brief Implementation of the ODESolver class
29
30 #include "ODESolver.hh"
31 #include <iostream>
32 #include <cmath>
33
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
35
36 std::vector<double> operator*(const std::vector<double> v, double alfa)
37 {
38 std::vector<double> vout;
39 for (auto const val : v) vout.push_back(val*alfa);
40 return vout;
41 }
42
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
44
45 std::vector<double> operator+(const std::vector<double> v, double alfa)
46 {
47 std::vector<double> vout;
48 for (auto const val : v) vout.push_back(val + alfa);
49 return vout;
50 }
51
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
53
54 std::vector<double> operator+(const std::vector<double> v1, const std::vector<double> v2)
55 {
56 std::vector<double> vout;
57 for (size_t i=0;i<v1.size();i++) vout.push_back(v1.at(i) + v2.at(i));
58 return vout;
59 }
60
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
62
63 ODESolver::ODESolver(): fNstepsForObserver(1)
64 {}
65
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
67
68 double ODESolver::RungeKutta_Fehlberg( std::function<std::vector<double>(double,std::vector<double>)>
69 func,std::vector<double> &y, double t, double stepsize)
70 {
71 //based on https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method
72 const int nk=6;
73 double h = stepsize;
74 double CH[nk]={47./450.,0,12./25.,32./255.,1./30.,6./25.};
75 double CT[nk]={-1./150.,0.,3./100.,-16./75.,-1./20.,6./25.};
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., B51=-17./12., B61=65./432.;
78 double B32=1./4., B42=-243./128., B52=27./5., B62=13./16.;
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 + k1*B21)*h;
85 std::vector<double> k3 = func(t+A[2]*h,y + k1*B31 + k2*B32)*h;
86 std::vector<double> k4 = func(t+A[3]*h,y + k1*B41 + k2*B42 + k3*B43)*h;
87 std::vector<double> k5 = func(t+A[4]*h,y + k1*B51 + k2*B52 + k3*B53 + k4*B54)*h;
88 std::vector<double> k6 = func(t+A[5]*h,y + k1*B61 + k2*B62 + k3*B63 + k4*B64 + k5*B65)*h;
89 y = y + k1*CH[0] + k2*CH[1] + k3*CH[2] + k4*CH[3] + k5*CH[4] + k6*CH[5];
90 auto TE = k1*CT[0] + k2*CT[1] + k3*CT[2] + k4*CT[3] + k5*CT[4] + k6*CT[5];
91 absValuesVector(TE);
92 maxError = *std::max_element(TE.begin(),TE.end());
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........oooOO0OOooo........oooOO0OOooo.....
105
106 void ODESolver::Embedded_RungeKutta_Fehlberg(
107 std::function<std::vector<double>(double,std::vector<double>)> func, std::vector<double> &y,
108 double start,double end,double stepsize,double epsilon,
109 std::vector<double> *time_observer,std::vector<std::vector<double>> *state_observer)
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(t);
116 if (state_observer) state_observer->push_back(y);
117 auto ytemp = y;
118 while (t < end)
119 {
120 ytemp = y;
121 double maxerror = RungeKutta_Fehlberg(func,ytemp,t,h);
122 double scale = 0.9*std::pow(epsilon/maxerror,1./5.);
123
124 double hnew = h*scale;
125 while (maxerror > epsilon)
126 {
127 ytemp = y;
128 maxerror = RungeKutta_Fehlberg(func,ytemp,t,hnew);
129 scale = 0.9*std::pow(epsilon/maxerror,1./5.);
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_back(t);
140 if (state_observer) state_observer->push_back(y);
141 nsteps = 0;
142 }
143 }
144 }
145
146 ytemp.clear(); ytemp.shrink_to_fit();
147 }
148
149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
150
151 void ODESolver::RungeKutta4(
152 std::function<std::vector<double>(double,std::vector<double>)> func, std::vector<double> &y,
153 double start,double end,double stepsize,
154 std::vector<double> *time_observer,std::vector<std::vector<double>> *state_observer)
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(t);
161 if (state_observer) state_observer->push_back(y);
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 + k1*0.5)*h;
166 std::vector<double> k3 = func(t+0.5*h,y + k2*0.5)*h;
167 std::vector<double> k4 = func(t+h,y + k3)*h;
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_back(t);
175 if (state_observer) state_observer->push_back(y);
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........oooOO0OOooo........oooOO0OOooo.....