00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include<iostream>
00019
00020
00021
00022 #include "numerical_integration.h"
00023
00024
00025
00026 #include "../../parsers/tinyxml/tinyxml.h"
00027
00028
00029 namespace OpenNN
00030 {
00031
00032
00033
00035
00036 NumericalIntegration::NumericalIntegration(void)
00037 {
00038
00039 }
00040
00041
00042
00043
00045
00046 NumericalIntegration::~NumericalIntegration(void)
00047 {
00048
00049 }
00050
00051
00052
00053
00054
00055
00057
00058 const NumericalIntegration::NumericalIntegrationMethod& NumericalIntegration::get_numerical_integration_method(void) const
00059 {
00060 return(numerical_integration_method);
00061 }
00062
00063
00064
00065
00067
00068 std::string NumericalIntegration::write_numerical_integration_method(void) const
00069 {
00070 switch(numerical_integration_method)
00071 {
00072 case TrapezoidMethod:
00073 {
00074 return("TrapezoidMethod");
00075 }
00076 break;
00077
00078 case SimpsonMethod:
00079 {
00080 return("SimpsonMethod");
00081 }
00082 break;
00083
00084 default:
00085 {
00086 std::ostringstream buffer;
00087
00088 buffer << "OpenNN Exception: NumericalIntegration class.\n"
00089 << "std::string write_numerical_integration_method(void) const method.\n"
00090 << "Unknown numerical integration method.\n";
00091
00092 throw std::logic_error(buffer.str());
00093 }
00094 break;
00095 }
00096 }
00097
00098
00099
00100
00102
00103 const bool& NumericalIntegration::get_display(void) const
00104 {
00105 return(display);
00106 }
00107
00108
00109
00110
00113
00114 void NumericalIntegration::set(const NumericalIntegration& other_numerical_integration)
00115 {
00116 numerical_integration_method = other_numerical_integration.numerical_integration_method;
00117
00118 display = other_numerical_integration.display;
00119 }
00120
00121
00122
00123
00126
00127 void NumericalIntegration::set_numerical_integration_method
00128 (const NumericalIntegration::NumericalIntegrationMethod& new_numerical_integration_method)
00129 {
00130 numerical_integration_method = new_numerical_integration_method;
00131 }
00132
00133
00134
00135
00139
00140 void NumericalIntegration::set_numerical_integration_method(const std::string& new_numerical_integration_method)
00141 {
00142 if(new_numerical_integration_method == "TrapezoidMethod")
00143 {
00144 numerical_integration_method = TrapezoidMethod;
00145 }
00146 else if(new_numerical_integration_method == "SimpsonMethod")
00147 {
00148 numerical_integration_method = SimpsonMethod;
00149 }
00150 else
00151 {
00152 std::ostringstream buffer;
00153
00154 buffer << "OpenNN Exception: NumericalIntegration class.\n"
00155 << "void set_numerical_integration_method(const std::string&) method.\n"
00156 << "Unknown numerical integration method: " << new_numerical_integration_method << ".\n";
00157
00158 throw std::logic_error(buffer.str());
00159 }
00160 }
00161
00162
00163
00164
00167
00168 void NumericalIntegration::set_display(const bool& new_display)
00169 {
00170 display = new_display;
00171 }
00172
00173
00174
00175
00181
00182 void NumericalIntegration::set_default(void)
00183 {
00184 numerical_integration_method = SimpsonMethod;
00185
00186 display = true;
00187 }
00188
00189
00190
00191
00196
00197 double NumericalIntegration::calculate_trapezoid_integral(const Vector<double>& x, const Vector<double>& y) const
00198 {
00199
00200
00201 unsigned int n = x.size();
00202
00203
00204
00205 double trapezoid_integral = 0;
00206
00207 for(unsigned int i = 0; i < n-1; i++)
00208 {
00209 trapezoid_integral += 0.5*(x[i+1]-x[i])*(y[i+1]+y[i]);
00210 }
00211
00212
00213
00214 return(trapezoid_integral);
00215 }
00216
00217
00218
00219
00224
00225 double NumericalIntegration::calculate_Simpson_integral(const Vector<double>& x, const Vector<double>& y) const
00226 {
00227 double Simpson_integral = 0.0;
00228
00229
00230
00231 unsigned int n = x.size();
00232
00233 unsigned int m = 0;
00234
00235 double a = 0.0;
00236 double fa = 0.0;
00237 double b = 0.0;
00238 double fb = 0.0;
00239 double c = 0.0;
00240 double fc = 0.0;
00241 double wa = 0.0;
00242 double wb = 0.0;
00243 double wc = 0.0;
00244 double h = 0.0;
00245
00246 double sum = 0.0;
00247
00248 if(n%2 != 0)
00249 {
00250 m=(n-1)/2;
00251
00252 for(unsigned int i = 0 ; i < m ; i++ )
00253 {
00254 a = x[2*i];
00255 b = x[2*i+1];
00256 c = x[2*i+2];
00257
00258 fa = y[2*i];
00259 fb = y[2*i+1];
00260 fc = y[2*i+2];
00261
00262 wa = (c-a)/((a-b)*(a-c))*(1.0/3.0*(a*a+c*c+a*c)-0.5*(a+c)*(b+c)+b*c);
00263 wb = (c-a)/((b-a)*(b-c))*(1.0/3.0*(a*a+c*c+a*c)-0.5*(a+c)*(a+c)+a*c);
00264 wc = (c-a)/((c-a)*(c-b))*(1.0/3.0*(a*a+c*c+a*c)-0.5*(a+c)*(a+b)+a*b);
00265
00266 sum += wa*fa+wb*fb+wc*fc;
00267 }
00268 }
00269 else
00270 {
00271 m=(n-2)/2;
00272
00273 for(unsigned int i = 0; i < m; i++ )
00274 {
00275 a = x[2*i];
00276 b = x[2*i+1];
00277 c = x[2*i+2];
00278
00279 fa = y[2*i];
00280 fb = y[2*i+1];
00281 fc = y[2*i+2];
00282
00283 wa = (c-a)/((a-b)*(a-c))*(1.0/3.0*(a*a+c*c+a*c)-0.5*(a+c)*(b+c)+b*c);
00284 wb = (c-a)/((b-a)*(b-c))*(1.0/3.0*(a*a+c*c+a*c)-0.5*(a+c)*(a+c)+a*c);
00285 wc = (c-a)/((c-a)*(c-b))*(1.0/3.0*(a*a+c*c+a*c)-0.5*(a+c)*(a+b)+a*b);
00286
00287 sum += wa*fa+wb*fb+wc*fc;
00288 }
00289
00290
00291
00292 h = x[n-1]-x[n-2];
00293
00294 sum += h*(y[n-1]+y[n-2])/2.0;
00295 }
00296
00297 Simpson_integral = sum ;
00298
00299 return(Simpson_integral);
00300 }
00301
00302
00303
00304
00308
00309 double NumericalIntegration::calculate_integral(const Vector<double>& x, const Vector<double>& y) const
00310 {
00311 switch(numerical_integration_method)
00312 {
00313 case TrapezoidMethod:
00314 {
00315 return(calculate_trapezoid_integral(x, y));
00316 }
00317 break;
00318
00319 case SimpsonMethod:
00320 {
00321 return(calculate_Simpson_integral(x, y));
00322 }
00323 break;
00324
00325 default:
00326 {
00327 std::ostringstream buffer;
00328
00329 buffer << "OpenNN Exception: NumericalIntegration class.\n"
00330 << "double calculate_integral(const Vector<double>&, const Vector<double>&) const method.\n"
00331 << "Unknown numerical integration method.\n";
00332
00333 throw std::logic_error(buffer.str());
00334 }
00335 break;
00336 }
00337
00338 }
00339
00340
00341
00343
00344 TiXmlElement* NumericalIntegration::to_XML(void) const
00345 {
00346 std::ostringstream buffer;
00347
00348
00349
00350 TiXmlElement* numerical_integration_element = new TiXmlElement("NumericalIntegration");
00351 numerical_integration_element->SetAttribute("Version", 4);
00352
00353
00354 {
00355 TiXmlElement* element = new TiXmlElement("NumericalIntegrationMethod");
00356 numerical_integration_element->LinkEndChild(element);
00357
00358 TiXmlText* text = new TiXmlText(write_numerical_integration_method().c_str());
00359 element->LinkEndChild(text);
00360 }
00361
00362
00363 {
00364 TiXmlElement* element = new TiXmlElement("Display");
00365 numerical_integration_element->LinkEndChild(element);
00366
00367 buffer.str("");
00368 buffer << display;
00369
00370 TiXmlText* text = new TiXmlText(buffer.str().c_str());
00371 element->LinkEndChild(text);
00372 }
00373
00374 return(numerical_integration_element);
00375 }
00376
00377
00378
00379
00382
00383 void NumericalIntegration::from_XML(TiXmlElement* numerical_integration_element)
00384 {
00385 if(numerical_integration_element)
00386 {
00387
00388 {
00389 TiXmlElement* element = numerical_integration_element->FirstChildElement("NumericalIntegrationMethod");
00390
00391 if(element)
00392 {
00393 std::string new_numerical_integration_method = element->GetText();
00394
00395 try
00396 {
00397 set_numerical_integration_method(new_numerical_integration_method);
00398 }
00399 catch(std::exception& e)
00400 {
00401 std::cout << e.what() << std::endl;
00402 }
00403 }
00404 }
00405
00406
00407 {
00408 TiXmlElement* element = numerical_integration_element->FirstChildElement("Display");
00409
00410 if(element)
00411 {
00412 std::string new_display = element->GetText();
00413
00414 try
00415 {
00416 set_display(new_display != "0");
00417 }
00418 catch(std::exception& e)
00419 {
00420 std::cout << e.what() << std::endl;
00421 }
00422 }
00423 }
00424 }
00425 }
00426
00427
00428 }
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446