Flood: OrdinaryDifferentialEquations.h Source File

00001 /****************************************************************************************************************/
00002 /*                                                                                                              */
00003 /*   Flood: An Open Source Neural Networks C++ Library                                                          */
00004 /*   www.cimne.com/flood                                                                                        */
00005 /*                                                                                                              */
00006 /*   O R D I N A R Y   D I F F E R E N T I A L   E Q U A T I O N S   C L A S S   H E A D E R                    */
00007 /*                                                                                                              */
00008 /*   Roberto Lopez                                                                                              */
00009 /*   International Center for Numerical Methods in Engineering (CIMNE)                                          */
00010 /*   Technical University of Catalonia (UPC)                                                                    */
00011 /*   Barcelona, Spain                                                                                           */
00012 /*   E-mail: rlopez@cimne.upc.edu                                                                               */
00013 /*                                                                                                              */
00014 /****************************************************************************************************************/
00015 
00016 #ifndef __ORDINARYDIFFERENTIALEQUATIONS_H__
00017 #define __ORDINARYDIFFERENTIALEQUATIONS_H__
00018 
00019 #include <iostream>
00020 #include <fstream>
00021 #include <limits>
00022 #include <cmath>
00023 
00024 #include "Vector.h"
00025 #include "Matrix.h"
00026 
00027 namespace Flood
00028 {
00029 
00032 
00033 class OrdinaryDifferentialEquations 
00034 {
00035 public:
00036 
00037    // CONSTRUCTOR
00038 
00039    explicit OrdinaryDifferentialEquations(void);
00040 
00041 
00042    // DESTRUCTOR
00043 
00044    virtual ~OrdinaryDifferentialEquations(void);
00045 
00046 
00047    // METHODS
00048 
00049    int get_points_number(void);
00050 
00051    double get_tolerance(void);
00052    int get_initial_size(void);
00053    int get_warning_size(void);
00054    int get_error_size(void);
00055 
00056    bool get_display(void);
00057 
00058    void set_default(void);
00059 
00060    void set_points_number(int);
00061 
00062    void set_tolerance(double);
00063    void set_initial_size(int);
00064    void set_warning_size(int);
00065    void set_error_size(int);
00066 
00067    void set_display(bool);
00068 
00069    // C-style functions (state equations) 
00070 
00071    void calculate_Runge_Kutta_integral(
00072    Vector<double>&, Vector<double>&,
00073    double (*f)(double, double), double, double, double);
00074 
00075    void calculate_Runge_Kutta_integral(
00076    Vector<double>&, Vector<double>&, Vector<double>&, 
00077    double (*f_1)(double, double,double), 
00078    double (*f_2)(double, double,double), 
00079    double, double, double, double);
00080 
00081    void calculate_Runge_Kutta_integral
00082    (Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&, 
00083    double (*f_1)(double, double,double, double), 
00084    double (*f_2)(double, double,double, double),
00085    double (*f_3)(double, double,double, double),   
00086    double, double, double, double, double);
00087 
00088    void calculate_Runge_Kutta_integral
00089    (Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&, 
00090    double (*f_1)(double, double, double, double, double), 
00091    double (*f_2)(double, double, double, double, double),
00092    double (*f_3)(double, double, double, double, double),   
00093    double (*f_4)(double, double, double, double, double),   
00094    double, double, double, double, double, double);
00095 
00096    void calculate_Runge_Kutta_integral
00097    (Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&,  Vector<double>&, 
00098    double (*f_1)(double, double, double, double, double, double), 
00099    double (*f_2)(double, double, double, double, double, double),
00100    double (*f_3)(double, double, double, double, double, double),   
00101    double (*f_4)(double, double, double, double, double, double),   
00102    double (*f_5)(double, double, double, double, double, double),   
00103    double, double, double, double, double, double, double);
00104 
00105    int calculate_Runge_Kutta_Fehlberg_integral(Vector<double>&, Vector<double>&,
00106    double (*f)(double, double), double, double, double);
00107 
00108    int calculate_Runge_Kutta_Fehlberg_integral
00109    (Vector<double>&, Vector<double>&, Vector<double>&, 
00110    double (*f_1)(double, double, double), 
00111    double (*f_2)(double, double, double), 
00112    double, double, double, double);
00113 
00114    int calculate_Runge_Kutta_Fehlberg_integral
00115    (Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&, 
00116    double (*f_1)(double, double, double, double), 
00117    double (*f_2)(double, double, double, double), 
00118    double (*f_3)(double, double, double, double), 
00119    double, double, double, double, double);
00120 
00121    int calculate_Runge_Kutta_Fehlberg_integral
00122    (Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&,
00123    double (*f_1)(double, double, double, double, double), 
00124    double (*f_2)(double, double, double, double, double), 
00125    double (*f_3)(double, double, double, double, double), 
00126    double (*f_4)(double, double, double, double, double), 
00127    double, double, double, double, double, double);
00128 
00129 
00130    // Class member methods (state equations)
00131 
00132    // template<class T> void calculate_Runge_Kutta_integral(T, 
00133    // Vector<double>&, Vector<double>&,  
00134    // double (T::*f)(double, double), 
00135    // double, double, double, int) method
00136 
00153 
00154    template<class T> void calculate_Runge_Kutta_integral(T& t, 
00155    Vector<double>& x, Vector<double>& y,  
00156    double (T::*f)(double, double), 
00157    double a, double b, double ya) 
00158    {
00159       // Integration step
00160 
00161       double h = (b-a)/(points_number-1.0);      
00162 
00163       // Fourth order Runge-Kutta coefficients
00164 
00165       Vector<double> k(4);
00166 
00167       // Independent variable
00168 
00169       x[0] = a;
00170 
00171       for(int i = 0;  i < points_number-2; i++)
00172       {
00173          x[i+1] = x[i] + h;
00174      }
00175 
00176       x[points_number-1] = b;
00177 
00178       // Dependent variable
00179 
00180       y[0] = ya;
00181 
00182       for(int i = 0;  i < points_number-1; i++)
00183       {
00184          // Runge-Kutta coefficients
00185 
00186          k[0] = (t.*f)(x[i], y[i]);
00187          k[1] = (t.*f)(x[i]+h/2.0, y[i]+h*k[0]/2.0);
00188          k[2] = (t.*f)(x[i]+h/2.0, y[i]+h*k[1]/2.0);
00189          k[3] = (t.*f)(x[i+1], y[i]+h*k[2]);
00190 
00191          // Dependent variable derivative
00192 
00193          y[i+1] = y[i] + h*(k[0] + 2.0*k[1] + 2.0*k[2] + k[3])/6.0;
00194       }
00195    }
00196 
00197 
00198    // template<class T> void calculate_Runge_Kutta_integral(T, 
00199    // Vector<double>&, Vector<double>&, Vector<double>&, 
00200    // double (T::*f_1)(double, double, double), 
00201    // double (T::*f_2)(double, double, double), 
00202    // double, double, double, double, int) method 
00203 
00225 
00226    template<class T> void calculate_Runge_Kutta_integral(T& t, 
00227    Vector<double>& x, Vector<double>& y_1, Vector<double>& y_2, 
00228    double (T::*f_1)(double, double, double), 
00229    double (T::*f_2)(double, double, double), 
00230    double a, double b, double y1a, double y2a) 
00231    {
00232       // Integration step
00233 
00234       double h = (b-a)/(points_number-1.0);      
00235 
00236       // Fourth order Runge-Kutta coefficients
00237 
00238       Vector<double> k1(4);
00239       Vector<double> k2(4);
00240 
00241       // Independent variable
00242 
00243       x[0] = a;
00244 
00245       for(int i = 0; i < points_number-2; i++)
00246       {
00247          x[i+1] = x[i] + h;
00248      }
00249 
00250       x[points_number-1] = b;
00251 
00252       // Dependent variables
00253 
00254       y_1[0] = y1a;
00255       y_2[0] = y2a;
00256 
00257       for(int i = 0; i < points_number-1; i++)
00258       {
00259          // Runge-Kutta coefficients
00260 
00261          k1[0] = (t.*f_1)(x[i], y_1[i], y_2[i]);
00262          k2[0] = (t.*f_2)(x[i], y_1[i], y_2[i]);
00263 
00264          k1[1] = (t.*f_1)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0);
00265          k2[1] = (t.*f_2)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0);
00266 
00267          k1[2] = (t.*f_1)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0);
00268          k2[2] = (t.*f_2)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0);
00269 
00270          k1[3] = (t.*f_1)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2]);
00271          k2[3] = (t.*f_2)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2]);
00272 
00273          // Dependent variables
00274 
00275          y_1[i+1] = y_1[i] + h*(k1[0] + 2.0*k1[1] + 2.0*k1[2] + k1[3])/6.0;
00276          y_2[i+1] = y_2[i] + h*(k2[0] + 2.0*k2[1] + 2.0*k2[2] + k2[3])/6.0;
00277       }
00278    }
00279 
00280 
00281    // template<class T> void calculate_Runge_Kutta_integral(T, 
00282    // Vector<double>&, 
00283    // Vector<double>&, Vector<double>&, Vector<double>&, 
00284    // double (T::*f_1)(double, double, double, double), 
00285    // double (T::*f_2)(double, double, double, double), 
00286    // double (T::*f_3)(double, double, double, double), 
00287    // double, double, double, double, double, int) method
00288 
00314 
00315    template<class T> void calculate_Runge_Kutta_integral(T& t, 
00316    Vector<double>& x, Vector<double>& y_1, Vector<double>& y_2, Vector<double>& y_3, 
00317    double (T::*f_1)(double, double, double, double), 
00318    double (T::*f_2)(double, double, double, double), 
00319    double (T::*f_3)(double, double, double, double), 
00320    double a, double b, double y1a, double y2a, double y3a)
00321    {
00322       // Integration step
00323 
00324       double h = (b-a)/(points_number-1.0);      
00325 
00326       // Fourth order Runge-Kutta coefficients
00327 
00328       Vector<double> k1(4);
00329       Vector<double> k2(4);
00330       Vector<double> k3(4);
00331 
00332       // Independent variable
00333 
00334       x[0] = a;
00335 
00336       for(int i = 0; i < points_number-2; i++)
00337       {
00338          x[i+1] = x[i] + h;
00339       }
00340 
00341       x[points_number-1] = b;
00342 
00343       // Dependent variables
00344 
00345       y_1[0] = y1a;
00346       y_2[0] = y2a;
00347       y_3[0] = y3a;
00348 
00349 
00350       for(int i = 0;  i < points_number-1; i++)
00351       {
00352          // Runge-Kutta coefficients
00353 
00354          k1[0] = (t.*f_1)(x[i], y_1[i], y_2[i], y_3[i]);
00355          k2[0] = (t.*f_2)(x[i], y_1[i], y_2[i], y_3[i]);
00356          k3[0] = (t.*f_3)(x[i], y_1[i], y_2[i], y_3[i]);
00357 
00358          k1[1] = (t.*f_1)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0);
00359          k2[1] = (t.*f_2)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0);
00360          k3[1] = (t.*f_3)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0);
00361 
00362          k1[2] = (t.*f_1)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0);
00363          k2[2] = (t.*f_2)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0);
00364          k3[2] = (t.*f_3)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0);
00365 
00366          k1[3] = (t.*f_1)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2]);
00367          k2[3] = (t.*f_2)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2]);
00368          k3[3] = (t.*f_3)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2]);
00369 
00370          // Dependent variables
00371 
00372          y_1[i+1] = y_1[i] + h*(k1[0] + 2.0*k1[1] + 2.0*k1[2] + k1[3])/6.0;
00373          y_2[i+1] = y_2[i] + h*(k2[0] + 2.0*k2[1] + 2.0*k2[2] + k2[3])/6.0;
00374          y_3[i+1] = y_3[i] + h*(k3[0] + 2.0*k3[1] + 2.0*k3[2] + k3[3])/6.0;
00375       }
00376    }
00377 
00378 
00379    // template<class T> void calculate_Runge_Kutta_integral(T, 
00380    // Vector<double>&, 
00381    // Vector<double>&, Vector<double>&, 
00382    // Vector<double>&, Vector<double>&, 
00383    // double (T::*f_1)(double, double, double, double, double), 
00384    // double (T::*f_2)(double, double, double, double, double), 
00385    // double (T::*f_3)(double, double, double, double, double), 
00386    // double (T::*f_4)(double, double, double, double, double), 
00387    // double, double, double, double, double, double, int) method
00388 
00418    
00419    template<class T> void calculate_Runge_Kutta_integral(T& t, 
00420    Vector<double>& x, 
00421    Vector<double>& y_1, 
00422    Vector<double>& y_2, 
00423    Vector<double>& y_3, 
00424    Vector<double>& y_4, 
00425    double (T::*f_1)(double, double, double, double, double), 
00426    double (T::*f_2)(double, double, double, double, double), 
00427    double (T::*f_3)(double, double, double, double, double), 
00428    double (T::*f_4)(double, double, double, double, double), 
00429    double a, double b, double y1a, double y2a, double y3a, double y4a)
00430    {
00431       // Integration step
00432 
00433       double h = (b-a)/(points_number-1.0);      
00434 
00435       // Fourth order Runge-Kutta coefficients
00436 
00437       Vector<double> k1(4);
00438       Vector<double> k2(4);
00439       Vector<double> k3(4);
00440       Vector<double> k4(4);
00441 
00442       // Independent variable
00443 
00444       x[0] = a;
00445 
00446       for(int i = 0; i < points_number-2; i++)
00447       {
00448          x[i+1] = x[i] + h;
00449       }
00450 
00451       x[points_number-1] = b;
00452 
00453       // Dependent variables
00454 
00455       y_1[0] = y1a;
00456       y_2[0] = y2a;
00457       y_3[0] = y3a;
00458       y_4[0] = y4a;
00459 
00460       for(int i = 0;  i < points_number-1; i++)
00461       {
00462          // Runge-Kutta coefficients
00463 
00464          k1[0] = (t.*f_1)(x[i], y_1[i], y_2[i], y_3[i], y_4[i]);
00465          k2[0] = (t.*f_2)(x[i], y_1[i], y_2[i], y_3[i], y_4[i]);
00466          k3[0] = (t.*f_3)(x[i], y_1[i], y_2[i], y_3[i], y_4[i]);
00467          k4[0] = (t.*f_4)(x[i], y_1[i], y_2[i], y_3[i], y_4[i]);
00468 
00469          k1[1] = (t.*f_1)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0, y_4[i]+h*k4[0]/2.0);
00470          k2[1] = (t.*f_2)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0, y_4[i]+h*k4[0]/2.0);
00471          k3[1] = (t.*f_3)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0, y_4[i]+h*k4[0]/2.0);
00472          k4[1] = (t.*f_4)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0, y_4[i]+h*k4[0]/2.0);
00473 
00474          k1[2] = (t.*f_1)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0, y_4[i]+h*k4[1]/2.0);
00475          k2[2] = (t.*f_2)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0, y_4[i]+h*k4[1]/2.0);
00476          k3[2] = (t.*f_3)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0, y_4[i]+h*k4[1]/2.0);
00477          k4[2] = (t.*f_4)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0, y_4[i]+h*k4[1]/2.0);
00478 
00479          k1[3] = (t.*f_1)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2], y_4[i]+h*k4[2]);
00480          k2[3] = (t.*f_2)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2], y_4[i]+h*k4[2]);
00481          k3[3] = (t.*f_3)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2], y_4[i]+h*k4[2]);
00482          k4[3] = (t.*f_4)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2], y_4[i]+h*k4[2]);
00483 
00484          // Dependent variables
00485 
00486          y_1[i+1] = y_1[i] + h*(k1[0] + 2.0*k1[1] + 2.0*k1[2] + k1[3])/6.0;
00487          y_2[i+1] = y_2[i] + h*(k2[0] + 2.0*k2[1] + 2.0*k2[2] + k2[3])/6.0;
00488          y_3[i+1] = y_3[i] + h*(k3[0] + 2.0*k3[1] + 2.0*k3[2] + k3[3])/6.0;
00489          y_4[i+1] = y_4[i] + h*(k4[0] + 2.0*k4[1] + 2.0*k4[2] + k4[3])/6.0;
00490       }   
00491    }
00492 
00493 
00494    // template<class T> void calculate_Runge_Kutta_integral(T, 
00495    // Vector<double>&, Vector<double>&, Vector<double>&, 
00496    // Vector<double>&, Vector<double>&, Vector<double>&, 
00497    // double (T::*f_1)(double, double, double, double, double, double), 
00498    // double (T::*f_2)(double, double, double, double, double, double), 
00499    // double (T::*f_3)(double, double, double, double, double, double), 
00500    // double (T::*f_4)(double, double, double, double, double, double), 
00501    // double (T::*f_5)(double, double, double, double, double, double), 
00502    // double, double, double, double, double, double, double, int) method
00503 
00537    
00538    template<class T> void calculate_Runge_Kutta_integral(T& t, 
00539    Vector<double>& x, 
00540    Vector<double>& y_1, 
00541    Vector<double>& y_2, 
00542    Vector<double>& y_3, 
00543    Vector<double>& y_4, 
00544    Vector<double>& y_5, 
00545    double (T::*f_1)(double, double, double, double, double, double), 
00546    double (T::*f_2)(double, double, double, double, double, double), 
00547    double (T::*f_3)(double, double, double, double, double, double), 
00548    double (T::*f_4)(double, double, double, double, double, double), 
00549    double (T::*f_5)(double, double, double, double, double, double), 
00550    double a, double b, 
00551    double y1a, double y2a, double y3a, double y4a, double y5a)
00552    {
00553       // Integration step
00554 
00555       double h = (b-a)/(points_number-1.0);      
00556 
00557       // Fourth order Runge-Kutta coefficients
00558 
00559       Vector<double> k1(4);
00560       Vector<double> k2(4);
00561       Vector<double> k3(4);
00562       Vector<double> k4(4);
00563       Vector<double> k5(4);
00564 
00565       // Independent variable
00566 
00567       x[0] = a;
00568 
00569       for(int i = 0; i < points_number-2; i++)
00570       {
00571          x[i+1] = x[i] + h;
00572       }
00573 
00574       x[points_number-1] = b;
00575 
00576       // Dependent variables
00577 
00578       y_1[0] = y1a;
00579       y_2[0] = y2a;
00580       y_3[0] = y3a;
00581       y_4[0] = y4a;
00582       y_5[0] = y5a;
00583 
00584       for(int i = 0;  i < points_number-1; i++)
00585       {     
00586          // Runge-Kutta coefficients
00587 
00588          k1[0] = (t.*f_1)(x[i], y_1[i], y_2[i], y_3[i], y_4[i], y_5[i]);
00589          k2[0] = (t.*f_2)(x[i], y_1[i], y_2[i], y_3[i], y_4[i], y_5[i]);
00590          k3[0] = (t.*f_3)(x[i], y_1[i], y_2[i], y_3[i], y_4[i], y_5[i]);
00591          k4[0] = (t.*f_4)(x[i], y_1[i], y_2[i], y_3[i], y_4[i], y_5[i]);
00592          k5[0] = (t.*f_5)(x[i], y_1[i], y_2[i], y_3[i], y_4[i], y_5[i]);
00593 
00594          k1[1] = (t.*f_1)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0, y_4[i]+h*k4[0]/2.0, y_5[i]+h*k5[0]/2.0);
00595          k2[1] = (t.*f_2)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0, y_4[i]+h*k4[0]/2.0, y_5[i]+h*k5[0]/2.0);
00596          k3[1] = (t.*f_3)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0, y_4[i]+h*k4[0]/2.0, y_5[i]+h*k5[0]/2.0);
00597          k4[1] = (t.*f_4)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0, y_4[i]+h*k4[0]/2.0, y_5[i]+h*k5[0]/2.0);
00598          k5[1] = (t.*f_5)(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0, y_3[i]+h*k3[0]/2.0, y_4[i]+h*k4[0]/2.0, y_5[i]+h*k5[0]/2.0);
00599 
00600          k1[2] = (t.*f_1)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0, y_4[i]+h*k4[1]/2.0, y_5[i]+h*k5[1]/2.0);
00601          k2[2] = (t.*f_2)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0, y_4[i]+h*k4[1]/2.0, y_5[i]+h*k5[1]/2.0);
00602          k3[2] = (t.*f_3)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0, y_4[i]+h*k4[1]/2.0, y_5[i]+h*k5[1]/2.0);
00603          k4[2] = (t.*f_4)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0, y_4[i]+h*k4[1]/2.0, y_5[i]+h*k5[1]/2.0);
00604          k5[2] = (t.*f_5)(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0, y_3[i]+h*k3[1]/2.0, y_4[i]+h*k4[1]/2.0, y_5[i]+h*k5[1]/2.0);
00605 
00606          k1[3] = (t.*f_1)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2], y_4[i]+h*k4[2], y_5[i]+h*k5[2]);
00607          k2[3] = (t.*f_2)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2], y_4[i]+h*k4[2], y_5[i]+h*k5[2]);
00608          k3[3] = (t.*f_3)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2], y_4[i]+h*k4[2], y_5[i]+h*k5[2]);
00609          k4[3] = (t.*f_4)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2], y_4[i]+h*k4[2], y_5[i]+h*k5[2]);
00610          k5[3] = (t.*f_5)(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2], y_4[i]+h*k4[2], y_5[i]+h*k5[2]);
00611 
00612          // Dependent variables
00613 
00614          y_1[i+1] = y_1[i] + h*(k1[0] + 2.0*k1[1] + 2.0*k1[2] + k1[3])/6.0;
00615          y_2[i+1] = y_2[i] + h*(k2[0] + 2.0*k2[1] + 2.0*k2[2] + k2[3])/6.0;
00616          y_3[i+1] = y_3[i] + h*(k3[0] + 2.0*k3[1] + 2.0*k3[2] + k3[3])/6.0;
00617          y_4[i+1] = y_4[i] + h*(k4[0] + 2.0*k4[1] + 2.0*k4[2] + k4[3])/6.0;
00618          y_5[i+1] = y_5[i] + h*(k5[0] + 2.0*k5[1] + 2.0*k5[2] + k5[3])/6.0;
00619       }
00620    }
00621 
00622    // template<class T> int calculate_Runge_Kutta_Fehlberg_integral(T, 
00623    // Vector<double>&, Vector<double>&, 
00624    // double (T::*f)(double, double), 
00625    // double, double, double, double, int) method
00626  
00643 
00644    template<class T> int calculate_Runge_Kutta_Fehlberg_integral(T& t, 
00645    Vector<double>& x_out, Vector<double>& yOut,  
00646    double (T::*f)(double, double), 
00647    double a, double b, double ya)
00648    {
00649       int count = 0;
00650       int size = initial_size;
00651 
00652       Vector<double> tempX(size);
00653       Vector<double> tempY(size);
00654 
00655       tempX[0] = a;
00656       tempY[0] = ya;
00657 
00658       count++;
00659 
00660       double epsilon = std::numeric_limits<double>::epsilon();
00661 
00662       double a2 = 1.0/5.0;
00663       double a3 = 3.0/10.0;
00664       double a4 = 3.0/5.0;
00665       double a5 = 1.0;
00666       double a6 = 7.0/8.0;
00667 
00668       double b21 = 1.0/5.0;   
00669       double b31 = 3.0/40.0;
00670       double b32 = 9.0/40.0;   
00671       double b41 = 3.0/10.0;
00672       double b42 = -9.0/10.0;   
00673       double b43 = 6.0/5.0; 
00674       double b51 = -11.0/54.0;
00675       double b52 = 5.0/2.0;
00676       double b53 = -70.0/27.0;
00677       double b54 = 35.0/27.0; 
00678       double b61 = 1631.0/55296.0;
00679       double b62 = 175.0/512.0;   
00680       double b63 = 575.0/13824.0;   
00681       double b64 = 44275.0/110592.0;   
00682       double b65 = 253.0/4096.0;
00683 
00684       double c41 = 37.0/378.0;
00685       double c42 = 0.0;
00686       double c43 = 250.0/621.0;
00687       double c44 = 125.0/594.0;
00688       double c45 = 0.0;
00689       double c46 = 512.0/1771.0;
00690 
00691       double c51 = 2825.0/27648.0;  
00692       double c52 = 0.0;  
00693       double c53 = 18575.0/48384.0;   
00694       double c54 = 13525.0/55296.0;
00695       double c55 = 277.0/14336.0;
00696       double c56 = 1.0/4.0;
00697 
00698       double d1 = c41 - c51;
00699       double d2 = c42 - c52;
00700       double d3 = c43 - c53;
00701       double d4 = c44 - c54;
00702       double d5 = c45 - c55;
00703       double d6 = c46 - c56;
00704 
00705       double k11 = 0.0;
00706       double k12 = 0.0;
00707       double k13 = 0.0;
00708       double k14 = 0.0;
00709       double k15 = 0.0;
00710       double k16 = 0.0;
00711 
00712       // Main loop
00713 
00714       double x = a;
00715       double y = ya;
00716 
00717       double error = 0.0;
00718 
00719       double hmin = 16.0*epsilon*fabs(x);
00720 
00721       double h = (b-a)*1.0e-3;
00722 
00723       while(x < b)
00724       {   
00725          // Set smallest allowable stepsize
00726 
00727          hmin = 32.0*epsilon*fabs(x);
00728 
00729          if(h < hmin)
00730          {
00731             if(display)
00732          {
00733                std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
00734                          << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
00735                          << "Step size is less than smallest allowable." << std::endl;
00736          }
00737 
00738             h = hmin;
00739          }
00740 
00741          if(x + h > b)
00742          {
00743             h = b - x;
00744          }
00745 
00746          // Runge-Kutta-Fehlberg coefficients
00747 
00748          k11 = h*(t.*f)(x, y);
00749 
00750          k12 = h*(t.*f)(x+a2*h, y+b21*k11);
00751 
00752          k13 = h*(t.*f)(x+a3*h, y+b31*k11+b32*k12);
00753 
00754          k14 = h*(t.*f)(x+a4*h, y+b41*k11+b42*k12+b43*k13); 
00755 
00756          k15 = h*(t.*f)(x+a5*h, y+b51*k11+b52*k12+b53*k13+b54*k14); 
00757 
00758          k16 = h*(t.*f)(x+a6*h, y+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15);
00759 
00760          // Error estimate
00761 
00762          error = fabs(d1*k11+d2*k12+d3*k13+d4*k14+d5*k15+d6*k16);
00763 
00764          if(error <= tolerance) 
00765          {
00766             // Successful step
00767 
00768             x += h;
00769 
00770             y += c51*k11+c52*k12+c53*k13+c54*k14+c55*k15+c56*k16;
00771 
00772             tempX[count] = x;
00773             tempY[count] = y;
00774 
00775             count++;
00776 
00777             if(error != 0.0)
00778             {
00779                h *= 0.9*pow(fabs(tolerance/error),0.2);
00780             }
00781 
00782             // Computer's memory stuff
00783 
00784             if(count >= size)
00785             {
00786                size *= 2;
00787 
00788                // Control sentence
00789                               
00790                if(size > error_size)
00791                {
00792                   std::cerr << "Flood Error: OrdinaryDifferentialEquations class." << std::endl
00793                             << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
00794                             << "Step size is bigger than greatest allowable." << std::endl;
00795                             
00796                   exit(1);
00797                }
00798                else if(display && (size > warning_size))
00799                {
00800                   std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
00801                             << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
00802                             << "Size is " << size << std::endl;
00803                }
00804                 
00805                tempX.resize(size); 
00806                tempY.resize(size);               
00807             }
00808          }
00809          else
00810          {
00811             // Failed step 
00812 
00813             h *= 0.9*pow(fabs(tolerance/error),0.25);
00814          }
00815 
00816       } // end while loop    
00817 
00818       Vector<double> new_x(count);
00819       Vector<double> newY(count);
00820 
00821       for(int i = 0; i < count; i++)
00822       {
00823          new_x[i] = tempX[i];
00824          newY[i] = tempY[i];
00825       }
00826 
00827       x_out = new_x;
00828       yOut = newY;
00829 
00830       return(count);
00831    }
00832 
00833 
00834    // template<class T> int calculate_Runge_Kutta_Fehlberg_integral(T, 
00835    // Vector<double>&, 
00836    // Vector<double>&, Vector<double>&, 
00837    // Vector<double>&, Vector<double>&, 
00838    // double (T::*f_1)(double, double, double, double, double), 
00839    // double (T::*f_2)(double, double, double, double, double), 
00840    // double, double, double, double, double, double, int) method
00841  
00863 
00864    template<class T> int calculate_Runge_Kutta_Fehlberg_integral(T& t, 
00865    Vector<double>& x_out, Vector<double>& y1_out, Vector<double>& y2_out, 
00866    double (T::*f_1)(double, double, double), 
00867    double (T::*f_2)(double, double, double), 
00868    double a, double b, double y1a, double y2a)
00869    {
00870       double error = 0.0, error1 = 0.0, error2 = 0.0;
00871 
00872       int count = 0;
00873       int size = initial_size;
00874 
00875       Vector<double> tempX(size);
00876       Vector<double> tempY1(size);
00877       Vector<double> tempY2(size);
00878 
00879       tempX[0] = a;
00880       tempY1[0] = y1a;
00881       tempY2[0] = y2a;
00882       count++;
00883 
00884       double epsilon = std::numeric_limits<double>::epsilon();
00885 
00886       double a2 = 1.0/5.0;
00887       double a3 = 3.0/10.0;
00888       double a4 = 3.0/5.0;
00889       double a5 = 1.0;
00890       double a6 = 7.0/8.0;
00891 
00892       double b21 = 1.0/5.0;   
00893       double b31 = 3.0/40.0;
00894       double b32 = 9.0/40.0;   
00895       double b41 = 3.0/10.0;
00896       double b42 = -9.0/10.0;   
00897       double b43 = 6.0/5.0; 
00898       double b51 = -11.0/54.0;
00899       double b52 = 5.0/2.0;
00900       double b53 = -70.0/27.0;
00901       double b54 = 35.0/27.0; 
00902       double b61 = 1631.0/55296.0;
00903       double b62 = 175.0/512.0;   
00904       double b63 = 575.0/13824.0;   
00905       double b64 = 44275.0/110592.0;   
00906       double b65 = 253.0/4096.0;
00907 
00908       double c41 = 37.0/378.0;
00909       double c42 = 0.0;
00910       double c43 = 250.0/621.0;
00911       double c44 = 125.0/594.0;
00912       double c45 = 0.0;
00913       double c46 = 512.0/1771.0;
00914 
00915       double c51 = 2825.0/27648.0;  
00916       double c52 = 0.0;  
00917       double c53 = 18575.0/48384.0;   
00918       double c54 = 13525.0/55296.0;
00919       double c55 = 277.0/14336.0;
00920       double c56 = 1.0/4.0;
00921 
00922       double d1 = c41 - c51;
00923       double d2 = c42 - c52;
00924       double d3 = c43 - c53;
00925       double d4 = c44 - c54;
00926       double d5 = c45 - c55;
00927       double d6 = c46 - c56;
00928 
00929       double k11 = 0.0, k21 = 0.0;
00930       double k12 = 0.0, k22 = 0.0;
00931       double k13 = 0.0, k23 = 0.0;
00932       double k14 = 0.0, k24 = 0.0;
00933       double k15 = 0.0, k25 = 0.0;
00934       double k16 = 0.0, k26 = 0.0;
00935 
00936       // Main loop
00937 
00938       double x = a;
00939       double y_1 = y1a;
00940       double y_2 = y2a;
00941 
00942       double hmin = 16.0*epsilon*fabs(x);
00943 
00944       double h = (b-a)*1.0e-3;
00945 
00946       while(x < b)
00947       {   
00948          // Set smallest allowable stepsize
00949 
00950          hmin = 32.0*epsilon*fabs(x);
00951 
00952          if(h < hmin)
00953          {
00954             std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
00955                       << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
00956                       << "Step size is less than smallest allowable." << std::endl
00957                       << std::endl;
00958 
00959             h = hmin;
00960          }
00961 
00962          if(x + h > b)
00963          {
00964             h = b - x;
00965          }
00966 
00967          // Runge-Kutta-Fehlberg coefficients
00968 
00969          k11 = h*(t.*f_1)(x, y_1, y_2);
00970          k21 = h*(t.*f_2)(x, y_1, y_2);
00971 
00972          k12 = h*(t.*f_1)(x+a2*h, y_1+b21*k11, y_2+b21*k21);
00973          k22 = h*(t.*f_2)(x+a2*h, y_1+b21*k11, y_2+b21*k21);
00974 
00975          k13 = h*(t.*f_1)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22);
00976          k23 = h*(t.*f_2)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22);
00977 
00978          k14 = h*(t.*f_1)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23); 
00979          k24 = h*(t.*f_2)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23); 
00980 
00981          k15 = h*(t.*f_1)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14,
00982          y_2+b51*k21+b52*k22+b53*k23+b54*k24); 
00983          k25 = h*(t.*f_2)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, 
00984          y_2+b51*k21+b52*k22+b53*k23+b54*k24); 
00985 
00986          k16 = h*(t.*f_1)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, 
00987          y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25);
00988          k26 = h*(t.*f_2)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, 
00989          y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25);
00990 
00991          // Error estimate
00992 
00993          error1 = fabs(d1*k11+d2*k12+d3*k13+d4*k14+d5*k15+d6*k16);
00994          error2 = fabs(d1*k21+d2*k22+d3*k23+d4*k24+d5*k25+d6*k26);
00995 
00996          error = 0.0;
00997 
00998          if(error1 > error)
00999          {
01000             error = error1;
01001          }
01002          if(error2 > error)
01003          {
01004             error = error2;
01005          }
01006 
01007          if(error <= tolerance)
01008          {
01009             x += h;
01010 
01011             y_1 += c51*k11+c52*k12+c53*k13+c54*k14+c55*k15+c56*k16;
01012             y_2 += c51*k21+c52*k22+c53*k23+c54*k24+c55*k25+c56*k26;
01013 
01014             tempX[count] = x;
01015             tempY1[count] = y_1;
01016             tempY2[count] = y_2;
01017 
01018             count++;
01019 
01020             if(error != 0.0)
01021             {
01022                h *= 0.9*pow(fabs(tolerance/error),0.2);
01023             }
01024 
01025             if(count >= size)
01026             {
01027                size *= 2;
01028 
01029                if(size > error_size)
01030                {
01031                   std::cerr << "Flood Error: OrdinaryDifferentialEquations class." << std::endl
01032                             << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01033                             << "Step size is bigger than greatest allowable." << std::endl;
01034                             
01035                   exit(1);
01036                }
01037                else if(display && (size > warning_size))
01038                {
01039                   std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
01040                             << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01041                             << "Size is " << size << std::endl;
01042                }
01043 
01044                tempX.resize(size);
01045                tempY1.resize(size);
01046                tempY2.resize(size);
01047             }
01048          }
01049          else
01050          {
01051             h *= 0.9*pow(fabs(tolerance/error),0.25);
01052          }
01053       } // end while loop   
01054 
01055       Vector<double> new_x(count);
01056       Vector<double> new_y1(count);
01057       Vector<double> new_y2(count);
01058 
01059       for(int i = 0; i < count; i++)
01060       {
01061          new_x[i] = tempX[i];
01062          new_y1[i] = tempY1[i];
01063          new_y2[i] = tempY2[i];
01064       }
01065 
01066       x_out = new_x;
01067       y1_out = new_y1;
01068       y2_out = new_y2;
01069 
01070      return(count);
01071    }
01072 
01073 
01074    // template<class T> int calculate_Runge_Kutta_Fehlberg_integral(T, 
01075    // Vector<double>&, 
01076    // Vector<double>&, Vector<double>&, 
01077    // Vector<double>&, Vector<double>&, 
01078    // double (T::*f_1)(double, double, double, double, double), 
01079    // double (T::*f_2)(double, double, double, double, double), 
01080    // double (T::*f_3)(double, double, double, double, double), 
01081    // double, double, double, double, double, double, int) method
01082  
01108 
01109    template<class T> int calculate_Runge_Kutta_Fehlberg_integral(T& t, 
01110    Vector<double>& x_out, 
01111    Vector<double>& y1_out, 
01112    Vector<double>& y2_out, 
01113    Vector<double>& y3_out, 
01114    double (T::*f_1)(double, double, double, double), 
01115    double (T::*f_2)(double, double, double, double), 
01116    double (T::*f_3)(double, double, double, double), 
01117    double a, double b, 
01118    double y1a, double y2a, double y3a)
01119    {
01120       double error = 0.0, error1 = 0.0, error2 = 0.0;
01121 
01122       int count = 0;
01123 
01124       int size = initial_size;
01125 
01126       Vector<double> tempX(size);
01127       Vector<double> tempY1(size);
01128       Vector<double> tempY2(size);
01129 
01130       tempX[0] = a;
01131       tempY1[0] = y1a;
01132       tempY2[0] = y2a;
01133       count++;
01134 
01135       double epsilon = std::numeric_limits<double>::epsilon();
01136 
01137       double a2 = 1.0/5.0;
01138       double a3 = 3.0/10.0;
01139       double a4 = 3.0/5.0;
01140       double a5 = 1.0;
01141       double a6 = 7.0/8.0;
01142 
01143       double b21 = 1.0/5.0;   
01144       double b31 = 3.0/40.0;
01145       double b32 = 9.0/40.0;   
01146       double b41 = 3.0/10.0;
01147       double b42 = -9.0/10.0;   
01148       double b43 = 6.0/5.0; 
01149       double b51 = -11.0/54.0;
01150       double b52 = 5.0/2.0;
01151       double b53 = -70.0/27.0;
01152       double b54 = 35.0/27.0; 
01153       double b61 = 1631.0/55296.0;
01154       double b62 = 175.0/512.0;   
01155       double b63 = 575.0/13824.0;   
01156       double b64 = 44275.0/110592.0;   
01157       double b65 = 253.0/4096.0;
01158 
01159       double c41 = 37.0/378.0;
01160       double c42 = 0.0;
01161       double c43 = 250.0/621.0;
01162       double c44 = 125.0/594.0;
01163       double c45 = 0.0;
01164       double c46 = 512.0/1771.0;
01165 
01166       double c51 = 2825.0/27648.0;  
01167       double c52 = 0.0;  
01168       double c53 = 18575.0/48384.0;   
01169       double c54 = 13525.0/55296.0;
01170       double c55 = 277.0/14336.0;
01171       double c56 = 1.0/4.0;
01172 
01173       double d1 = c41 - c51;
01174       double d2 = c42 - c52;
01175       double d3 = c43 - c53;
01176       double d4 = c44 - c54;
01177       double d5 = c45 - c55;
01178       double d6 = c46 - c56;
01179 
01180       double k11 = 0.0, k21 = 0.0;
01181       double k12 = 0.0, k22 = 0.0;
01182       double k13 = 0.0, k23 = 0.0;
01183       double k14 = 0.0, k24 = 0.0;
01184       double k15 = 0.0, k25 = 0.0;
01185       double k16 = 0.0, k26 = 0.0;
01186 
01187       // Main loop
01188 
01189       double x = a;
01190       double y_1 = y1a;
01191       double y_2 = y2a;
01192 
01193       double hmin = 16.0*epsilon*fabs(x);
01194 
01195       double h = (b-a)*1.0e-3;
01196 
01197       while(x < b)
01198       {   
01199          // Set smallest allowable stepsize
01200 
01201          hmin = 32.0*epsilon*fabs(x);
01202 
01203          if(h < hmin)
01204          {
01205             std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
01206                       << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01207                       << "Step size is less than smallest allowable." << std::endl;
01208 
01209             h = hmin;
01210          }
01211 
01212          if(x + h > b)
01213          {
01214             h = b - x;
01215          }
01216 
01217          // Runge-Kutta-Fehlberg coefficients
01218 
01219          k11 = h*(t.*f_1)(x, y_1, y_2);
01220          k21 = h*(t.*f_2)(x, y_1, y_2);
01221 
01222          k12 = h*(t.*f_1)(x+a2*h, y_1+b21*k11, y_2+b21*k21);
01223          k22 = h*(t.*f_2)(x+a2*h, y_1+b21*k11, y_2+b21*k21);
01224 
01225          k13 = h*(t.*f_1)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22);
01226          k23 = h*(t.*f_2)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22);
01227 
01228          k14 = h*(t.*f_1)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23); 
01229          k24 = h*(t.*f_2)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23); 
01230 
01231          k15 = h*(t.*f_1)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24); 
01232          k25 = h*(t.*f_2)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24); 
01233 
01234          k16 = h*(t.*f_1)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25);
01235          k26 = h*(t.*f_2)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25);
01236 
01237          // Error estimate
01238 
01239          error1 = fabs(d1*k11+d2*k12+d3*k13+d4*k14+d5*k15+d6*k16);
01240          error2 = fabs(d1*k21+d2*k22+d3*k23+d4*k24+d5*k25+d6*k26);
01241 
01242          error = 0.0;
01243 
01244          if(error1 > error)
01245          {
01246             error = error1;
01247          }
01248          if(error2 > error)
01249          {
01250             error = error2;
01251          }
01252 
01253          if(error <= tolerance)
01254          {
01255             x += h;
01256 
01257             y_1 += c51*k11+c52*k12+c53*k13+c54*k14+c55*k15+c56*k16;
01258             y_2 += c51*k21+c52*k22+c53*k23+c54*k24+c55*k25+c56*k26;
01259 
01260             tempX[count] = x;
01261             tempY1[count] = y_1;
01262             tempY2[count] = y_2;
01263 
01264             count++;
01265 
01266             if(error != 0.0)
01267             {
01268                h *= 0.9*pow(fabs(tolerance/error),0.2);
01269             }
01270 
01271             if(count >= size)
01272             {
01273                size *= 2;
01274 
01275             if(size > error_size)
01276             {
01277                   std::cerr << "Flood Error: OrdinaryDifferentialEquations class." << std::endl
01278                             << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01279                             << "Step size is bigger than greatest allowable." << std::endl;
01280                             
01281                   exit(1);
01282             }
01283             else if(display && (size > warning_size))
01284             {
01285                   std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
01286                             << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01287                             << "Size is " << size << std::endl;
01288             }
01289 
01290             tempX.resize(size);
01291                tempY1.resize(size);
01292                tempY2.resize(size);
01293             }
01294          }
01295          else
01296          {
01297             h *= 0.9*pow(fabs(tolerance/error),0.25);
01298          }
01299       } // end while loop   
01300 
01301       Vector<double> new_x(count);
01302       Vector<double> new_y1(count);
01303       Vector<double> new_y2(count);
01304 
01305       for(int i = 0; i < count; i++)
01306       {
01307          new_x[i] = tempX[i];
01308          new_y1[i] = tempY1[i];
01309          new_y2[i] = tempY2[i];
01310       }
01311 
01312       x_out = new_x;
01313       y1_out = new_y1;
01314       y2_out = new_y2;
01315    }
01316 
01317 
01318    // template<class T> int calculate_Runge_Kutta_Fehlberg_integral(T& t, 
01319    // Vector<double>& x_out, 
01320    // Vector<double>& y1_out, 
01321    // Vector<double>& y2_out, 
01322    // Vector<double>& y3_out, 
01323    // Vector<double>& y4_out, 
01324    // double (T::*f_1)(double, double, double, double, double), 
01325    // double (T::*f_2)(double, double, double, double, double), 
01326    // double (T::*f_3)(double, double, double, double, double), 
01327    // double (T::*f_4)(double, double, double, double, double), 
01328    // double a, double b, double y1a, double y2a, double y3a, double y4a) method
01329 
01359 
01360    template<class T> int calculate_Runge_Kutta_Fehlberg_integral(T& t, 
01361    Vector<double>& x_out, 
01362    Vector<double>& y1_out, 
01363    Vector<double>& y2_out, 
01364    Vector<double>& y3_out, 
01365    Vector<double>& y4_out, 
01366    double (T::*f_1)(double, double, double, double, double), 
01367    double (T::*f_2)(double, double, double, double, double), 
01368    double (T::*f_3)(double, double, double, double, double), 
01369    double (T::*f_4)(double, double, double, double, double), 
01370    double a, double b, double y1a, double y2a, double y3a, double y4a)
01371    {
01372       double error = 0.0, error1 = 0.0, error2 = 0.0, error3 = 0.0, error4 = 0.0;
01373 
01374       int count = 0;
01375 
01376      int size = initial_size;
01377 
01378       Vector<double> tempX(size);
01379       Vector<double> tempY1(size);
01380       Vector<double> tempY2(size);
01381       Vector<double> tempY3(size);
01382       Vector<double> tempY4(size);
01383 
01384       tempX[0] = a;
01385 
01386       tempY1[0] = y1a;
01387       tempY2[0] = y2a;
01388       tempY3[0] = y3a;
01389       tempY4[0] = y4a;
01390 
01391       count++;
01392 
01393       double epsilon = std::numeric_limits<double>::epsilon();
01394 
01395       double a2 = 1.0/5.0;
01396       double a3 = 3.0/10.0;
01397       double a4 = 3.0/5.0;
01398       double a5 = 1.0;
01399       double a6 = 7.0/8.0;
01400 
01401       double b21 = 1.0/5.0;   
01402       double b31 = 3.0/40.0;
01403       double b32 = 9.0/40.0;   
01404       double b41 = 3.0/10.0;
01405       double b42 = -9.0/10.0;   
01406       double b43 = 6.0/5.0; 
01407       double b51 = -11.0/54.0;
01408       double b52 = 5.0/2.0;
01409       double b53 = -70.0/27.0;
01410       double b54 = 35.0/27.0; 
01411       double b61 = 1631.0/55296.0;
01412       double b62 = 175.0/512.0;   
01413       double b63 = 575.0/13824.0;   
01414       double b64 = 44275.0/110592.0;   
01415       double b65 = 253.0/4096.0;
01416 
01417       double c41 = 37.0/378.0;
01418       double c42 = 0.0;
01419       double c43 = 250.0/621.0;
01420       double c44 = 125.0/594.0;
01421       double c45 = 0.0;
01422       double c46 = 512.0/1771.0;
01423 
01424       double c51 = 2825.0/27648.0;  
01425       double c52 = 0.0;  
01426       double c53 = 18575.0/48384.0;   
01427       double c54 = 13525.0/55296.0;
01428       double c55 = 277.0/14336.0;
01429       double c56 = 1.0/4.0;
01430 
01431       double d1 = c41 - c51;
01432       double d2 = c42 - c52;
01433       double d3 = c43 - c53;
01434       double d4 = c44 - c54;
01435       double d5 = c45 - c55;
01436       double d6 = c46 - c56;
01437 
01438       double k11 = 0.0, k21 = 0.0, k31 = 0.0, k41 = 0.0;
01439       double k12 = 0.0, k22 = 0.0, k32 = 0.0, k42 = 0.0;
01440       double k13 = 0.0, k23 = 0.0, k33 = 0.0, k43 = 0.0;
01441       double k14 = 0.0, k24 = 0.0, k34 = 0.0, k44 = 0.0;
01442       double k15 = 0.0, k25 = 0.0, k35 = 0.0, k45 = 0.0;
01443       double k16 = 0.0, k26 = 0.0, k36 = 0.0, k46 = 0.0;
01444 
01445       // Main loop
01446 
01447       double x = a;
01448       double y_1 = y1a;
01449       double y_2 = y2a;
01450       double y_3 = y3a;
01451       double y_4 = y4a;
01452 
01453       double hmin = 16.0*epsilon*fabs(x);
01454 
01455       double h = (b-a)*1.0e-3;
01456 
01457       while(x < b)
01458       {
01459          // Set smallest allowable stepsize
01460 
01461          hmin = 32.0*epsilon*fabs(x);
01462 
01463          if(h < hmin)
01464          {
01465             std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
01466                       << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01467                       << "Step size is less than smallest allowable." << std::endl;
01468 
01469             h = hmin;
01470          }
01471 
01472          if(x + h > b)
01473          {
01474             h = b - x;
01475          }
01476 
01477          // Runge-Kutta-Fehlberg coefficients
01478 
01479          k11 = h*(t.*f_1)(x, y_1, y_2, y_3, y_4);
01480          k21 = h*(t.*f_2)(x, y_1, y_2, y_3, y_4);
01481          k31 = h*(t.*f_3)(x, y_1, y_2, y_3, y_4);
01482          k41 = h*(t.*f_4)(x, y_1, y_2, y_3, y_4);
01483 
01484          k12 = h*(t.*f_1)(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41);
01485          k22 = h*(t.*f_2)(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41);
01486          k32 = h*(t.*f_3)(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41);
01487          k42 = h*(t.*f_4)(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41);
01488 
01489          k13 = h*(t.*f_1)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32, y_4+b31*k41+b32*k42);
01490          k23 = h*(t.*f_2)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32, y_4+b31*k41+b32*k42);
01491          k33 = h*(t.*f_3)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32, y_4+b31*k41+b32*k42);
01492          k43 = h*(t.*f_4)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32, y_4+b31*k41+b32*k42);
01493 
01494          k14 = h*(t.*f_1)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23, y_3+b41*k31+b42*k32+b43*k33, y_4+b41*k41+b42*k42+b43*k43); 
01495          k24 = h*(t.*f_2)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23, y_3+b41*k31+b42*k32+b43*k33, y_4+b41*k41+b42*k42+b43*k43); 
01496          k34 = h*(t.*f_3)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23, y_3+b41*k31+b42*k32+b43*k33, y_4+b41*k41+b42*k42+b43*k43); 
01497          k44 = h*(t.*f_4)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23, y_3+b41*k31+b42*k32+b43*k33, y_4+b41*k41+b42*k42+b43*k43); 
01498 
01499          k15 = h*(t.*f_1)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24, y_3+b51*k31+b52*k32+b53*k33+b54*k34, y_4+b51*k41+b52*k42+b53*k43+b54*k44); 
01500          k25 = h*(t.*f_2)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24, y_3+b51*k31+b52*k32+b53*k33+b54*k34, y_4+b51*k41+b52*k42+b53*k43+b54*k44); 
01501          k35 = h*(t.*f_3)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24, y_3+b51*k31+b52*k32+b53*k33+b54*k34, y_4+b51*k41+b52*k42+b53*k43+b54*k44); 
01502          k45 = h*(t.*f_4)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24, y_3+b51*k31+b52*k32+b53*k33+b54*k34, y_4+b51*k41+b52*k42+b53*k43+b54*k44); 
01503 
01504          k16 = h*(t.*f_1)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25, y_3+b61*k31+b62*k32+b63*k33+b64*k34+b65*k35, y_4+b61*k41+b62*k42+b63*k43+b64*k44+b65*k45);
01505          k26 = h*(t.*f_2)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25, y_3+b61*k31+b62*k32+b63*k33+b64*k34+b65*k35, y_4+b61*k41+b62*k42+b63*k43+b64*k44+b65*k45);
01506          k36 = h*(t.*f_3)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25, y_3+b61*k31+b62*k32+b63*k33+b64*k34+b65*k35, y_4+b61*k41+b62*k42+b63*k43+b64*k44+b65*k45);
01507          k46 = h*(t.*f_4)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25, y_3+b61*k31+b62*k32+b63*k33+b64*k34+b65*k35, y_4+b61*k41+b62*k42+b63*k43+b64*k44+b65*k45);
01508 
01509          // Error estimate
01510 
01511          error1 = fabs(d1*k11+d2*k12+d3*k13+d4*k14+d5*k15+d6*k16);
01512          error2 = fabs(d1*k21+d2*k22+d3*k23+d4*k24+d5*k25+d6*k26);
01513          error3 = fabs(d1*k31+d2*k32+d3*k33+d4*k34+d5*k35+d6*k36);
01514          error4 = fabs(d1*k41+d2*k42+d3*k43+d4*k44+d5*k45+d6*k46);
01515 
01516          error = 0.0;
01517 
01518          if(error1 > error)
01519          {
01520             error = error1;
01521          }
01522          if(error2 > error)
01523          {
01524             error = error2;
01525          }
01526          if(error3 > error)
01527          {
01528             error = error3;
01529          }
01530          if(error4 > error)
01531          {
01532             error = error4;
01533          }
01534 
01535          if(error <= tolerance)
01536          {
01537             // Successful step
01538 
01539             x += h;
01540 
01541             y_1 += c51*k11+c52*k12+c53*k13+c54*k14+c55*k15+c56*k16;
01542             y_2 += c51*k21+c52*k22+c53*k23+c54*k24+c55*k25+c56*k26;
01543             y_3 += c51*k31+c52*k32+c53*k33+c54*k34+c55*k35+c56*k36;
01544             y_4 += c51*k41+c52*k42+c53*k43+c54*k44+c55*k45+c56*k46;
01545 
01546             tempX[count] = x;
01547             tempY1[count] = y_1;
01548             tempY2[count] = y_2;
01549             tempY3[count] = y_3;
01550             tempY4[count] = y_4;
01551 
01552             count++;
01553 
01554             if(error != 0.0)
01555             {
01556                h *= 0.9*pow(fabs(tolerance/error),0.2);
01557             }
01558 
01559             if(count >= size)
01560             {
01561                // Not enouth memory reserved
01562 
01563                size *= 2;
01564 
01565             if(size > error_size)
01566             {
01567                   std::cerr << "Flood Error: OrdinaryDifferentialEquations class." << std::endl
01568                             << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01569                             << "Step size is bigger than greatest allowable." << std::endl;
01570                             
01571                   exit(1);
01572             }
01573             else if(display && (size > warning_size))
01574             {
01575                   std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
01576                             << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01577                             << "Size is " << size << std::endl;
01578             }
01579 
01580                tempX.resize(size);
01581                tempY1.resize(size);
01582                tempY2.resize(size);
01583                tempY3.resize(size);
01584                tempY4.resize(size);
01585          }
01586          }
01587          else
01588          {
01589             // Failed step
01590 
01591             h *= 0.9*pow(fabs(tolerance/error),0.25);
01592          }
01593       } // end while loop   
01594 
01595       Vector<double> new_x(count);
01596       Vector<double> new_y1(count);
01597       Vector<double> new_y2(count);
01598       Vector<double> new_y3(count);
01599       Vector<double> new_y4(count);
01600 
01601       for(int i = 0; i < count; i++)
01602       {
01603          new_x[i] = tempX[i];
01604          new_y1[i] = tempY1[i];
01605          new_y2[i] = tempY2[i];
01606          new_y3[i] = tempY3[i];
01607          new_y4[i] = tempY4[i];
01608       }
01609 
01610       x_out = new_x;
01611       y1_out = new_y1;
01612       y2_out = new_y2;
01613       y3_out = new_y3;
01614       y4_out = new_y4;
01615 
01616       return(count);
01617    }
01618 
01619 
01620    // template<class T> int calculate_Runge_Kutta_Fehlberg_integral(T& t, 
01621    // Vector<double>& x_out, 
01622    // Vector<double>& y1_out, 
01623    // Vector<double>& y2_out, 
01624    // Vector<double>& y3_out, 
01625    // Vector<double>& y4_out, 
01626    // Vector<double>& y5_out, 
01627    // double (T::*f_1)(double, double, double, double, double, double), 
01628    // double (T::*f_2)(double, double, double, double, double, double), 
01629    // double (T::*f_3)(double, double, double, double, double, double), 
01630    // double (T::*f_4)(double, double, double, double, double, double), 
01631    // double (T::*f_5)(double, double, double, double, double, double), 
01632    // double a, double b, double y1a, double y2a, double y3a, double y4a, double y5a) method
01633 
01667 
01668    template<class T> int calculate_Runge_Kutta_Fehlberg_integral(T& t, 
01669    Vector<double>& x_out, 
01670    Vector<double>& y1_out, 
01671    Vector<double>& y2_out, 
01672    Vector<double>& y3_out, 
01673    Vector<double>& y4_out, 
01674    Vector<double>& y5_out, 
01675    double (T::*f_1)(double, double, double, double, double, double), 
01676    double (T::*f_2)(double, double, double, double, double, double), 
01677    double (T::*f_3)(double, double, double, double, double, double), 
01678    double (T::*f_4)(double, double, double, double, double, double), 
01679    double (T::*f_5)(double, double, double, double, double, double), 
01680    double a, double b, double y1a, double y2a, double y3a, double y4a, double y5a)
01681    {
01682       double error = 0.0;
01683 
01684      double error1 = 0.0, error2 = 0.0, error3 = 0.0, error4 = 0.0, error5 = 0.0;
01685 
01686       int count = 0;
01687 
01688      int size = initial_size;
01689 
01690       Vector<double> tempX(size);
01691       Vector<double> tempY1(size);
01692       Vector<double> tempY2(size);
01693       Vector<double> tempY3(size);
01694       Vector<double> tempY4(size);
01695       Vector<double> tempY5(size);
01696 
01697       tempX[0] = a;
01698 
01699       tempY1[0] = y1a;
01700       tempY2[0] = y2a;
01701       tempY3[0] = y3a;
01702       tempY4[0] = y4a;
01703       tempY5[0] = y5a;
01704 
01705       count++;
01706 
01707       double epsilon = std::numeric_limits<double>::epsilon();
01708 
01709       double a2 = 1.0/5.0;
01710       double a3 = 3.0/10.0;
01711       double a4 = 3.0/5.0;
01712       double a5 = 1.0;
01713       double a6 = 7.0/8.0;
01714 
01715       double b21 = 1.0/5.0;   
01716       double b31 = 3.0/40.0;
01717       double b32 = 9.0/40.0;   
01718       double b41 = 3.0/10.0;
01719       double b42 = -9.0/10.0;   
01720       double b43 = 6.0/5.0; 
01721       double b51 = -11.0/54.0;
01722       double b52 = 5.0/2.0;
01723       double b53 = -70.0/27.0;
01724       double b54 = 35.0/27.0; 
01725       double b61 = 1631.0/55296.0;
01726       double b62 = 175.0/512.0;   
01727       double b63 = 575.0/13824.0;   
01728       double b64 = 44275.0/110592.0;   
01729       double b65 = 253.0/4096.0;
01730 
01731       double c41 = 37.0/378.0;
01732       double c42 = 0.0;
01733       double c43 = 250.0/621.0;
01734       double c44 = 125.0/594.0;
01735       double c45 = 0.0;
01736       double c46 = 512.0/1771.0;
01737 
01738       double c51 = 2825.0/27648.0;  
01739       double c52 = 0.0;  
01740       double c53 = 18575.0/48384.0;   
01741       double c54 = 13525.0/55296.0;
01742       double c55 = 277.0/14336.0;
01743       double c56 = 1.0/4.0;
01744 
01745       double d1 = c41 - c51;
01746       double d2 = c42 - c52;
01747       double d3 = c43 - c53;
01748       double d4 = c44 - c54;
01749       double d5 = c45 - c55;
01750       double d6 = c46 - c56;
01751 
01752       double k11 = 0.0, k21 = 0.0, k31 = 0.0, k41 = 0.0, k51 = 0.0;
01753       double k12 = 0.0, k22 = 0.0, k32 = 0.0, k42 = 0.0, k52 = 0.0;
01754       double k13 = 0.0, k23 = 0.0, k33 = 0.0, k43 = 0.0, k53 = 0.0;
01755       double k14 = 0.0, k24 = 0.0, k34 = 0.0, k44 = 0.0, k54 = 0.0;
01756       double k15 = 0.0, k25 = 0.0, k35 = 0.0, k45 = 0.0, k55 = 0.0;
01757       double k16 = 0.0, k26 = 0.0, k36 = 0.0, k46 = 0.0, k56 = 0.0;
01758 
01759       // Main loop
01760 
01761       double x = a;
01762       double y_1 = y1a;
01763       double y_2 = y2a;
01764       double y_3 = y3a;
01765       double y_4 = y4a;
01766       double y_5 = y5a;
01767 
01768       double hmin = 16.0*epsilon*fabs(x);
01769 
01770       double h = (b-a)*1.0e-3;
01771 
01772       while(x < b)
01773       {
01774          // Set smallest allowable stepsize
01775 
01776          hmin = 32.0*epsilon*fabs(x);
01777 
01778          if(h < hmin)
01779          {
01780             std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
01781                       << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01782                       << "Step size is less than smallest allowable." << std::endl
01783                       << std::endl;
01784 
01785             h = hmin;
01786          }
01787 
01788          if(x + h > b)
01789          {
01790             h = b - x;
01791          }
01792 
01793          // Runge-Kutta-Fehlberg coefficients
01794 
01795          k11 = h*(t.*f_1)(x, y_1, y_2, y_3, y_4, y_5);
01796          k21 = h*(t.*f_2)(x, y_1, y_2, y_3, y_4, y_5);
01797          k31 = h*(t.*f_3)(x, y_1, y_2, y_3, y_4, y_5);
01798          k41 = h*(t.*f_4)(x, y_1, y_2, y_3, y_4, y_5);
01799          k51 = h*(t.*f_5)(x, y_1, y_2, y_3, y_4, y_5);
01800 
01801          k12 = h*(t.*f_1)(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41, y_5+b21*k51);
01802          k22 = h*(t.*f_2)(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41, y_5+b21*k51);
01803          k32 = h*(t.*f_3)(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41, y_5+b21*k51);
01804          k42 = h*(t.*f_4)(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41, y_5+b21*k51);
01805          k52 = h*(t.*f_5)(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41, y_5+b21*k51);
01806 
01807          k13 = h*(t.*f_1)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32, y_4+b31*k41+b32*k42, y_5+b31*k51+b32*k52);
01808          k23 = h*(t.*f_2)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32, y_4+b31*k41+b32*k42, y_5+b31*k51+b32*k52);
01809          k33 = h*(t.*f_3)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32, y_4+b31*k41+b32*k42, y_5+b31*k51+b32*k52);
01810          k43 = h*(t.*f_4)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32, y_4+b31*k41+b32*k42, y_5+b31*k51+b32*k52);
01811          k53 = h*(t.*f_5)(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32, y_4+b31*k41+b32*k42, y_5+b31*k51+b32*k52);
01812 
01813          k14 = h*(t.*f_1)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23, y_3+b41*k31+b42*k32+b43*k33, y_4+b41*k41+b42*k42+b43*k43, y_5+b41*k51+b42*k52+b43*k53); 
01814          k24 = h*(t.*f_2)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23, y_3+b41*k31+b42*k32+b43*k33, y_4+b41*k41+b42*k42+b43*k43, y_5+b41*k51+b42*k52+b43*k53); 
01815          k34 = h*(t.*f_3)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23, y_3+b41*k31+b42*k32+b43*k33, y_4+b41*k41+b42*k42+b43*k43, y_5+b41*k51+b42*k52+b43*k53); 
01816          k44 = h*(t.*f_4)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23, y_3+b41*k31+b42*k32+b43*k33, y_4+b41*k41+b42*k42+b43*k43, y_5+b41*k51+b42*k52+b43*k53); 
01817          k54 = h*(t.*f_5)(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23, y_3+b41*k31+b42*k32+b43*k33, y_4+b41*k41+b42*k42+b43*k43, y_5+b41*k51+b42*k52+b43*k53); 
01818 
01819          k15 = h*(t.*f_1)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24, y_3+b51*k31+b52*k32+b53*k33+b54*k34, y_4+b51*k41+b52*k42+b53*k43+b54*k44, y_5+b51*k51+b52*k52+b53*k53+b54*k54); 
01820          k25 = h*(t.*f_2)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24, y_3+b51*k31+b52*k32+b53*k33+b54*k34, y_4+b51*k41+b52*k42+b53*k43+b54*k44, y_5+b51*k51+b52*k52+b53*k53+b54*k54); 
01821          k35 = h*(t.*f_3)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24, y_3+b51*k31+b52*k32+b53*k33+b54*k34, y_4+b51*k41+b52*k42+b53*k43+b54*k44, y_5+b51*k51+b52*k52+b53*k53+b54*k54); 
01822          k45 = h*(t.*f_4)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24, y_3+b51*k31+b52*k32+b53*k33+b54*k34, y_4+b51*k41+b52*k42+b53*k43+b54*k44, y_5+b51*k51+b52*k52+b53*k53+b54*k54); 
01823          k55 = h*(t.*f_5)(x+a5*h, y_1+b51*k11+b52*k12+b53*k13+b54*k14, y_2+b51*k21+b52*k22+b53*k23+b54*k24, y_3+b51*k31+b52*k32+b53*k33+b54*k34, y_4+b51*k41+b52*k42+b53*k43+b54*k44, y_5+b51*k51+b52*k52+b53*k53+b54*k54); 
01824 
01825          k16 = h*(t.*f_1)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25, y_3+b61*k31+b62*k32+b63*k33+b64*k34+b65*k35, y_4+b61*k41+b62*k42+b63*k43+b64*k44+b65*k45, y_5+b61*k51+b62*k52+b63*k53+b64*k54+b65*k55);
01826          k26 = h*(t.*f_2)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25, y_3+b61*k31+b62*k32+b63*k33+b64*k34+b65*k35, y_4+b61*k41+b62*k42+b63*k43+b64*k44+b65*k45, y_5+b61*k51+b62*k52+b63*k53+b64*k54+b65*k55);
01827          k36 = h*(t.*f_3)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25, y_3+b61*k31+b62*k32+b63*k33+b64*k34+b65*k35, y_4+b61*k41+b62*k42+b63*k43+b64*k44+b65*k45, y_5+b61*k51+b62*k52+b63*k53+b64*k54+b65*k55);
01828          k46 = h*(t.*f_4)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25, y_3+b61*k31+b62*k32+b63*k33+b64*k34+b65*k35, y_4+b61*k41+b62*k42+b63*k43+b64*k44+b65*k45, y_5+b61*k51+b62*k52+b63*k53+b64*k54+b65*k55);
01829          k56 = h*(t.*f_5)(x+a6*h, y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25, y_3+b61*k31+b62*k32+b63*k33+b64*k34+b65*k35, y_4+b61*k41+b62*k42+b63*k43+b64*k44+b65*k45, y_5+b61*k51+b62*k52+b63*k53+b64*k54+b65*k55);
01830 
01831          // Error estimate
01832 
01833          error1 = fabs(d1*k11+d2*k12+d3*k13+d4*k14+d5*k15+d6*k16);
01834          error2 = fabs(d1*k21+d2*k22+d3*k23+d4*k24+d5*k25+d6*k26);
01835          error3 = fabs(d1*k31+d2*k32+d3*k33+d4*k34+d5*k35+d6*k36);
01836          error4 = fabs(d1*k41+d2*k42+d3*k43+d4*k44+d5*k45+d6*k46);
01837 
01838          error = 0.0;
01839 
01840          if(error1 > error)
01841          {
01842             error = error1;
01843          }
01844          if(error2 > error)
01845          {
01846             error = error2;
01847          }
01848          if(error3 > error)
01849          {
01850             error = error3;
01851          }
01852          if(error4 > error)
01853          {
01854             error = error4;
01855          }
01856          if(error5 > error)
01857          {
01858             error = error5;
01859          }
01860 
01861          if(error <= tolerance)
01862          {
01863             x += h;
01864 
01865             y_1 += c51*k11+c52*k12+c53*k13+c54*k14+c55*k15+c56*k16;
01866             y_2 += c51*k21+c52*k22+c53*k23+c54*k24+c55*k25+c56*k26;
01867             y_3 += c51*k31+c52*k32+c53*k33+c54*k34+c55*k35+c56*k36;
01868             y_4 += c51*k41+c52*k42+c53*k43+c54*k44+c55*k45+c56*k46;
01869             y_5 += c51*k51+c52*k52+c53*k53+c54*k54+c55*k55+c56*k56;
01870 
01871             tempX[count] = x;
01872 
01873             tempY1[count] = y_1;
01874             tempY2[count] = y_2;
01875             tempY3[count] = y_3;
01876             tempY4[count] = y_4;
01877             tempY5[count] = y_5;
01878 
01879             count++;
01880 
01881             if(error != 0.0)
01882             {
01883                h *= 0.9*pow(fabs(tolerance/error),0.2);
01884             }
01885 
01886             if(count >= size)
01887             {
01888                // Not enought memory reserved 
01889                
01890                size *= 2;
01891 
01892             if(size > error_size)
01893             {
01894                   std::cerr << "Flood Error: OrdinaryDifferentialEquations class." << std::endl
01895                             << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01896                             << "Step size is bigger than greatest allowable." << std::endl;
01897                             
01898                   exit(1);
01899             }
01900             else if(display && (size > warning_size))
01901             {
01902                   std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
01903                             << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01904                             << "Size is " << size << std::endl;
01905             }
01906 
01907                tempX.resize(size);
01908                tempY1.resize(size);
01909                tempY2.resize(size);
01910                tempY3.resize(size);
01911                tempY4.resize(size);
01912                tempY5.resize(size);
01913          }
01914          }
01915          else
01916          {
01917             h *= 0.9*pow(fabs(tolerance/error),0.25);
01918          }
01919       } // end while loop   
01920 
01921       Vector<double> new_x(count);
01922 
01923       Vector<double> new_y1(count);
01924       Vector<double> new_y2(count);
01925       Vector<double> new_y3(count);
01926       Vector<double> new_y4(count);
01927       Vector<double> new_y5(count);
01928 
01929       for(int i = 0; i < count; i++)
01930       {
01931          new_x[i] = tempX[i];
01932          new_y1[i] = tempY1[i];
01933          new_y2[i] = tempY2[i];
01934          new_y3[i] = tempY3[i];
01935          new_y4[i] = tempY4[i];
01936          new_y5[i] = tempY5[i];
01937       }
01938 
01939       x_out = new_x;
01940       y1_out = new_y1;
01941       y2_out = new_y2;
01942       y3_out = new_y3;
01943       y4_out = new_y4;
01944       y5_out = new_y5;
01945 
01946      return(count);
01947    }
01948 
01949 private:
01950 
01951    // Runge-Kutta method parameters
01952 
01954 
01955    int points_number;
01956 
01957    // Runge-Kutta-Fehlberg method parameters
01958 
01961 
01962    double tolerance;
01963 
01965 
01966    int initial_size;
01967 
01969 
01970    int warning_size;
01971    
01973    
01974    int error_size;
01975 
01977 
01978    bool display;
01979 
01980 };
01981 
01982 }
01983 
01984 #endif
01985 
01986 
01987 // Flood: An Open Source Neural Networks C++ Library.
01988 // Copyright (C) 2005-2010 Roberto Lopez 
01989 //
01990 // This library is free software; you can redistribute it and/or
01991 // modify it under the s of the GNU Lesser General Public
01992 // License as published by the Free Software Foundation; either
01993 // version 2.1 of the License, or any later version.
01994 //
01995 // This library is distributed in the hope that it will be useful,
01996 // but WITHOUT ANY WARRANTY; without even the implied warranty of
01997 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
01998 // Lesser General Public License for more details.
01999 
02000 // You should have received a copy of the GNU Lesser General Public
02001 // License along with this library; if not, write to the Free Software
02002 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA