Flood: OrdinaryDifferentialEquations.cpp 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                                  */
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 #include "OrdinaryDifferentialEquations.h"
00017 
00018 #include<iostream>
00019 #include<fstream>
00020 #include<limits>
00021 #include<cmath>
00022 
00023 namespace Flood
00024 {
00025 
00026 // GENERAL CONSTRUCTOR
00027 
00029 
00030 OrdinaryDifferentialEquations::OrdinaryDifferentialEquations(void)
00031 {   
00032    set_default();
00033 }
00034   
00035 
00036 // DESTRUCTOR
00037 
00039 
00040 OrdinaryDifferentialEquations::~OrdinaryDifferentialEquations(void)
00041 {
00042 }
00043 
00044 
00045 // METHODS
00046 
00047 // GET METHODS
00048 
00049 // int get_points_number(void) method
00050 
00052 
00053 int OrdinaryDifferentialEquations::get_points_number(void)
00054 {
00055    return(points_number);
00056 }
00057 
00058 
00059 // double get_tolerance(void) method
00060 
00062 
00063 double OrdinaryDifferentialEquations::get_tolerance(void)
00064 {
00065    return(tolerance);
00066 }
00067 
00068 
00069 // int get_initial_size(void) method
00070 
00073 
00074 int OrdinaryDifferentialEquations::get_initial_size(void)
00075 {
00076    return(initial_size);
00077 }
00078 
00079 
00080 // int get_warning_size(void) method
00081 
00084 
00085 int OrdinaryDifferentialEquations::get_warning_size(void)
00086 {
00087    return(warning_size);
00088 }
00089 
00090 
00091 // int get_error_size(void) method
00092 
00094 
00095 int OrdinaryDifferentialEquations::get_error_size(void)
00096 {
00097    return(error_size);
00098 }
00099 
00100 
00101 // bool get_display(void) method
00102 
00105 
00106 bool OrdinaryDifferentialEquations::get_display(void)
00107 {
00108    return(display);
00109 }
00110 
00111 
00112 // SET METHODS
00113 
00114 // void set_default(void) method
00115 
00132 
00133 void OrdinaryDifferentialEquations::set_default(void)
00134 {
00135    // Runge-Kutta method parameters
00136 
00137    points_number = 1001;
00138 
00139    // Runge-Kutta-Fehlberg method parameters
00140 
00141    tolerance = 1.0e-9;
00142 
00143    initial_size = 1000;
00144    warning_size = 1000000;
00145    error_size = 1000000000;
00146 
00147    // Utilities
00148 
00149    display = true;
00150 }
00151 
00152 // void set_points_number(int) method
00153 
00156 
00157 void OrdinaryDifferentialEquations::set_points_number(int new_points_number)
00158 {
00159    points_number = new_points_number;
00160 }
00161 
00162 
00163 // void set_tolerance(double) method
00164 
00167 
00168 void OrdinaryDifferentialEquations::set_tolerance(double new_tolerance)
00169 {
00170    tolerance = new_tolerance;
00171 }
00172 
00173 
00174 // void set_initial_size(int) method
00175 
00179 
00180 void OrdinaryDifferentialEquations::set_initial_size(int new_initial_size)
00181 {
00182    initial_size = new_initial_size;
00183 }
00184 
00185 
00186 // void set_warning_size(int) method
00187 
00191 
00192 void OrdinaryDifferentialEquations::set_warning_size(int new_warning_size)
00193 {
00194    warning_size = new_warning_size;
00195 }
00196 
00197 
00198 // void set_error_size(int) method
00199 
00203 
00204 void OrdinaryDifferentialEquations::set_error_size(int new_error_size)
00205 {
00206    error_size = new_error_size;
00207 }
00208 
00209 
00210 // void set_display(bool) method
00211 
00216 
00217 void OrdinaryDifferentialEquations::set_display(bool new_display)
00218 {
00219    display = new_display;
00220 }
00221 
00222 
00223 // void OrdinaryDifferentialEquations::calculate_Runge_Kutta_integral(
00224 // Vector<double>& x, Vector<double>& y, 
00225 // double (*f)(double, double), 
00226 // double a, double b, double ya, int n)
00227 
00243 
00244 void OrdinaryDifferentialEquations::calculate_Runge_Kutta_integral(
00245 Vector<double>& x, 
00246 Vector<double>& y, 
00247 double (*f)(double, double), 
00248 double a, double b, double ya)
00249 {
00250    // Integration step
00251 
00252    double h = (b-a)/(points_number-1.0);      
00253 
00254    // Fourth order Runge-Kutta coefficients
00255 
00256    Vector<double> k(4);
00257 
00258    // Independent variable
00259 
00260    x[0] = a;
00261 
00262    for(int i = 0;  i < points_number-2; i++)
00263    {
00264       x[i+1] = x[i] + h;
00265    }
00266 
00267    x[points_number-1] = b;
00268 
00269    // Dependent variable
00270 
00271    y[0] = ya;
00272 
00273    for(int i = 0;  i < points_number-1; i++)
00274    {   
00275       // Runge-Kutta coefficients
00276 
00277       k[0] = f(x[i], y[i]);
00278       k[1] = f(x[i]+h/2.0, y[i]+h*k[0]/2.0);
00279       k[2] = f(x[i]+h/2.0, y[i]+h*k[1]/2.0);
00280       k[3] = f(x[i+1], y[i]+h*k[2]);
00281 
00282       // Dependent variable
00283 
00284       y[i+1] = y[i] + h*(k[0] + 2.0*k[1] + 2.0*k[2] + k[3])/6.0;
00285    }
00286 }
00287 
00288 
00289 // void calculate_Runge_Kutta_integral(
00290 // double (*f_1)(double, double, double), 
00291 // double (*f_2)(double, double,double),
00292 // Vector<double>&, Vector<double>&, Vector<double>&, 
00293 // double, double, double) method
00294 //
00315 
00316 void OrdinaryDifferentialEquations::calculate_Runge_Kutta_integral(
00317 Vector<double>& x, 
00318 Vector<double>& y_1, 
00319 Vector<double>& y_2, 
00320 double (*f_1)(double, double,double), 
00321 double (*f_2)(double, double,double), 
00322 double a, double b, 
00323 double y1a, double y2a)
00324 {
00325    // Integration step
00326 
00327    double h = (b-a)/(points_number-1.0);      
00328 
00329    // Fourth order Runge-Kutta coefficients
00330 
00331    Vector<double> k1(4);
00332    Vector<double> k2(4);
00333 
00334    // Independent variable
00335 
00336    x[0] = a;
00337 
00338    for(int i = 0;  i < points_number-2; i++)
00339    {
00340       x[i+1] = x[i] + h;
00341    }
00342 
00343    x[points_number-1] = b;
00344 
00345    // Dependent variables
00346   
00347    y_1[0] = y1a;
00348    y_2[0] = y2a;
00349 
00350    for(int i = 0;  i < points_number-1; i++)
00351    {
00352       // Runge-Kutta coefficients
00353 
00354       k1[0] = f_1(x[i], y_1[i], y_2[i]);
00355       k2[0] = f_2(x[i], y_1[i], y_2[i]);
00356 
00357       k1[1] = f_1(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0);
00358       k2[1] = f_2(x[i]+h/2.0, y_1[i]+h*k1[0]/2.0, y_2[i]+h*k2[0]/2.0);
00359 
00360       k1[2] = f_1(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0);
00361       k2[2] = f_2(x[i]+h/2.0, y_1[i]+h*k1[1]/2.0, y_2[i]+h*k2[1]/2.0);
00362 
00363       k1[3] = f_1(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2]);
00364       k2[3] = f_2(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2]);
00365 
00366       // Dependent variables
00367 
00368       y_1[i+1] = y_1[i] + h*(k1[0] + 2.0*k1[1] + 2.0*k1[2] + k1[3])/6.0;
00369       y_2[i+1] = y_2[i] + h*(k2[0] + 2.0*k2[1] + 2.0*k2[2] + k2[3])/6.0;
00370    }
00371 }
00372 
00373 
00374 // void calculate_Runge_Kutta_integral(
00375 // double (*f_1)(double, double, double, double), 
00376 // double (*f_2)(double, double, double, double),
00377 // double (*f_3)(double, double, double, double),   
00378 // Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&, 
00379 // double, double, double, double, double, int) method
00380 //
00407 
00408 void OrdinaryDifferentialEquations::calculate_Runge_Kutta_integral(
00409 Vector<double> &x, 
00410 Vector<double>& y_1, 
00411 Vector<double>& y_2, 
00412 Vector<double>& y_3, 
00413 double (*f_1)(double, double, double, double), 
00414 double (*f_2)(double, double, double, double), 
00415 double (*f_3)(double, double, double, double), 
00416 double a, double b, 
00417 double y1a, double y2a, double y3a)
00418 {
00419    // Integration step
00420 
00421    double h = (b-a)/(points_number-1.0);      
00422 
00423    // Fourth order Runge-Kutta coefficients
00424 
00425    Vector<double> k1(4);
00426    Vector<double> k2(4);
00427    Vector<double> k3(4);
00428 
00429    // Independent variable
00430 
00431    x[0] = a;
00432 
00433    for(int i = 0;  i < points_number-2; i++)
00434    {
00435       x[i+1] = x[i] + h;
00436    }
00437 
00438    x[points_number-1] = b;
00439 
00440    // Dependent variables
00441 
00442    y_1[0] = y1a;
00443    y_2[0] = y2a;
00444    y_3[0] = y3a;
00445 
00446    for(int i = 0;  i < points_number-1; i++)
00447    {
00448       // Runge-Kutta coefficients
00449 
00450       k1[0] = f_1(x[i], y_1[i], y_2[i], y_3[i]);
00451       k2[0] = f_2(x[i], y_1[i], y_2[i], y_3[i]);
00452       k3[0] = f_3(x[i], y_1[i], y_2[i], y_3[i]);
00453 
00454       k1[1] = 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);
00455       k2[1] = 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);
00456       k3[1] = 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);
00457 
00458       k1[2] = 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);
00459       k2[2] = 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);
00460       k3[2] = 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);
00461 
00462       k1[3] = f_1(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2]);
00463       k2[3] = f_2(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2]);
00464       k3[3] = f_2(x[i+1], y_1[i]+h*k1[2], y_2[i]+h*k2[2], y_3[i]+h*k3[2]);
00465 
00466       // Dependent variables
00467 
00468       y_1[i+1] = y_1[i] + h*(k1[0] + 2.0*k1[1] + 2.0*k1[2] + k1[3])/6.0;
00469       y_2[i+1] = y_2[i] + h*(k2[0] + 2.0*k2[1] + 2.0*k2[2] + k2[3])/6.0;
00470       y_3[i+1] = y_3[i] + h*(k3[0] + 2.0*k3[1] + 2.0*k3[2] + k3[3])/6.0;
00471    }
00472 }
00473 
00474 
00475 // void calculate_Runge_Kutta_integral(
00476 // double (*f_1)(double, double, double, double, double), 
00477 // double (*f_2)(double, double, double, double, double),
00478 // double (*f_3)(double, double, double, double, double),   
00479 // double (*f_4)(double, double, double, double, double),   
00480 // Vector<double>&, 
00481 // Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&, 
00482 // double, double, double, double, double, double, int) method
00483 
00513 
00514 void OrdinaryDifferentialEquations::calculate_Runge_Kutta_integral(
00515 Vector<double>& x, 
00516 Vector<double>& y_1, 
00517 Vector<double>& y_2, 
00518 Vector<double>& y_3, 
00519 Vector<double>& y_4, 
00520 double (*f_1)(double, double, double, double, double), 
00521 double (*f_2)(double, double, double, double, double), 
00522 double (*f_3)(double, double, double, double, double), 
00523 double (*f_4)(double, double, double, double, double), 
00524 double a, double b, 
00525 double y1a, double y2a, double y3a, double y4a)
00526 {
00527    // Integration step
00528 
00529    double h = (b-a)/(points_number-1.0);
00530 
00531    // Fourth order Runge-Kutta coefficients
00532 
00533    Vector<double> k1(4);
00534    Vector<double> k2(4);
00535    Vector<double> k3(4);
00536    Vector<double> k4(4);
00537 
00538    // Independent variable
00539 
00540    x[0] = a;
00541 
00542    for(int i = 0;  i < points_number-2; i++)
00543    {
00544       x[i+1] = x[i] + h;
00545    }
00546 
00547    x[points_number-1] = b;
00548 
00549    // Dependent variables
00550 
00551    y_1[0] = y1a;
00552    y_2[0] = y2a;
00553    y_3[0] = y3a;
00554    y_4[0] = y4a;
00555 
00556    for(int i = 0;  i < points_number-1; i++)
00557    {
00558       // Runge-Kutta coefficients
00559 
00560       k1[0] = f_1(x[i], y_1[i], y_2[i], y_3[i], y_4[i]);
00561       k2[0] = f_2(x[i], y_1[i], y_2[i], y_3[i], y_4[i]);
00562       k3[0] = f_3(x[i], y_1[i], y_2[i], y_3[i], y_4[i]);
00563       k4[0] = f_4(x[i], y_1[i], y_2[i], y_3[i], y_4[i]);
00564 
00565       k1[1] = f_1(x[i]+h/2.0, 
00566       y_1[i]+h*k1[0]/2.0, 
00567       y_2[i]+h*k2[0]/2.0, 
00568       y_3[i]+h*k3[0]/2.0, 
00569       y_4[i]+h*k4[0]/2.0);
00570 
00571       k2[1] = f_2(x[i]+h/2.0, 
00572       y_1[i]+h*k1[0]/2.0, 
00573       y_2[i]+h*k2[0]/2.0, 
00574       y_3[i]+h*k3[0]/2.0, 
00575       y_4[i]+h*k4[0]/2.0);
00576 
00577       k3[1] = f_3(x[i]+h/2.0, 
00578       y_1[i]+h*k1[0]/2.0, 
00579       y_2[i]+h*k2[0]/2.0, 
00580       y_3[i]+h*k3[0]/2.0, 
00581       y_4[i]+h*k4[0]/2.0);
00582 
00583       k4[1] = f_4(x[i]+h/2.0, 
00584       y_1[i]+h*k1[0]/2.0, 
00585       y_2[i]+h*k2[0]/2.0, 
00586       y_3[i]+h*k3[0]/2.0, 
00587       y_4[i]+h*k4[0]/2.0);
00588 
00589       k1[2] = f_1(x[i]+h/2.0, 
00590       y_1[i]+h*k1[1]/2.0, 
00591       y_2[i]+h*k2[1]/2.0, 
00592       y_3[i]+h*k3[1]/2.0, 
00593       y_4[i]+h*k4[1]/2.0);
00594 
00595       k2[2] = f_2(x[i]+h/2.0, 
00596       y_1[i]+h*k1[1]/2.0, 
00597       y_2[i]+h*k2[1]/2.0, 
00598       y_3[i]+h*k3[1]/2.0, 
00599       y_4[i]+h*k4[1]/2.0);
00600 
00601       k3[2] = f_3(x[i]+h/2.0, 
00602       y_1[i]+h*k1[1]/2.0, 
00603       y_2[i]+h*k2[1]/2.0, 
00604       y_3[i]+h*k3[1]/2.0, 
00605       y_4[i]+h*k4[1]/2.0);
00606 
00607       k4[2] = f_4(x[i]+h/2.0, 
00608       y_1[i]+h*k1[1]/2.0, 
00609       y_2[i]+h*k2[1]/2.0, 
00610       y_3[i]+h*k3[1]/2.0, 
00611       y_4[i]+h*k4[1]/2.0);
00612 
00613       k1[3] = 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]);
00614       k2[3] = 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]);
00615       k3[3] = 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]);
00616       k4[3] = 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]);
00617 
00618       // Dependent variables
00619 
00620       y_1[i+1] = y_1[i] + h*(k1[0] + 2.0*k1[1] + 2.0*k1[2] + k1[3])/6.0;
00621       y_2[i+1] = y_2[i] + h*(k2[0] + 2.0*k2[1] + 2.0*k2[2] + k2[3])/6.0;
00622       y_3[i+1] = y_3[i] + h*(k3[0] + 2.0*k3[1] + 2.0*k3[2] + k3[3])/6.0;
00623       y_4[i+1] = y_4[i] + h*(k4[0] + 2.0*k4[1] + 2.0*k4[2] + k4[3])/6.0;
00624    }
00625 }
00626 
00627 
00628 // void calculate_Runge_Kutta_integral(
00629 // double (*f_1)(double, double, double, double, double), 
00630 // double (*f_2)(double, double, double, double, double),
00631 // double (*f_3)(double, double, double, double, double),   
00632 // double (*f_4)(double, double, double, double, double),   
00633 // Vector<double>&, 
00634 // Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&, 
00635 // double, double, double, double, double, double) method
00636 
00670 
00671 void OrdinaryDifferentialEquations::calculate_Runge_Kutta_integral(
00672 Vector<double>& x, 
00673 Vector<double>& y_1, 
00674 Vector<double>& y_2, 
00675 Vector<double>& y_3, 
00676 Vector<double>& y_4, 
00677 Vector<double>& y_5, 
00678 double (*f_1)(double, double, double, double, double, double), 
00679 double (*f_2)(double, double, double, double, double, double), 
00680 double (*f_3)(double, double, double, double, double, double), 
00681 double (*f_4)(double, double, double, double, double, double), 
00682 double (*f_5)(double, double, double, double, double, double), 
00683 double a, double b, 
00684 double y1a, double y2a, double y3a, double y4a, double y5a)
00685 {
00686    // Integration step
00687 
00688    double h = (b-a)/(points_number-1.0);
00689 
00690    // Fourth order Runge-Kutta coefficients
00691 
00692    Vector<double> k1(4);
00693    Vector<double> k2(4);
00694    Vector<double> k3(4);
00695    Vector<double> k4(4);
00696    Vector<double> k5(4);
00697 
00698    // Independent variable
00699 
00700    x[0] = a;
00701 
00702    for(int i = 0;  i < points_number-2; i++)
00703    {
00704       x[i+1] = x[i] + h;
00705    }
00706 
00707    x[points_number-1] = b;
00708 
00709    // Dependent variables
00710 
00711    y_1[0] = y1a;
00712    y_2[0] = y2a;
00713    y_3[0] = y3a;
00714    y_4[0] = y4a;
00715    y_5[0] = y5a;
00716 
00717    for(int i = 0;  i < points_number-1; i++)
00718    {
00719       // Runge-Kutta coefficients
00720 
00721       k1[0] = f_1(x[i], y_1[i], y_2[i], y_3[i], y_4[i], y_5[i]);
00722       k2[0] = f_2(x[i], y_1[i], y_2[i], y_3[i], y_4[i], y_5[i]);
00723       k3[0] = f_3(x[i], y_1[i], y_2[i], y_3[i], y_4[i], y_5[i]);
00724       k4[0] = f_4(x[i], y_1[i], y_2[i], y_3[i], y_4[i], y_5[i]);
00725       k5[0] = f_5(x[i], y_1[i], y_2[i], y_3[i], y_4[i], y_5[i]);
00726 
00727       k1[1] = 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);
00728       k2[1] = 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);
00729       k3[1] = 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);
00730       k4[1] = 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);
00731       k5[1] = 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);
00732 
00733       k1[2] = 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);
00734       k2[2] = 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);
00735       k3[2] = 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);
00736       k4[2] = 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);
00737       k5[2] = 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);
00738 
00739       k1[3] = 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]);
00740       k2[3] = 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]);
00741       k3[3] = 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]);
00742       k4[3] = 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]);
00743       k5[3] = 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]);
00744 
00745       // Dependent variables
00746 
00747       y_1[i+1] = y_1[i] + h*(k1[0] + 2.0*k1[1] + 2.0*k1[2] + k1[3])/6.0;
00748       y_2[i+1] = y_2[i] + h*(k2[0] + 2.0*k2[1] + 2.0*k2[2] + k2[3])/6.0;
00749       y_3[i+1] = y_3[i] + h*(k3[0] + 2.0*k3[1] + 2.0*k3[2] + k3[3])/6.0;
00750       y_4[i+1] = y_4[i] + h*(k4[0] + 2.0*k4[1] + 2.0*k4[2] + k4[3])/6.0;
00751       y_5[i+1] = y_5[i] + h*(k5[0] + 2.0*k5[1] + 2.0*k5[2] + k5[3])/6.0;
00752    }
00753 }
00754 
00755 
00756 // int calculate_Runge_Kutta_Fehlberg_integral(Vector<double>& x, Vector<double>& y, 
00757 // double (*f)(double, double), 
00758 // double a, double b, double ya) method
00759 
00776 
00777 int OrdinaryDifferentialEquations
00778 ::calculate_Runge_Kutta_Fehlberg_integral(Vector<double>& x_out, Vector<double>& yOut, 
00779 double (*f)(double, double), double a, double b, double ya)
00780 {
00781    double error = 0.0;
00782 
00783    int count = 0;
00784    int size = initial_size;
00785 
00786    Vector<double> tempX(size);
00787    Vector<double> tempY(size);
00788 
00789    tempX[0] = a;
00790    tempY[0] = ya;
00791 
00792    count++;
00793 
00794    double epsilon = std::numeric_limits<double>::epsilon();
00795 
00796    double a2 = 1.0/5.0;
00797    double a3 = 3.0/10.0;
00798    double a4 = 3.0/5.0;
00799    double a5 = 1.0;
00800    double a6 = 7.0/8.0;
00801 
00802    double b21 = 1.0/5.0;   
00803    double b31 = 3.0/40.0;
00804    double b32 = 9.0/40.0;   
00805    double b41 = 3.0/10.0;
00806    double b42 = -9.0/10.0;   
00807    double b43 = 6.0/5.0; 
00808    double b51 = -11.0/54.0;
00809    double b52 = 5.0/2.0;
00810    double b53 = -70.0/27.0;
00811    double b54 = 35.0/27.0; 
00812    double b61 = 1631.0/55296.0;
00813    double b62 = 175.0/512.0;   
00814    double b63 = 575.0/13824.0;   
00815    double b64 = 44275.0/110592.0;   
00816    double b65 = 253.0/4096.0;
00817    
00818    double c41 = 37.0/378.0;
00819    double c42 = 0.0;
00820    double c43 = 250.0/621.0;
00821    double c44 = 125.0/594.0;
00822    double c45 = 0.0;
00823    double c46 = 512.0/1771.0;
00824 
00825    double c51 = 2825.0/27648.0;  
00826    double c52 = 0.0;  
00827    double c53 = 18575.0/48384.0;   
00828    double c54 = 13525.0/55296.0;
00829    double c55 = 277.0/14336.0;
00830    double c56 = 1.0/4.0;
00831 
00832    double d1 = c41 - c51;
00833    double d2 = c42 - c52;
00834    double d3 = c43 - c53;
00835    double d4 = c44 - c54;
00836    double d5 = c45 - c55;
00837    double d6 = c46 - c56;
00838 
00839    double k1 = 0.0;
00840    double k2 = 0.0;
00841    double k3 = 0.0;
00842    double k4 = 0.0;
00843    double k5 = 0.0;
00844    double k6 = 0.0;
00845 
00846    // Main loop
00847 
00848    double x = a;
00849    double y = ya;
00850 
00851    double hmin = 16.0*epsilon*fabs(x);
00852 
00853    double h = (b-a)*1.0e-3;
00854 
00855    while(x < b)
00856    { 
00857       // Set smallest allowable stepsize
00858 
00859       hmin = 32.0*epsilon*fabs(x);
00860 
00861       if(h < hmin)
00862       {
00863          std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
00864                    << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
00865                    << "Step size is less than smallest allowable." << std::endl;
00866 
00867          h = hmin;
00868       }
00869 
00870       if(x + h > b)
00871       {
00872          h = b - x;
00873       }
00874 
00875       // Runge-Kutta-Fehlberg coefficients
00876 
00877       k1 = h*f(x, y);
00878       k2 = h*f(x+a2*h, y+b21*k1);
00879       k3 = h*f(x+a3*h, y+b31*k1+b32*k2);
00880       k4 = h*f(x+a4*h, y+b41*k1+b42*k2+b43*k3); 
00881       k5 = h*f(x+a5*h, y+b51*k1+b52*k2+b53*k3+b54*k4); 
00882       k6 = h*f(x+a6*h, y+b61*k1+b62*k2+b63*k3+b64*k4+b65*k5);
00883 
00884       // Error estimate
00885 
00886       error = fabs(d1*k1+d2*k2+d3*k3+d4*k4+d5*k5+d6*k6);
00887 
00888       if(error <= tolerance)
00889       {
00890          x += h;
00891          y += c51*k1+c52*k2+c53*k3+c54*k4+c55*k5+c56*k6;
00892 
00893          tempX[count] = x;
00894          tempY[count] = y;
00895 
00896          count++;
00897 
00898          if(error != 0.0)
00899          {
00900             h *= 0.9*pow(fabs(tolerance/error),0.2);      
00901          }
00902 
00903          if(count >= size)
00904          {
00905             std::cerr << "Flood Error: calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
00906                       << "Insufficient memory reserved." << std::endl;
00907 
00908             exit(1);
00909          }
00910       }
00911       else
00912       {
00913          h *= 0.9*pow(fabs(tolerance/error),0.25);
00914       }
00915    }// end while loop   
00916 
00917    Vector<double> new_x(count);
00918    Vector<double> newY(count);
00919 
00920    for(int i = 0; i < count; i++)
00921    {
00922       new_x[i] = tempX[i];
00923       newY[i] = tempY[i];
00924    }
00925 
00926    x_out = new_x;
00927    yOut = newY;
00928 
00929    return(count);
00930 }
00931 
00932 
00933 // int calculate_Runge_Kutta_Fehlberg_integral(Vector<double>&, Vector<double>&, Vector<double>&, 
00934 // double (*f_1)(double, double, double), 
00935 // double (*f_2)(double, double, double), 
00936 // double, double, double, double)
00937 
00959 
00960 int OrdinaryDifferentialEquations::calculate_Runge_Kutta_Fehlberg_integral(
00961 Vector<double>& x_out, Vector<double>& y1_out, Vector<double>& y2_out, 
00962 double (*f_1)(double, double, double), 
00963 double (*f_2)(double, double, double), 
00964 double a, double b, double y1a, double y2a)
00965 {
00966    double error = 0.0, error1 = 0.0, error2 = 0.0;
00967 
00968    int count = 0;
00969    int size = initial_size;
00970 
00971    Vector<double> tempX(size);
00972    Vector<double> tempY1(size);
00973    Vector<double> tempY2(size);
00974 
00975    tempX[0] = a;
00976 
00977    tempY1[0] = y1a;
00978    tempY2[0] = y2a;
00979 
00980    count++;
00981 
00982    double epsilon = std::numeric_limits<double>::epsilon();
00983 
00984    double a2 = 1.0/5.0;
00985    double a3 = 3.0/10.0;
00986    double a4 = 3.0/5.0;
00987    double a5 = 1.0;
00988    double a6 = 7.0/8.0;
00989 
00990    double b21 = 1.0/5.0;   
00991    double b31 = 3.0/40.0;
00992    double b32 = 9.0/40.0;   
00993    double b41 = 3.0/10.0;
00994    double b42 = -9.0/10.0;   
00995    double b43 = 6.0/5.0; 
00996    double b51 = -11.0/54.0;
00997    double b52 = 5.0/2.0;
00998    double b53 = -70.0/27.0;
00999    double b54 = 35.0/27.0; 
01000    double b61 = 1631.0/55296.0;
01001    double b62 = 175.0/512.0;   
01002    double b63 = 575.0/13824.0;   
01003    double b64 = 44275.0/110592.0;   
01004    double b65 = 253.0/4096.0;
01005 
01006    double c41 = 37.0/378.0;
01007    double c42 = 0.0;
01008    double c43 = 250.0/621.0;
01009    double c44 = 125.0/594.0;
01010    double c45 = 0.0;
01011    double c46 = 512.0/1771.0;
01012 
01013    double c51 = 2825.0/27648.0;  
01014    double c52 = 0.0;  
01015    double c53 = 18575.0/48384.0;   
01016    double c54 = 13525.0/55296.0;
01017    double c55 = 277.0/14336.0;
01018    double c56 = 1.0/4.0;
01019 
01020    double d1 = c41 - c51;
01021    double d2 = c42 - c52;
01022    double d3 = c43 - c53;
01023    double d4 = c44 - c54;
01024    double d5 = c45 - c55;
01025    double d6 = c46 - c56;
01026 
01027    double k11 = 0.0, k21 = 0.0;
01028    double k12 = 0.0, k22 = 0.0;
01029    double k13 = 0.0, k23 = 0.0;
01030    double k14 = 0.0, k24 = 0.0;
01031    double k15 = 0.0, k25 = 0.0;
01032    double k16 = 0.0, k26 = 0.0;
01033 
01034    // Main loop
01035 
01036    double x = a;
01037    double y_1 = y1a;
01038    double y_2 = y2a;
01039 
01040    double hmin = 16.0*epsilon*fabs(x);
01041 
01042    double h = (b-a)*1.0e-3;
01043 
01044    while(x < b)
01045    {
01046       // Set smallest allowable stepsize
01047 
01048       hmin = 32.0*epsilon*fabs(x);
01049 
01050       if(h < hmin)
01051       {
01052          std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
01053                    << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01054                    << "Step size is less than smallest allowable." << std::endl;
01055 
01056        h = hmin;
01057       }
01058 
01059       if(x + h > b)
01060       {
01061          h = b - x;
01062       }
01063 
01064       // Runge-Kutta-Fehlberg coefficients
01065 
01066       k11 = h*f_1(x, y_1, y_2);
01067       k21 = h*f_2(x, y_1, y_2);
01068 
01069       k12 = h*f_1(x+a2*h, y_1+b21*k11, y_2+b21*k21);
01070       k22 = h*f_2(x+a2*h, y_1+b21*k11, y_2+b21*k21);
01071 
01072       k13 = h*f_1(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22);
01073       k23 = h*f_2(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22);
01074 
01075       k14 = h*f_1(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23); 
01076       k24 = h*f_2(x+a4*h, y_1+b41*k11+b42*k12+b43*k13, y_2+b41*k21+b42*k22+b43*k23); 
01077 
01078       k15 = h*f_1(x+a5*h, 
01079       y_1+b51*k11+b52*k12+b53*k13+b54*k14, 
01080       y_2+b51*k21+b52*k22+b53*k23+b54*k24); 
01081 
01082       k25 = h*f_2(x+a5*h, 
01083       y_1+b51*k11+b52*k12+b53*k13+b54*k14, 
01084       y_2+b51*k21+b52*k22+b53*k23+b54*k24); 
01085 
01086       k16 = h*f_1(x+a6*h, 
01087       y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, 
01088       y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25);
01089 
01090       k26 = h*f_2(x+a6*h, 
01091       y_1+b61*k11+b62*k12+b63*k13+b64*k14+b65*k15, 
01092       y_2+b61*k21+b62*k22+b63*k23+b64*k24+b65*k25);
01093 
01094       // Error estimate
01095 
01096       error1 = fabs(d1*k11+d2*k12+d3*k13+d4*k14+d5*k15+d6*k16);
01097       error2 = fabs(d1*k21+d2*k22+d3*k23+d4*k24+d5*k25+d6*k26);
01098 
01099       error = 0.0;
01100       if(error1 > error)
01101       {
01102          error = error1;
01103       }
01104       if(error2 > error)
01105       {
01106          error = error2;
01107       }
01108 
01109       if(error <= tolerance)
01110       {
01111          x += h;
01112          y_1 += c51*k11+c52*k12+c53*k13+c54*k14+c55*k15+c56*k16;
01113          y_2 += c51*k21+c52*k22+c53*k23+c54*k24+c55*k25+c56*k26;
01114 
01115          tempX[count] = x;
01116          tempY1[count] = y_1;
01117          tempY2[count] = y_2;
01118          count++;
01119 
01120          if(error != 0.0)
01121          {
01122             h *= 0.9*pow(fabs(tolerance/error),0.2);
01123          }
01124 
01125          if(count >= size)
01126          {
01127             std::cerr << "Flood Error: calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01128                       << "Insufficient memory reserved." << std::endl;
01129 
01130             exit(1);
01131          }
01132       }
01133       else
01134       {
01135          h *= 0.9*pow(fabs(tolerance/error),0.25);
01136       }
01137    } // end while loop   
01138 
01139    Vector<double> new_x(count);
01140    Vector<double> new_y1(count);
01141    Vector<double> new_y2(count);
01142 
01143    for(int i = 0; i < count; i++)
01144    {
01145       new_x[i] = tempX[i];
01146       new_y1[i] = tempY1[i];
01147       new_y2[i] = tempY2[i];
01148    }
01149 
01150    x_out = new_x;
01151    y1_out = new_y1;
01152    y2_out = new_y2;
01153 
01154    return(count);
01155 }
01156 
01157 
01158 // int calculate_Runge_Kutta_Fehlberg_integral(Vector<double>&,
01159 // Vector<double>&, Vector<double>&, Vector<double>&, 
01160 // double (*f_1)(double, double, double, double), 
01161 // double (*f_2)(double, double, double, double), 
01162 // double (*f_3)(double, double, double, double), 
01163 // double, double, double, double, double) method
01164 
01190 
01191 int OrdinaryDifferentialEquations::calculate_Runge_Kutta_Fehlberg_integral(
01192 Vector<double>& x_out, Vector<double>& y1_out, Vector<double>& y2_out, Vector<double>& y3_out,
01193 double (*f_1)(double, double, double, double), 
01194 double (*f_2)(double, double, double, double), 
01195 double (*f_3)(double, double, double, double), 
01196 double a, double b, double y1a, double y2a, double y3a)
01197 {
01198    double error = 0.0, error1 = 0.0, error2 = 0.0, error3 = 0.0;
01199 
01200    int size = initial_size;
01201 
01202    int count = 0;
01203 
01204    Vector<double> tempX(size);
01205    Vector<double> tempY1(size);
01206    Vector<double> tempY2(size);
01207    Vector<double> tempY3(size);
01208 
01209    tempX[0] = a;
01210 
01211    tempY1[0] = y1a;
01212    tempY2[0] = y2a;
01213    tempY3[0] = y3a;
01214 
01215    count++;
01216 
01217    double epsilon = std::numeric_limits<double>::epsilon();
01218 
01219    double a2 = 1.0/5.0;
01220    double a3 = 3.0/10.0;
01221    double a4 = 3.0/5.0;
01222    double a5 = 1.0;
01223    double a6 = 7.0/8.0;
01224 
01225    double b21 = 1.0/5.0;   
01226    double b31 = 3.0/40.0;
01227    double b32 = 9.0/40.0;   
01228    double b41 = 3.0/10.0;
01229    double b42 = -9.0/10.0;   
01230    double b43 = 6.0/5.0; 
01231    double b51 = -11.0/54.0;
01232    double b52 = 5.0/2.0;
01233    double b53 = -70.0/27.0;
01234    double b54 = 35.0/27.0; 
01235    double b61 = 1631.0/55296.0;
01236    double b62 = 175.0/512.0;   
01237    double b63 = 575.0/13824.0;   
01238    double b64 = 44275.0/110592.0;   
01239    double b65 = 253.0/4096.0;
01240    
01241    double c41 = 37.0/378.0;
01242    double c42 = 0.0;
01243    double c43 = 250.0/621.0;
01244    double c44 = 125.0/594.0;
01245    double c45 = 0.0;
01246    double c46 = 512.0/1771.0;
01247 
01248    double c51 = 2825.0/27648.0;  
01249    double c52 = 0.0;  
01250    double c53 = 18575.0/48384.0;   
01251    double c54 = 13525.0/55296.0;
01252    double c55 = 277.0/14336.0;
01253    double c56 = 1.0/4.0;
01254 
01255    double d1 = c41 - c51;
01256    double d2 = c42 - c52;
01257    double d3 = c43 - c53;
01258    double d4 = c44 - c54;
01259    double d5 = c45 - c55;
01260    double d6 = c46 - c56;
01261 
01262    double k11 = 0.0, k21 = 0.0, k31 = 0.0;
01263    double k12 = 0.0, k22 = 0.0, k32 = 0.0;
01264    double k13 = 0.0, k23 = 0.0, k33 = 0.0;
01265    double k14 = 0.0, k24 = 0.0, k34 = 0.0;
01266    double k15 = 0.0, k25 = 0.0, k35 = 0.0;
01267    double k16 = 0.0, k26 = 0.0, k36 = 0.0;
01268 
01269    // Main loop
01270 
01271    double x = a;
01272    double y_1 = y1a;
01273    double y_2 = y2a;
01274    double y_3 = y3a;
01275 
01276    double hmin = 16.0*epsilon*fabs(x);
01277 
01278    double h = (b-a)*1.0e-3;
01279 
01280    while(x < b)
01281    {   
01282       // Set smallest allowable stepsize
01283 
01284       hmin = 32.0*epsilon*fabs(x);
01285 
01286       if(h < hmin)
01287       {
01288          std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
01289                    << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01290                    << "Step size is less than smallest allowable." << std::endl;
01291 
01292        h = hmin;
01293       }
01294 
01295       if(x + h > b)
01296       {
01297          h = b - x;
01298       }
01299 
01300       // Runge-Kutta-Fehlberg coefficients
01301 
01302       k11 = h*f_1(x, y_1, y_2, y_3);
01303       k21 = h*f_2(x, y_1, y_2, y_3);
01304       k31 = h*f_3(x, y_1, y_2, y_3);
01305 
01306       k12 = h*f_1(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31);
01307       k22 = h*f_2(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31);
01308       k32 = h*f_3(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31);
01309 
01310       k13 = h*f_1(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32);
01311       k23 = h*f_2(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32);
01312       k33 = h*f_3(x+a3*h, y_1+b31*k11+b32*k12, y_2+b31*k21+b32*k22, y_3+b31*k31+b32*k32);
01313 
01314       k14 = h*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); 
01315       k24 = h*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); 
01316       k34 = h*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); 
01317 
01318       k15 = h*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); 
01319       k25 = h*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); 
01320       k35 = h*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); 
01321 
01322       k16 = h*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);
01323       k26 = h*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);
01324       k36 = h*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);
01325 
01326       // Error estimate
01327 
01328       error1 = fabs(d1*k11+d2*k12+d3*k13+d4*k14+d5*k15+d6*k16);
01329       error2 = fabs(d1*k21+d2*k22+d3*k23+d4*k24+d5*k25+d6*k26);
01330       error3 = fabs(d1*k31+d2*k32+d3*k33+d4*k34+d5*k35+d6*k36);
01331 
01332       error = 0.0;
01333       if(error1 > error)
01334       {
01335          error = error1;
01336       }
01337       if(error2 > error)
01338       {
01339          error = error2;
01340       }
01341       if(error3 > error)
01342       {
01343          error = error3;
01344       }
01345       if(error <= tolerance)
01346       {
01347          x += h;
01348 
01349          y_1 += c51*k11+c52*k12+c53*k13+c54*k14+c55*k15+c56*k16;
01350          y_2 += c51*k21+c52*k22+c53*k23+c54*k24+c55*k25+c56*k26;
01351          y_3 += c51*k31+c52*k32+c53*k33+c54*k34+c55*k35+c56*k36;
01352 
01353          tempX[count] = x;
01354          tempY1[count] = y_1;
01355          tempY2[count] = y_2;
01356          tempY3[count] = y_3;
01357 
01358          count++;
01359 
01360          if(error != 0.0)
01361          {
01362             h *= 0.9*pow(fabs(tolerance/error),0.2);
01363          }
01364 
01365          if(count >= size)
01366          {
01367             std::cout << "Flood Error: calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01368                       << "Insufficient memory reserved." << std::endl;
01369 
01370             exit(1);
01371          }
01372       }
01373       else
01374       {
01375          h *= 0.9*pow(fabs(tolerance/error),0.25);      
01376       }
01377    } // end while loop   
01378 
01379    Vector<double> new_x(count);
01380    Vector<double> new_y1(count);
01381    Vector<double> new_y2(count);
01382    Vector<double> new_y3(count);
01383 
01384    for(int i = 0; i < count; i++)
01385    {
01386       new_x[i] = tempX[i];
01387       new_y1[i] = tempY1[i];
01388       new_y2[i] = tempY2[i];
01389       new_y3[i] = tempY3[i];
01390    }
01391 
01392    x_out = new_x;
01393    y1_out = new_y1;
01394    y2_out = new_y2;
01395    y3_out = new_y3;
01396 
01397    return(count);
01398 }
01399 
01400 
01401 // int calculate_Runge_Kutta_Fehlberg_integral(Vector<double>&,
01402 // Vector<double>&, Vector<double>&, Vector<double>&, Vector<double>&,
01403 // double (*f_1)(double, double, double, double, double), 
01404 // double (*f_2)(double, double, double, double, double), 
01405 // double (*f_3)(double, double, double, double, double), 
01406 // double (*f_4)(double, double, double, double, double), 
01407 // double, double, double, double, double, double) method
01408 
01438 
01439 int OrdinaryDifferentialEquations::calculate_Runge_Kutta_Fehlberg_integral(
01440 Vector<double>& x_out, Vector<double>& y1_out, Vector<double>& y2_out, Vector<double>& y3_out, Vector<double>& y4_out,
01441 double (*f_1)(double, double, double, double, double), 
01442 double (*f_2)(double, double, double, double, double), 
01443 double (*f_3)(double, double, double, double, double), 
01444 double (*f_4)(double, double, double, double, double), 
01445 double a, double b, double y1a, double y2a, double y3a, double y4a)
01446 {
01447    double error = 0.0, error1 = 0.0, error2 = 0.0, error3 = 0.0, error4 = 0.0;
01448 
01449    int size = initial_size;
01450 
01451    int count = 0;
01452 
01453    Vector<double> tempX(size);
01454    Vector<double> tempY1(size);
01455    Vector<double> tempY2(size);
01456    Vector<double> tempY3(size);
01457    Vector<double> tempY4(size);
01458 
01459    tempX[0] = a;
01460 
01461    tempY1[0] = y1a;
01462    tempY2[0] = y2a;
01463    tempY3[0] = y3a;
01464    tempY4[0] = y4a;
01465 
01466    count++;
01467 
01468    double epsilon = std::numeric_limits<double>::epsilon();
01469 
01470    double a2 = 1.0/5.0;
01471    double a3 = 3.0/10.0;
01472    double a4 = 3.0/5.0;
01473    double a5 = 1.0;
01474    double a6 = 7.0/8.0;
01475 
01476    double b21 = 1.0/5.0;   
01477    double b31 = 3.0/40.0;
01478    double b32 = 9.0/40.0;   
01479    double b41 = 3.0/10.0;
01480    double b42 = -9.0/10.0;   
01481    double b43 = 6.0/5.0; 
01482    double b51 = -11.0/54.0;
01483    double b52 = 5.0/2.0;
01484    double b53 = -70.0/27.0;
01485    double b54 = 35.0/27.0; 
01486    double b61 = 1631.0/55296.0;
01487    double b62 = 175.0/512.0;   
01488    double b63 = 575.0/13824.0;   
01489    double b64 = 44275.0/110592.0;   
01490    double b65 = 253.0/4096.0;
01491    
01492    double c41 = 37.0/378.0;
01493    double c42 = 0.0;
01494    double c43 = 250.0/621.0;
01495    double c44 = 125.0/594.0;
01496    double c45 = 0.0;
01497    double c46 = 512.0/1771.0;
01498 
01499    double c51 = 2825.0/27648.0;  
01500    double c52 = 0.0;  
01501    double c53 = 18575.0/48384.0;   
01502    double c54 = 13525.0/55296.0;
01503    double c55 = 277.0/14336.0;
01504    double c56 = 1.0/4.0;
01505 
01506    double d1 = c41 - c51;
01507    double d2 = c42 - c52;
01508    double d3 = c43 - c53;
01509    double d4 = c44 - c54;
01510    double d5 = c45 - c55;
01511    double d6 = c46 - c56;
01512 
01513    double k11 = 0.0, k21 = 0.0, k31 = 0.0, k41 = 0.0;
01514    double k12 = 0.0, k22 = 0.0, k32 = 0.0, k42 = 0.0;
01515    double k13 = 0.0, k23 = 0.0, k33 = 0.0, k43 = 0.0;
01516    double k14 = 0.0, k24 = 0.0, k34 = 0.0, k44 = 0.0;
01517    double k15 = 0.0, k25 = 0.0, k35 = 0.0, k45 = 0.0;
01518    double k16 = 0.0, k26 = 0.0, k36 = 0.0, k46 = 0.0;
01519 
01520    // Main loop
01521 
01522    double x = a;
01523    double y_1 = y1a;
01524    double y_2 = y2a;
01525    double y_3 = y3a;
01526    double y_4 = y4a;
01527 
01528    double hmin = 16.0*epsilon*fabs(x);
01529 
01530    double h = (b-a)*1.0e-3;
01531 
01532    while(x < b)
01533    {   
01534       // Set smallest allowable stepsize
01535 
01536       hmin = 32.0*epsilon*fabs(x);
01537 
01538       if(h < hmin)
01539       {
01540          std::cout << "Flood Warning: OrdinaryDifferentialEquations class." << std::endl
01541                    << "calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01542                    << "Step size is less than smallest allowable." << std::endl;
01543 
01544        h = hmin;
01545       }
01546 
01547       if(x + h > b)
01548       {
01549          h = b - x;
01550       }
01551 
01552       // Runge-Kutta-Fehlberg coefficients
01553 
01554       k11 = h*f_1(x, y_1, y_2, y_3, y_4);
01555       k21 = h*f_2(x, y_1, y_2, y_3, y_4);
01556       k31 = h*f_3(x, y_1, y_2, y_3, y_4);
01557       k41 = h*f_4(x, y_1, y_2, y_3, y_4);
01558 
01559       k12 = h*f_1(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41);
01560       k22 = h*f_2(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41);
01561       k32 = h*f_3(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41);
01562       k42 = h*f_4(x+a2*h, y_1+b21*k11, y_2+b21*k21, y_3+b21*k31, y_4+b21*k41);
01563 
01564       k13 = h*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);
01565       k23 = h*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);
01566       k33 = h*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);
01567       k43 = h*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);
01568 
01569       k14 = h*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); 
01570       k24 = h*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); 
01571       k34 = h*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); 
01572       k44 = h*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); 
01573 
01574       k15 = h*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); 
01575       k25 = h*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); 
01576       k35 = h*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); 
01577       k45 = h*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); 
01578 
01579       k16 = h*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);
01580       k26 = h*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);
01581       k36 = h*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);
01582       k46 = h*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);
01583 
01584       // Error estimate
01585 
01586       error1 = fabs(d1*k11+d2*k12+d3*k13+d4*k14+d5*k15+d6*k16);
01587       error2 = fabs(d1*k21+d2*k22+d3*k23+d4*k24+d5*k25+d6*k26);
01588       error3 = fabs(d1*k31+d2*k32+d3*k33+d4*k34+d5*k35+d6*k36);
01589       error4 = fabs(d1*k41+d2*k42+d3*k43+d4*k44+d5*k45+d6*k46);
01590 
01591       error = 0.0;
01592       if(error1 > error)
01593       {
01594          error = error1;
01595       }
01596       if(error2 > error)
01597       {
01598          error = error2;
01599       }
01600       if(error3 > error)
01601       {
01602          error = error3;
01603       }
01604       if(error4 > error)
01605       {
01606          error = error4;
01607       }
01608       if(error <= tolerance)
01609       {
01610          x += h;
01611 
01612          y_1 += c51*k11+c52*k12+c53*k13+c54*k14+c55*k15+c56*k16;
01613          y_2 += c51*k21+c52*k22+c53*k23+c54*k24+c55*k25+c56*k26;
01614          y_3 += c51*k31+c52*k32+c53*k33+c54*k34+c55*k35+c56*k36;
01615          y_4 += c51*k41+c52*k42+c53*k43+c54*k44+c55*k45+c56*k46;
01616 
01617          tempX[count] = x;
01618          tempY1[count] = y_1;
01619          tempY2[count] = y_2;
01620          tempY3[count] = y_3;
01621          tempY4[count] = y_4;
01622 
01623          count++;
01624 
01625          if(error != 0.0)
01626          {
01627             h *= 0.9*pow(fabs(tolerance/error),0.2);
01628          }
01629 
01630          if(count >= size)
01631          {
01632                          std::cerr << "Flood Error: calculate_Runge_Kutta_Fehlberg_integral() method." << std::endl
01633                        << "Insufficient memory reserved." << std::endl
01634                        << std::endl;
01635 
01636             exit(1);
01637          }
01638       }
01639       else
01640       {
01641          h *= 0.9*pow(fabs(tolerance/error),0.25);      
01642       }
01643    } // end while loop   
01644 
01645    Vector<double> new_x(count);
01646    Vector<double> new_y1(count);
01647    Vector<double> new_y2(count);
01648    Vector<double> new_y3(count);
01649    Vector<double> new_y4(count);
01650 
01651    for(int i = 0; i < count; i++)
01652    {
01653       new_x[i] = tempX[i];
01654       new_y1[i] = tempY1[i];
01655       new_y2[i] = tempY2[i];
01656       new_y3[i] = tempY3[i];
01657       new_y4[i] = tempY4[i];
01658    }
01659 
01660    x_out = new_x;
01661    y1_out = new_y1;
01662    y2_out = new_y2;
01663    y3_out = new_y3;
01664    y4_out = new_y4;
01665 
01666    return(count);
01667 }
01668 
01669 }
01670 
01671 
01672 // Flood: An Open Source Neural Networks C++ Library.
01673 // Copyright (C) 2005-2010 Roberto Lopez 
01674 //
01675 // This library is free software; you can redistribute it and/or
01676 // modify it under the s of the GNU Lesser General Public
01677 // License as published by the Free Software Foundation; either
01678 // version 2.1 of the License, or any later version.
01679 //
01680 // This library is distributed in the hope that it will be useful,
01681 // but WITHOUT ANY WARRANTY; without even the implied warranty of
01682 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
01683 // Lesser General Public License for more details.
01684 
01685 // You should have received a copy of the GNU Lesser General Public
01686 // License along with this library; if not, write to the Free Software
01687 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA