00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include<fstream>
00019 #include<iostream>
00020 #include<string>
00021 #include<sstream>
00022 #include <time.h>
00023
00024
00025
00026 #include "../utilities/vector.h"
00027 #include "../utilities/matrix.h"
00028
00029 #include "ordinary_differential_equations.h"
00030
00031 namespace OpenNN
00032 {
00033
00034
00035
00038
00039 OrdinaryDifferentialEquations::OrdinaryDifferentialEquations(void) : MathematicalModel()
00040 {
00041 set_default();
00042 }
00043
00044
00045
00046
00050
00051 OrdinaryDifferentialEquations::OrdinaryDifferentialEquations(TiXmlElement* ordinary_differential_equations_element)
00052 : MathematicalModel(ordinary_differential_equations_element)
00053 {
00054 set_default();
00055 }
00056
00057
00058
00059
00063
00064 OrdinaryDifferentialEquations::OrdinaryDifferentialEquations(const std::string& filename)
00065 : MathematicalModel(filename)
00066 {
00067 set_default();
00068 }
00069
00070
00071
00072
00076
00077 OrdinaryDifferentialEquations::OrdinaryDifferentialEquations(const OrdinaryDifferentialEquations& other_ordinary_differential_equations)
00078 : MathematicalModel()
00079 {
00080 set(other_ordinary_differential_equations);
00081 }
00082
00083
00084
00085
00088
00089 OrdinaryDifferentialEquations::~OrdinaryDifferentialEquations(void)
00090 {
00091 }
00092
00093
00094
00095
00096
00097
00101
00102 OrdinaryDifferentialEquations& OrdinaryDifferentialEquations::operator = (const OrdinaryDifferentialEquations& other_ordinary_differential_equations)
00103 {
00104 if(this != &other_ordinary_differential_equations)
00105 {
00106 initial_independent_variable = other_ordinary_differential_equations.initial_independent_variable;
00107 final_independent_variable = other_ordinary_differential_equations.final_independent_variable;
00108 initial_dependent_variables = other_ordinary_differential_equations.initial_dependent_variables;
00109 }
00110
00111 return(*this);
00112 }
00113
00114
00115
00116
00117
00118
00123
00124 bool OrdinaryDifferentialEquations::operator == (const OrdinaryDifferentialEquations& other_ordinary_differential_equations) const
00125 {
00126 if(initial_independent_variable == other_ordinary_differential_equations.initial_independent_variable
00127 && final_independent_variable == other_ordinary_differential_equations.final_independent_variable
00128 && initial_dependent_variables == other_ordinary_differential_equations.initial_dependent_variables)
00129 {
00130 return(true);
00131 }
00132 else
00133 {
00134 return(false);
00135 }
00136 }
00137
00138
00139
00140
00141
00142
00144
00145 const double& OrdinaryDifferentialEquations::get_initial_independent_variable(void) const
00146 {
00147 return(initial_independent_variable);
00148 }
00149
00150
00151
00152
00154
00155 const double& OrdinaryDifferentialEquations::get_final_independent_variable(void) const
00156 {
00157 return(final_independent_variable);
00158 }
00159
00160
00161
00162
00164
00165 const Vector<double>& OrdinaryDifferentialEquations::get_initial_dependent_variables(void) const
00166 {
00167 return(initial_dependent_variables);
00168 }
00169
00170
00171
00172
00175
00176 const double& OrdinaryDifferentialEquations::get_initial_dependent_variable(const unsigned int& i) const
00177 {
00178 return(initial_dependent_variables[i]);
00179 }
00180
00181
00182
00183
00185
00186 const OrdinaryDifferentialEquations::SolutionMethod& OrdinaryDifferentialEquations::get_solution_method(void) const
00187 {
00188 return(solution_method);
00189 }
00190
00191
00192
00193
00195
00196 std::string OrdinaryDifferentialEquations::write_solution_method(void) const
00197 {
00198 if(solution_method == RungeKutta)
00199 {
00200 return("RungeKutta");
00201 }
00202 else if(solution_method == RungeKuttaFehlberg)
00203 {
00204 return("RungeKuttaFehlberg");
00205 }
00206 else
00207 {
00208 std::ostringstream buffer;
00209
00210 buffer << "OpenNN Exception: OrdinaryDifferentialEquations class.\n"
00211 << "std::string write_solution_method(void) const method.\n"
00212 << "Unknown solution method.\n";
00213
00214 throw std::logic_error(buffer.str());
00215 }
00216 }
00217
00218
00219
00220
00222
00223 const unsigned int& OrdinaryDifferentialEquations::get_points_number(void) const
00224 {
00225 return(points_number);
00226 }
00227
00228
00229
00230
00232
00233 const double& OrdinaryDifferentialEquations::get_tolerance(void) const
00234 {
00235 return(tolerance);
00236 }
00237
00238
00239
00240
00242
00243 const unsigned int& OrdinaryDifferentialEquations::get_initial_size(void) const
00244 {
00245 return(initial_size);
00246 }
00247
00248
00249
00250
00252
00253 const unsigned int& OrdinaryDifferentialEquations::get_warning_size(void) const
00254 {
00255 return(warning_size);
00256 }
00257
00258
00259
00260
00262
00263 const unsigned int& OrdinaryDifferentialEquations::get_error_size(void) const
00264 {
00265 return(error_size);
00266 }
00267
00268
00269
00270
00273
00274 void OrdinaryDifferentialEquations::set(const OrdinaryDifferentialEquations& other_ordinary_differential_equations)
00275 {
00276 independent_variables_number = other_ordinary_differential_equations.independent_variables_number;
00277 dependent_variables_number = other_ordinary_differential_equations.dependent_variables_number;
00278
00279 initial_independent_variable = other_ordinary_differential_equations.initial_independent_variable;
00280 final_independent_variable = other_ordinary_differential_equations.final_independent_variable;
00281
00282 initial_dependent_variables = other_ordinary_differential_equations.initial_dependent_variables;
00283
00284 solution_method = other_ordinary_differential_equations.solution_method;
00285
00286 points_number = other_ordinary_differential_equations.points_number;
00287
00288 tolerance = other_ordinary_differential_equations.tolerance;
00289
00290 initial_size = other_ordinary_differential_equations.initial_size;
00291 warning_size = other_ordinary_differential_equations.warning_size;
00292 error_size = other_ordinary_differential_equations.error_size;
00293
00294 display = other_ordinary_differential_equations.display;
00295 }
00296
00297
00298
00299
00302
00303 void OrdinaryDifferentialEquations::set_initial_independent_variable(const double& new_initial_independent_variable)
00304 {
00305 initial_independent_variable = new_initial_independent_variable;
00306 }
00307
00308
00309
00310
00313
00314 void OrdinaryDifferentialEquations::set_final_independent_variable(const double& new_final_independent_variable)
00315 {
00316 final_independent_variable = new_final_independent_variable;
00317 }
00318
00319
00320
00321
00324
00325 void OrdinaryDifferentialEquations::set_initial_dependent_variables(const Vector<double>& new_initial_dependent_variables)
00326 {
00327 initial_dependent_variables = new_initial_dependent_variables;
00328 }
00329
00330
00331
00332
00336
00337 void OrdinaryDifferentialEquations::set_initial_dependent_variable(const unsigned int& i, const double& new_initial_dependent_variable)
00338 {
00339 initial_dependent_variables[i] = new_initial_dependent_variable;
00340 }
00341
00342
00343
00344
00348
00349 void OrdinaryDifferentialEquations::set_solution_method(const SolutionMethod& new_solution_method)
00350 {
00351 solution_method = new_solution_method;
00352 }
00353
00354
00355
00356
00360
00361 void OrdinaryDifferentialEquations::set_solution_method(const std::string& new_solution_method)
00362 {
00363 if(new_solution_method == "RungeKutta")
00364 {
00365 set_solution_method(RungeKutta);
00366 }
00367 else if(new_solution_method == "RungeKuttaFehlberg")
00368 {
00369 set_solution_method(RungeKuttaFehlberg);
00370 }
00371 else
00372 {
00373 std::ostringstream buffer;
00374
00375 buffer << "OpenNN Exception: OrdinaryDifferentialEquations class.\n"
00376 << "void set_solution_method(const std::string&) method.\n"
00377 << "Unknown solution method: " << new_solution_method << ".\n";
00378
00379 throw std::logic_error(buffer.str());
00380 }
00381 }
00382
00383
00384
00385
00388
00389 void OrdinaryDifferentialEquations::set_points_number(const unsigned int& new_points_number)
00390 {
00391 points_number = new_points_number;
00392 }
00393
00394
00395
00396
00399
00400 void OrdinaryDifferentialEquations::set_tolerance(const double& new_tolerance)
00401 {
00402 tolerance = new_tolerance;
00403 }
00404
00405
00406
00407
00410
00411 void OrdinaryDifferentialEquations::set_initial_size(const unsigned int& new_initial_size)
00412 {
00413 initial_size = new_initial_size;
00414 }
00415
00416
00417
00418
00421
00422 void OrdinaryDifferentialEquations::set_warning_size(const unsigned int& new_warning_size)
00423 {
00424 warning_size = new_warning_size;
00425 }
00426
00427
00428
00429
00432
00433 void OrdinaryDifferentialEquations::set_error_size(const unsigned int& new_error_size)
00434 {
00435 error_size = new_error_size;
00436 }
00437
00438
00439
00440
00453
00454 void OrdinaryDifferentialEquations::set_default(void)
00455 {
00456 independent_variables_number = 1;
00457 dependent_variables_number = 0;
00458
00459 solution_method = RungeKuttaFehlberg;
00460
00461 points_number = 101;
00462
00463 tolerance = 1.0e-6;
00464
00465 initial_size = (unsigned int)1.0e3;
00466 warning_size = (unsigned int)1.0e6;
00467 error_size = (unsigned int)1.0e9;
00468
00469 display = true;
00470
00471 }
00472
00473
00474
00475
00478
00479 Matrix<double> OrdinaryDifferentialEquations::calculate_Runge_Kutta_solution(const NeuralNetwork& neural_network) const
00480 {
00481 const unsigned int variables_number = count_variables_number();
00482
00483 const double h = (final_independent_variable - initial_independent_variable)/(points_number-1.0);
00484
00485
00486
00487 Vector< Vector<double> > c(4);
00488
00489 Matrix<double> solution(points_number, variables_number);
00490
00491 Vector<double> variables(variables_number);
00492
00493
00494
00495 solution[0][0] = initial_independent_variable;
00496
00497 for(unsigned int j = 0; j < dependent_variables_number; j++)
00498 {
00499 solution[0][1+j] = initial_dependent_variables[j];
00500 }
00501
00502
00503
00504 for(unsigned int i = 0; i < points_number-1; i++)
00505 {
00506 solution[i+1][0] = solution[i][0] + h;
00507
00508
00509
00510 variables[0] = solution[i][0];
00511
00512 for(unsigned int j = 0; j < dependent_variables_number; j++)
00513 {
00514 variables[1+j] = solution[i][1+j];
00515 }
00516
00517 c[0] = calculate_dependent_variables_dots(neural_network, variables);
00518
00519
00520
00521 variables[0] = solution[i][0] + h/2.0;
00522
00523 for(unsigned int j = 0; j < dependent_variables_number; j++)
00524 {
00525 variables[1+j] = solution[i][1+j] + h*c[0][j]/2.0;
00526 }
00527
00528 c[1] = calculate_dependent_variables_dots(neural_network, variables);
00529
00530
00531
00532 variables[0] = solution[i][0] + h/2.0;
00533
00534 for(unsigned int j = 0; j < dependent_variables_number; j++)
00535 {
00536 variables[1+j] = solution[i][1+j] + h*c[1][j]/2.0;
00537 }
00538
00539 c[2] = calculate_dependent_variables_dots(neural_network, variables);
00540
00541
00542
00543 variables[0] = solution[i+1][0];
00544
00545 for(unsigned int j = 0; j < dependent_variables_number; j++)
00546 {
00547 variables[1+j] = solution[i][1+j] + h*c[2][j]/2.0;
00548 }
00549
00550 c[3] = calculate_dependent_variables_dots(neural_network, variables);
00551
00552
00553
00554 for(unsigned int j = 0; j < dependent_variables_number; j++)
00555 {
00556 solution[i+1][1+j] = solution[i][1+j] + h*(c[0][j] + 2.0*c[1][j] + 2.0*c[2][j] + c[3][j])/6.0;
00557 }
00558 }
00559
00560 solution[points_number-1][0] = final_independent_variable;
00561
00562 return(solution);
00563 }
00564
00565
00566
00567
00570
00571 Vector<double> OrdinaryDifferentialEquations::calculate_Runge_Kutta_final_solution(const NeuralNetwork& neural_network) const
00572 {
00573 const unsigned int variables_number = count_variables_number();
00574
00575 const double h = (final_independent_variable - initial_independent_variable)/(points_number-1.0);
00576
00577
00578
00579 Vector< Vector<double> > c(4);
00580
00581 Vector<double> final_solution(variables_number);
00582
00583 Vector<double> variables(variables_number);
00584
00585
00586
00587 final_solution[0] = initial_independent_variable;
00588
00589 for(unsigned int j = 0; j < dependent_variables_number; j++)
00590 {
00591 final_solution[0] = initial_dependent_variables[j];
00592 }
00593
00594
00595
00596 for(unsigned int i = 0; i < points_number-1; i++)
00597 {
00598 final_solution[0] = final_solution[0] + h;
00599
00600
00601
00602 variables[0] = final_solution[0];
00603
00604 for(unsigned int j = 0; j < dependent_variables_number; j++)
00605 {
00606 variables[1+j] = final_solution[1+j];
00607 }
00608
00609 c[0] = calculate_dependent_variables_dots(neural_network, variables);
00610
00611
00612
00613 variables[0] = final_solution[0] + h/2.0;
00614
00615 for(unsigned int j = 0; j < dependent_variables_number; j++)
00616 {
00617 variables[1+j] = final_solution[1+j] + h*c[0][j]/2.0;
00618 }
00619
00620 c[1] = calculate_dependent_variables_dots(neural_network, variables);
00621
00622
00623
00624 variables[0] = final_solution[0] + h/2.0;
00625
00626 for(unsigned int j = 0; j < dependent_variables_number; j++)
00627 {
00628 variables[1+j] = final_solution[1+j] + h*c[1][j]/2.0;
00629 }
00630
00631 c[2] = calculate_dependent_variables_dots(neural_network, variables);
00632
00633
00634
00635 variables[0] = final_solution[0];
00636
00637 for(unsigned int j = 0; j < dependent_variables_number; j++)
00638 {
00639 variables[1+j] = final_solution[1+j] + h*c[2][j]/2.0;
00640 }
00641
00642 c[3] = calculate_dependent_variables_dots(neural_network, variables);
00643
00644
00645
00646 for(unsigned int j = 0; j < dependent_variables_number; j++)
00647 {
00648 final_solution[1+j] = final_solution[1+j] + h*(c[0][j] + 2.0*c[1][j] + 2.0*c[2][j] + c[3][j])/6.0;
00649 }
00650 }
00651
00652 return(final_solution);
00653 }
00654
00655
00656
00657
00660
00661 Matrix<double> OrdinaryDifferentialEquations::calculate_Runge_Kutta_Fehlberg_solution(const NeuralNetwork& neural_network) const
00662 {
00663 const double epsilon = 1.0e-12;
00664
00665 const double a2 = 1.0/5.0;
00666 const double a3 = 3.0/10.0;
00667 const double a4 = 3.0/5.0;
00668 const double a5 = 1.0;
00669 const double a6 = 7.0/8.0;
00670
00671 const double b21 = 1.0/5.0;
00672 const double b31 = 3.0/40.0;
00673 const double b32 = 9.0/40.0;
00674 const double b41 = 3.0/10.0;
00675 const double b42 = -9.0/10.0;
00676 const double b43 = 6.0/5.0;
00677 const double b51 = -11.0/54.0;
00678 const double b52 = 5.0/2.0;
00679 const double b53 = -70.0/27.0;
00680 const double b54 = 35.0/27.0;
00681 const double b61 = 1631.0/55296.0;
00682 const double b62 = 175.0/512.0;
00683 const double b63 = 575.0/13824.0;
00684 const double b64 = 44275.0/110592.0;
00685 const double b65 = 253.0/4096.0;
00686
00687 const double c41 = 37.0/378.0;
00688 const double c42 = 0.0;
00689 const double c43 = 250.0/621.0;
00690 const double c44 = 125.0/594.0;
00691 const double c45 = 0.0;
00692 const double c46 = 512.0/1771.0;
00693
00694 const double c51 = 2825.0/27648.0;
00695 const double c52 = 0.0;
00696 const double c53 = 18575.0/48384.0;
00697 const double c54 = 13525.0/55296.0;
00698 const double c55 = 277.0/14336.0;
00699 const double c56 = 1.0/4.0;
00700
00701 const double d1 = c41 - c51;
00702 const double d2 = c42 - c52;
00703 const double d3 = c43 - c53;
00704 const double d4 = c44 - c54;
00705 const double d5 = c45 - c55;
00706 const double d6 = c46 - c56;
00707
00708 const unsigned int variables_number = count_variables_number();
00709
00710 unsigned int size = initial_size;
00711
00712 Vector<double> errors(dependent_variables_number, 0.0);
00713
00714
00715 double error = 0.0;
00716
00717 Vector< Vector<double> > c(6);
00718
00719 double hmin = 0.0;
00720 double h = (final_independent_variable - initial_independent_variable)*1.0e-3;
00721
00722 Vector<double> variables(variables_number);
00723
00724 Matrix<double> solution(size, variables_number);
00725
00726 unsigned int point_index = 0;
00727
00728
00729
00730 solution[0][0] = initial_independent_variable;
00731
00732 for(unsigned int j = 0; j < dependent_variables_number; j++)
00733 {
00734 solution[0][1+j] = initial_dependent_variables[j];
00735 }
00736
00737
00738
00739 while(solution[point_index][0] < final_independent_variable)
00740 {
00741
00742
00743 hmin = 32.0*epsilon*fabs(solution[point_index][0]);
00744
00745 if(h < hmin)
00746 {
00747 if(display)
00748 {
00749 std::cout << "OpenNN Warning: OrdinaryDifferentialEquations class.\n"
00750 << "calculate_Runge_Kutta_Fehlberg_solution() method.\n"
00751 << "Step size is less than smallest allowable." << std::endl;
00752 }
00753
00754 h = hmin;
00755 }
00756
00757 if(solution[point_index][0] + h > final_independent_variable)
00758 {
00759 h = final_independent_variable - solution[point_index][0];
00760 }
00761
00762
00763
00764 variables[0] = solution[point_index][0];
00765
00766 for(unsigned int j = 0; j < dependent_variables_number; j++)
00767 {
00768 variables[1+j] = solution[point_index][1+j];
00769 }
00770
00771 c[0] = calculate_dependent_variables_dots(neural_network, variables)*h;
00772
00773
00774
00775 variables[0] = solution[point_index][0] + a2*h;
00776
00777 for(unsigned int j = 0; j < dependent_variables_number; j++)
00778 {
00779 variables[1+j] = solution[point_index][1+j]+ b21*c[0][j];
00780 }
00781
00782 c[1] = calculate_dependent_variables_dots(neural_network, variables)*h;
00783
00784
00785
00786 variables[0] = solution[point_index][0] + a3*h;
00787
00788 for(unsigned int j = 0; j < dependent_variables_number; j++)
00789 {
00790 variables[1+j] = solution[point_index][1+j] + b31*c[0][j] + b32*c[1][j];
00791 }
00792
00793 c[2] = calculate_dependent_variables_dots(neural_network, variables)*h;
00794
00795
00796
00797 variables[0] = solution[point_index][0] + a4*h;
00798
00799 for(unsigned int j = 0; j < dependent_variables_number; j++)
00800 {
00801 variables[1+j] = solution[point_index][1+j] + b41*c[0][j] + b42*c[1][j] + b43*c[2][j];
00802 }
00803
00804 c[3] = calculate_dependent_variables_dots(neural_network, variables)*h;
00805
00806
00807
00808 variables[0] = solution[point_index][0] + a5*h;
00809
00810 for(unsigned int j = 0; j < dependent_variables_number; j++)
00811 {
00812 variables[1+j] = solution[point_index][1+j] + b51*c[0][j] + b52*c[1][j] + b53*c[2][j] + b54*c[3][j];
00813 }
00814
00815 c[4] = calculate_dependent_variables_dots(neural_network, variables)*h;
00816
00817
00818
00819 variables[0] = solution[point_index][0] + a6*h;
00820
00821 for(unsigned int j = 0; j < dependent_variables_number; j++)
00822 {
00823 variables[1+j] = solution[point_index][1+j] + b61*c[0][j] + b62*c[1][j] + b63*c[2][j] + b64*c[3][j] + b65*c[4][j];
00824 }
00825
00826 c[5] = calculate_dependent_variables_dots(neural_network, variables)*h;
00827
00828
00829
00830 for(unsigned int j = 0; j < dependent_variables_number; j++)
00831 {
00832 errors[j] = fabs(d1*c[0][j] + d2*c[1][j] + d3*c[2][j] + d4*c[3][j] + d5*c[4][j] + d6*c[5][j]);
00833 }
00834
00835 error = errors.calculate_maximum();
00836
00837 if(error <= tolerance)
00838 {
00839 solution[point_index+1][0] = solution[point_index][0] + h;
00840
00841 for(unsigned int j = 0; j < dependent_variables_number; j++)
00842 {
00843 solution[point_index+1][1+j] = solution[point_index][1+j] + c51*c[0][j] + c52*c[1][j] + c53*c[2][j] + c54*c[3][j] + c55*c[4][j] + c56*c[5][j];
00844 }
00845
00846 point_index++;
00847
00848 if(error != 0.0)
00849 {
00850 h *= 0.9*pow(fabs(tolerance/error), 0.2);
00851 }
00852
00853 if(point_index >= size)
00854 {
00855 size *= 2;
00856
00857 if(display && size > warning_size)
00858 {
00859 std::cout << "OpenNN Warning: OrdinaryDifferentialEquations class." << std::endl
00860 << "calculate_Runge_Kutta_Fehlberg_solution() method." << std::endl
00861 << "Solution size is " << size << std::endl;
00862 }
00863 else if(size > error_size)
00864 {
00865 std::ostringstream buffer;
00866
00867 buffer << "OpenNN Exception: OrdinaryDifferentialEquations class." << std::endl
00868 << "calculate_Runge_Kutta_Fehlberg_solution() method." << std::endl
00869 << "Solution size is bigger than greatest allowable." << std::endl;
00870
00871 throw std::logic_error(buffer.str().c_str());
00872 }
00873
00874 solution.resize(size, variables_number);
00875 }
00876 }
00877 else
00878 {
00879 h *= 0.9*pow(fabs(tolerance/error), 0.25);
00880 }
00881 }
00882
00883 solution.resize(point_index+1, variables_number);
00884
00885 return(solution);
00886 }
00887
00888
00889
00890
00893
00894 Vector<double> OrdinaryDifferentialEquations::calculate_Runge_Kutta_Fehlberg_final_solution(const NeuralNetwork& neural_network) const
00895 {
00896 const double epsilon = 1.0e-12;
00897
00898 const double a2 = 1.0/5.0;
00899 const double a3 = 3.0/10.0;
00900 const double a4 = 3.0/5.0;
00901 const double a5 = 1.0;
00902 const double a6 = 7.0/8.0;
00903
00904 const double b21 = 1.0/5.0;
00905 const double b31 = 3.0/40.0;
00906 const double b32 = 9.0/40.0;
00907 const double b41 = 3.0/10.0;
00908 const double b42 = -9.0/10.0;
00909 const double b43 = 6.0/5.0;
00910 const double b51 = -11.0/54.0;
00911 const double b52 = 5.0/2.0;
00912 const double b53 = -70.0/27.0;
00913 const double b54 = 35.0/27.0;
00914 const double b61 = 1631.0/55296.0;
00915 const double b62 = 175.0/512.0;
00916 const double b63 = 575.0/13824.0;
00917 const double b64 = 44275.0/110592.0;
00918 const double b65 = 253.0/4096.0;
00919
00920 const double c41 = 37.0/378.0;
00921 const double c42 = 0.0;
00922 const double c43 = 250.0/621.0;
00923 const double c44 = 125.0/594.0;
00924 const double c45 = 0.0;
00925 const double c46 = 512.0/1771.0;
00926
00927 const double c51 = 2825.0/27648.0;
00928 const double c52 = 0.0;
00929 const double c53 = 18575.0/48384.0;
00930 const double c54 = 13525.0/55296.0;
00931 const double c55 = 277.0/14336.0;
00932 const double c56 = 1.0/4.0;
00933
00934 const double d1 = c41 - c51;
00935 const double d2 = c42 - c52;
00936 const double d3 = c43 - c53;
00937 const double d4 = c44 - c54;
00938 const double d5 = c45 - c55;
00939 const double d6 = c46 - c56;
00940
00941 const unsigned int variables_number = count_variables_number();
00942
00943 Vector<double> errors(dependent_variables_number, 0.0);
00944
00945 double error = 0.0;
00946
00947 Vector< Vector<double> > c(6);
00948
00949 double hmin = 0.0;
00950 double h = (final_independent_variable - initial_independent_variable)*1.0e-3;
00951
00952 Vector<double> variables(variables_number);
00953
00954 Vector<double> final_solution(variables_number);
00955
00956
00957
00958 final_solution[0] = initial_independent_variable;
00959
00960 for(unsigned int j = 0; j < dependent_variables_number; j++)
00961 {
00962 final_solution[1+j] = initial_dependent_variables[j];
00963 }
00964
00965
00966
00967 while(final_solution[0] < final_independent_variable)
00968 {
00969
00970
00971 hmin = 32.0*epsilon*fabs(final_solution[0]);
00972
00973 if(h < hmin)
00974 {
00975 if(display)
00976 {
00977 std::cout << "OpenNN Warning: OrdinaryDifferentialEquations class.\n"
00978 << "calculate_Runge_Kutta_Fehlberg_solution() method.\n"
00979 << "Step size is less than smallest allowable." << std::endl;
00980 }
00981
00982 h = hmin;
00983 }
00984
00985 if(final_solution[0] + h > final_independent_variable)
00986 {
00987 h = final_independent_variable - final_solution[0];
00988 }
00989
00990
00991
00992 variables[0] = final_solution[0];
00993
00994 for(unsigned int j = 0; j < dependent_variables_number; j++)
00995 {
00996 variables[1+j] = final_solution[1+j];
00997 }
00998
00999 c[0] = calculate_dependent_variables_dots(neural_network, variables)*h;
01000
01001
01002
01003 variables[0] = final_solution[0] + a2*h;
01004
01005 for(unsigned int j = 0; j < dependent_variables_number; j++)
01006 {
01007 variables[1+j] = final_solution[1+j]+ b21*c[0][j];
01008 }
01009
01010 c[1] = calculate_dependent_variables_dots(neural_network, variables)*h;
01011
01012
01013
01014 variables[0] = final_solution[0] + a3*h;
01015
01016 for(unsigned int j = 0; j < dependent_variables_number; j++)
01017 {
01018 variables[1+j] = final_solution[1+j] + b31*c[0][j] + b32*c[1][j];
01019 }
01020
01021 c[2] = calculate_dependent_variables_dots(neural_network, variables)*h;
01022
01023
01024
01025 variables[0] = final_solution[0] + a4*h;
01026
01027 for(unsigned int j = 0; j < dependent_variables_number; j++)
01028 {
01029 variables[1+j] = final_solution[1+j] + b41*c[0][j] + b42*c[1][j] + b43*c[2][j];
01030 }
01031
01032 c[3] = calculate_dependent_variables_dots(neural_network, variables)*h;
01033
01034
01035
01036 variables[0] = final_solution[0] + a5*h;
01037
01038 for(unsigned int j = 0; j < dependent_variables_number; j++)
01039 {
01040 variables[1+j] = final_solution[1+j] + b51*c[0][j] + b52*c[1][j] + b53*c[2][j] + b54*c[3][j];
01041 }
01042
01043 c[4] = calculate_dependent_variables_dots(neural_network, variables)*h;
01044
01045
01046
01047 variables[0] = final_solution[0] + a6*h;
01048
01049 for(unsigned int j = 0; j < dependent_variables_number; j++)
01050 {
01051 variables[1+j] = final_solution[1+j] + b61*c[0][j] + b62*c[1][j] + b63*c[2][j] + b64*c[3][j] + b65*c[4][j];
01052 }
01053
01054 c[5] = calculate_dependent_variables_dots(neural_network, variables)*h;
01055
01056
01057
01058 for(unsigned int j = 0; j < dependent_variables_number; j++)
01059 {
01060 errors[j] = fabs(d1*c[0][j] + d2*c[1][j] + d3*c[2][j] + d4*c[3][j] + d5*c[4][j] + d6*c[5][j]);
01061 }
01062
01063 error = errors.calculate_maximum();
01064
01065 if(error <= tolerance)
01066 {
01067 final_solution[0] += h;
01068
01069 for(unsigned int j = 0; j < dependent_variables_number; j++)
01070 {
01071 final_solution[1+j] += c51*c[0][j] + c52*c[1][j] + c53*c[2][j] + c54*c[3][j] + c55*c[4][j] + c56*c[5][j];
01072 }
01073
01074 if(error != 0.0)
01075 {
01076 h *= 0.9*pow(fabs(tolerance/error), 0.2);
01077 }
01078 }
01079 else
01080 {
01081 h *= 0.9*pow(fabs(tolerance/error), 0.25);
01082 }
01083 }
01084
01085
01086 return(final_solution);
01087 }
01088
01089
01090
01091
01092 Matrix<double> OrdinaryDifferentialEquations::calculate_solutions(const NeuralNetwork& neural_network) const
01093 {
01094 switch(solution_method)
01095 {
01096 case RungeKutta:
01097 {
01098 return(calculate_Runge_Kutta_solution(neural_network));
01099 }
01100 break;
01101
01102 case RungeKuttaFehlberg:
01103 {
01104 return(calculate_Runge_Kutta_Fehlberg_solution(neural_network));
01105 }
01106 break;
01107
01108 default:
01109 {
01110 std::ostringstream buffer;
01111
01112 buffer << "OpenNN Error: OrdinaryDifferentialEquations class\n"
01113 << "Vector<double> calculate_solutions(const NeuralNetwork&) const method.\n"
01114 << "Unknown solution method.\n";
01115
01116 throw std::logic_error(buffer.str());
01117 }
01118 break;
01119 }
01120 }
01121
01122
01123
01124
01125 Vector<double> OrdinaryDifferentialEquations::calculate_final_solutions(const NeuralNetwork& neural_network) const
01126 {
01127 switch(solution_method)
01128 {
01129 case RungeKutta:
01130 {
01131 return(calculate_Runge_Kutta_final_solution(neural_network));
01132 }
01133 break;
01134
01135 case RungeKuttaFehlberg:
01136 {
01137 return(calculate_Runge_Kutta_Fehlberg_final_solution(neural_network));
01138 }
01139 break;
01140
01141 default:
01142 {
01143 std::ostringstream buffer;
01144
01145 buffer << "OpenNN Exception: ScalingLayer class\n"
01146 << "Vector<double> calculate_final_solutions(const NeuralNetwork&) const method.\n"
01147 << "Unknown solution method.\n";
01148
01149 throw std::logic_error(buffer.str());
01150 }
01151 break;
01152 }
01153 }
01154
01155
01156
01157
01159
01160 std::string OrdinaryDifferentialEquations::to_string(void) const
01161 {
01162 std::ostringstream buffer;
01163
01164 buffer << "Mathematical model\n"
01165 << "Independent variables number: " << independent_variables_number << "\n"
01166 << "Dependent variables number: " << dependent_variables_number << "\n"
01167 << "Initial independent variable: " << initial_independent_variable << "\n"
01168 << "Final independent variable: " << final_independent_variable << "\n"
01169 << "Initial dependent variables: " << initial_dependent_variables << "\n"
01170 << "Solution method: " << write_solution_method() << "\n"
01171 << "Points number: " << points_number << "\n"
01172 << "Tolerance" << tolerance << "\n"
01173 << "Initial size: " << initial_size << "\n"
01174 << "Warning size: " << warning_size << "\n"
01175 << "Error size: " << error_size << "\n"
01176 << "Display: " << display << std::endl;
01177
01178 return(buffer.str());
01179 }
01180
01181
01182
01183
01186
01187 TiXmlElement* OrdinaryDifferentialEquations::to_XML(void) const
01188 {
01189 std::ostringstream buffer;
01190
01191 TiXmlElement* ordinary_differential_equations_element = new TiXmlElement("OrdinaryDifferentialEquations");
01192 ordinary_differential_equations_element->SetAttribute("Version", 4);
01193
01194
01195 {
01196 TiXmlElement* element = new TiXmlElement("IndependentVariablesNumber");
01197 ordinary_differential_equations_element->LinkEndChild(element);
01198
01199 buffer.str("");
01200 buffer << independent_variables_number;
01201
01202 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01203 element->LinkEndChild(text);
01204 }
01205
01206
01207 {
01208 TiXmlElement* element = new TiXmlElement("DependentVariablesNumber");
01209 ordinary_differential_equations_element->LinkEndChild(element);
01210
01211 buffer.str("");
01212 buffer << dependent_variables_number;
01213
01214 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01215 element->LinkEndChild(text);
01216 }
01217
01218
01219 {
01220 TiXmlElement* element = new TiXmlElement("InitialIndependentVariable");
01221 ordinary_differential_equations_element->LinkEndChild(element);
01222
01223 buffer.str("");
01224 buffer << initial_independent_variable;
01225
01226 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01227 element->LinkEndChild(text);
01228 }
01229
01230
01231 {
01232 TiXmlElement* element = new TiXmlElement("FinalIndependentVariable");
01233 ordinary_differential_equations_element->LinkEndChild(element);
01234
01235 buffer.str("");
01236 buffer << final_independent_variable;
01237
01238 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01239 element->LinkEndChild(text);
01240 }
01241
01242
01243 {
01244 TiXmlElement* element = new TiXmlElement("InitialDependentVariables");
01245 ordinary_differential_equations_element->LinkEndChild(element);
01246
01247 buffer.str("");
01248 buffer << initial_dependent_variables;
01249
01250 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01251 element->LinkEndChild(text);
01252 }
01253
01254
01255 {
01256 TiXmlElement* element = new TiXmlElement("SolutionMethod");
01257 ordinary_differential_equations_element->LinkEndChild(element);
01258
01259 TiXmlText* text = new TiXmlText(write_solution_method().c_str());
01260 element->LinkEndChild(text);
01261 }
01262
01263
01264 {
01265 TiXmlElement* element = new TiXmlElement("PointsNumber");
01266 ordinary_differential_equations_element->LinkEndChild(element);
01267
01268 buffer.str("");
01269 buffer << points_number;
01270
01271 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01272 element->LinkEndChild(text);
01273 }
01274
01275
01276 {
01277 TiXmlElement* element = new TiXmlElement("Tolerance");
01278 ordinary_differential_equations_element->LinkEndChild(element);
01279
01280 buffer.str("");
01281 buffer << tolerance;
01282
01283 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01284 element->LinkEndChild(text);
01285 }
01286
01287
01288 {
01289 TiXmlElement* element = new TiXmlElement("InitialSize");
01290 ordinary_differential_equations_element->LinkEndChild(element);
01291
01292 buffer.str("");
01293 buffer << initial_size;
01294
01295 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01296 element->LinkEndChild(text);
01297 }
01298
01299
01300 {
01301 TiXmlElement* element = new TiXmlElement("WarningSize");
01302 ordinary_differential_equations_element->LinkEndChild(element);
01303
01304 buffer.str("");
01305 buffer << warning_size;
01306
01307 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01308 element->LinkEndChild(text);
01309 }
01310
01311
01312 {
01313 TiXmlElement* element = new TiXmlElement("ErrorSize");
01314 ordinary_differential_equations_element->LinkEndChild(element);
01315
01316 buffer.str("");
01317 buffer << error_size;
01318
01319 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01320 element->LinkEndChild(text);
01321 }
01322
01323
01324 {
01325 TiXmlElement* element = new TiXmlElement("Display");
01326 ordinary_differential_equations_element->LinkEndChild(element);
01327
01328 buffer.str("");
01329 buffer << display;
01330
01331 TiXmlText* text = new TiXmlText(buffer.str().c_str());
01332 element->LinkEndChild(text);
01333 }
01334
01335 return(ordinary_differential_equations_element);
01336 }
01337
01338
01339
01340
01343
01344 void OrdinaryDifferentialEquations::from_XML(TiXmlElement* ordinary_differential_equations_element)
01345 {
01346 if(ordinary_differential_equations_element)
01347 {
01348
01349 {
01350 TiXmlElement* element = ordinary_differential_equations_element->FirstChildElement("DependentVariablesNumber");
01351
01352 if(element)
01353 {
01354 const char* text = element->GetText();
01355
01356 if(text)
01357 {
01358 try
01359 {
01360 set_dependent_variables_number(atoi(text));
01361 }
01362 catch(std::exception& e)
01363 {
01364 std::cout << e.what() << std::endl;
01365 }
01366 }
01367 }
01368 }
01369
01370
01371 {
01372 TiXmlElement* element = ordinary_differential_equations_element->FirstChildElement("InitialIndependentVariable");
01373
01374 if(element)
01375 {
01376 const char* text = element->GetText();
01377
01378 if(text)
01379 {
01380 try
01381 {
01382 set_initial_independent_variable(atof(text));
01383 }
01384 catch(std::exception& e)
01385 {
01386 std::cout << e.what() << std::endl;
01387 }
01388 }
01389 }
01390 }
01391
01392
01393 {
01394 TiXmlElement* element = ordinary_differential_equations_element->FirstChildElement("FinalIndependentVariable");
01395
01396 if(element)
01397 {
01398 const char* text = element->GetText();
01399
01400 if(text)
01401 {
01402 try
01403 {
01404 set_final_independent_variable(atof(text));
01405 }
01406 catch(std::exception& e)
01407 {
01408 std::cout << e.what() << std::endl;
01409 }
01410 }
01411 }
01412 }
01413
01414
01415 {
01416 TiXmlElement* element = ordinary_differential_equations_element->FirstChildElement("InitialDependentVariables");
01417
01418 if(element)
01419 {
01420 const char* text = element->GetText();
01421
01422 if(text)
01423 {
01424 try
01425 {
01426 Vector<double> new_initial_dependent_variables;
01427 new_initial_dependent_variables.parse(text);
01428
01429 set_initial_dependent_variables(new_initial_dependent_variables);
01430 }
01431 catch(std::exception& e)
01432 {
01433 std::cout << e.what() << std::endl;
01434 }
01435 }
01436 }
01437 }
01438
01439
01440 {
01441 TiXmlElement* element = ordinary_differential_equations_element->FirstChildElement("SolutionMethod");
01442
01443 if(element)
01444 {
01445 const char* text = element->GetText();
01446
01447 if(text)
01448 {
01449 try
01450 {
01451 std::string new_solution_method(text);
01452
01453 set_solution_method(new_solution_method);
01454 }
01455 catch(std::exception& e)
01456 {
01457 std::cout << e.what() << std::endl;
01458 }
01459 }
01460 }
01461 }
01462
01463
01464 {
01465 TiXmlElement* element = ordinary_differential_equations_element->FirstChildElement("PointsNumber");
01466
01467 if(element)
01468 {
01469 const char* text = element->GetText();
01470
01471 if(text)
01472 {
01473 try
01474 {
01475 set_points_number(atoi(text));
01476 }
01477 catch(std::exception& e)
01478 {
01479 std::cout << e.what() << std::endl;
01480 }
01481 }
01482 }
01483 }
01484
01485
01486 {
01487 TiXmlElement* element = ordinary_differential_equations_element->FirstChildElement("Tolerance");
01488
01489 if(element)
01490 {
01491 const char* text = element->GetText();
01492
01493 if(text)
01494 {
01495 try
01496 {
01497 set_tolerance(atof(text));
01498 }
01499 catch(std::exception& e)
01500 {
01501 std::cout << e.what() << std::endl;
01502 }
01503 }
01504 }
01505 }
01506
01507
01508 {
01509 TiXmlElement* element = ordinary_differential_equations_element->FirstChildElement("InitialSize");
01510
01511 if(element)
01512 {
01513 const char* text = element->GetText();
01514
01515 if(text)
01516 {
01517 try
01518 {
01519 set_initial_size(atoi(text));
01520 }
01521 catch(std::exception& e)
01522 {
01523 std::cout << e.what() << std::endl;
01524 }
01525 }
01526 }
01527 }
01528
01529
01530 {
01531 TiXmlElement* element = ordinary_differential_equations_element->FirstChildElement("WarningSize");
01532
01533 if(element)
01534 {
01535 const char* text = element->GetText();
01536
01537 if(text)
01538 {
01539 try
01540 {
01541 set_warning_size(atoi(text));
01542 }
01543 catch(std::exception& e)
01544 {
01545 std::cout << e.what() << std::endl;
01546 }
01547 }
01548 }
01549 }
01550
01551
01552 {
01553 TiXmlElement* element = ordinary_differential_equations_element->FirstChildElement("ErrorSize");
01554
01555 if(element)
01556 {
01557 const char* text = element->GetText();
01558
01559 if(text)
01560 {
01561 try
01562 {
01563 set_error_size(atoi(text));
01564 }
01565 catch(std::exception& e)
01566 {
01567 std::cout << e.what() << std::endl;
01568 }
01569 }
01570 }
01571 }
01572
01573
01574 {
01575 TiXmlElement* display_element = ordinary_differential_equations_element->FirstChildElement("Display");
01576
01577 if(display_element)
01578 {
01579 const char* display_text = display_element->GetText();
01580
01581 if(display_text)
01582 {
01583 try
01584 {
01585 std::string display_string(display_text);
01586
01587 set_display(display_string != "0");
01588 }
01589 catch(std::exception& e)
01590 {
01591 std::cout << e.what() << std::endl;
01592 }
01593 }
01594 }
01595 }
01596 }
01597 }
01598
01599
01600
01601
01602 void OrdinaryDifferentialEquations::save_data(const NeuralNetwork& neural_network, const std::string& filename) const
01603 {
01604 const Matrix<double> solution = calculate_Runge_Kutta_solution(neural_network);
01605
01606 solution.save(filename);
01607 }
01608
01609 }
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627