analysis/src/ODESolver.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// \file ODESolver.cc
0028 /// \brief Implementation of the ODESolver class
0029
0030 #include "ODESolver.hh"
0031 #include <iostream>
0032 #include <cmath>
0033
0034 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0035
0036 std::vector<double> operator*(const std::vector<double> v, double alfa)
0037 {
0038     std::vector<double> vout;
0039     for (auto const val : v) vout.push_back(val*alfa);
0040     return vout;
0041 }
0042
0043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0044
0045 std::vector<double> operator+(const std::vector<double> v, double alfa)
0046 {
0047     std::vector<double> vout;
0048     for (auto const val : v) vout.push_back(val + alfa);
0049     return vout;
0050 }
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0053
0054 std::vector<double> operator+(const std::vector<double> v1, const std::vector<double> v2)
0055 {
0056     std::vector<double> vout;
0057     for (size_t i=0;i<v1.size();i++) vout.push_back(v1.at(i) + v2.at(i));
0058     return vout;
0059 }
0060
0061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0062
0063 ODESolver::ODESolver(): fNstepsForObserver(1)
0064 {}
0065
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0067
0068 double ODESolver::RungeKutta_Fehlberg( std::function<std::vector<double>(double,std::vector<double>)>
0069 func,std::vector<double> &y, double t, double stepsize)
0070 {
0071     //based on https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method
0072     const int nk=6;
0073     double h = stepsize;
0074     double CH[nk]={47./450.,0,12./25.,32./255.,1./30.,6./25.};
0075     double CT[nk]={-1./150.,0.,3./100.,-16./75.,-1./20.,6./25.};
0076     double A[nk]={0.,2./9.,1./3.,3./4.,1.,5./6.};
0077     double B21=2./9., B31=1./12., B41=69./128., B51=-17./12., B61=65./432.;
0078     double B32=1./4., B42=-243./128., B52=27./5., B62=13./16.;
0079     double B43=135./64., B53=-27./5., B63=13./16.;
0080     double B54=16./15., B64=4./27.;
0081     double B65=5./144.;
0082     double maxError = 1.;
0083     std::vector<double> k1 = func(t+A[0]*h,y)*h;
0084     std::vector<double> k2 = func(t+A[1]*h,y + k1*B21)*h;
0085     std::vector<double> k3 = func(t+A[2]*h,y + k1*B31 + k2*B32)*h;
0086     std::vector<double> k4 = func(t+A[3]*h,y + k1*B41 + k2*B42 + k3*B43)*h;
0087     std::vector<double> k5 = func(t+A[4]*h,y + k1*B51 + k2*B52 + k3*B53 + k4*B54)*h;
0088     std::vector<double> k6 = func(t+A[5]*h,y + k1*B61 + k2*B62 + k3*B63 + k4*B64 + k5*B65)*h;
0089     y = y + k1*CH[0] + k2*CH[1] + k3*CH[2] + k4*CH[3] + k5*CH[4] + k6*CH[5];
0090     auto TE = k1*CT[0] + k2*CT[1] + k3*CT[2] + k4*CT[3] + k5*CT[4] + k6*CT[5];
0091     absValuesVector(TE);
0092     maxError = *std::max_element(TE.begin(),TE.end());
0093
0094     k1.clear(); k1.shrink_to_fit();
0095     k2.clear(); k2.shrink_to_fit();
0096     k3.clear(); k3.shrink_to_fit();
0097     k4.clear(); k4.shrink_to_fit();
0098     k5.clear(); k5.shrink_to_fit();
0099     k6.clear(); k6.shrink_to_fit();
0100     TE.clear(); TE.shrink_to_fit();
0101     return  maxError;
0102 }
0103
0104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0105
0106 void ODESolver::Embedded_RungeKutta_Fehlberg(
0107     std::function<std::vector<double>(double,std::vector<double>)> func, std::vector<double> &y,
0108     double start,double end,double stepsize,double epsilon,
0109     std::vector<double> *time_observer,std::vector<std::vector<double>> *state_observer)
0110 {
0111     double t = start;
0112     double h = stepsize;
0113     int nsteps = 0;
0114     if (h < 0) h = (end -  start)/(10000.);
0115     if (time_observer) time_observer->push_back(t);
0116     if (state_observer) state_observer->push_back(y);
0117     auto ytemp = y;
0118     while (t < end)
0119     {
0120         ytemp = y;
0121         double maxerror = RungeKutta_Fehlberg(func,ytemp,t,h);
0122         double scale = 0.9*std::pow(epsilon/maxerror,1./5.);
0123
0124         double hnew = h*scale;
0125         while (maxerror > epsilon)
0126         {
0127             ytemp = y;
0128             maxerror = RungeKutta_Fehlberg(func,ytemp,t,hnew);
0129             scale = 0.9*std::pow(epsilon/maxerror,1./5.);
0130             hnew = hnew*scale;
0131         }
0132         h = hnew;
0133         y = ytemp;
0134         t += h;
0135         if (t > end) break;
0136         if ( time_observer || state_observer) {
0137             nsteps++;
0138             if (nsteps%fNstepsForObserver == 0) {
0139                 if (time_observer) time_observer->push_back(t);
0140                 if (state_observer) state_observer->push_back(y);
0141                 nsteps = 0;
0142             }
0143         }
0144     }
0145
0146     ytemp.clear(); ytemp.shrink_to_fit();
0147 }
0148
0149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0150
0151 void ODESolver::RungeKutta4(
0152     std::function<std::vector<double>(double,std::vector<double>)> func, std::vector<double> &y,
0153     double start,double end,double stepsize,
0154     std::vector<double> *time_observer,std::vector<std::vector<double>> *state_observer)
0155 {
0156     double t = start;
0157     double h = stepsize;
0158     int nsteps = 0;
0159     if (h < 0) h = (end -  start)/(10000.);
0160     if (time_observer) time_observer->push_back(t);
0161     if (state_observer) state_observer->push_back(y);
0162     while (t < end)
0163     {
0164         std::vector<double> k1 = func(t,y)*h;
0165         std::vector<double> k2 = func(t+0.5*h,y + k1*0.5)*h;
0166         std::vector<double> k3 = func(t+0.5*h,y + k2*0.5)*h;
0167         std::vector<double> k4 = func(t+h,y + k3)*h;
0168         t += h;
0169         if (t > end) break;
0170         y = y +(k1 +k2*2+k3*2+k4)*(1./6.0);
0171         if ( time_observer || state_observer) {
0172             nsteps++;
0173             if (nsteps%fNstepsForObserver == 0) {
0174                 if (time_observer) time_observer->push_back(t);
0175                 if (state_observer) state_observer->push_back(y);
0176                 nsteps = 0;
0177             }
0178         }
0179         k1.clear(); k1.shrink_to_fit();
0180         k2.clear(); k2.shrink_to_fit();
0181         k3.clear(); k3.shrink_to_fit();
0182         k4.clear(); k4.shrink_to_fit();
0183     }
0184 }
0185
0186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....