tmainD.cpp - numeric - C++ library with numerical algorithms
(HTM) git clone git://src.adamsgaard.dk/numeric
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tmainD.cpp (1871B)
---
1 #include <iostream>
2 #include <vector>
3 #include <complex>
4 #include <cmath>
5 #include <fstream>
6 #include "typedefs.h"
7 #include "check.h"
8 #include "ode.h"
9 #include "functions.h"
10
11
12 int main()
13 {
14 // Namespace declarations
15 using std::cout;
16 using std::vector;
17 using std::complex;
18
19 // Calculate machine precision
20 Floattype eps_machine = 1.0f;
21 while (1.0f + eps_machine != 1.0f)
22 eps_machine /= 2.0f;
23
24 cout << "\n\033[1;33m--- Part D: Precision analysis ---\033[0m\n";
25 complex<Floattype> a(0.0f, 0.0f); // Lower limit
26 complex<Floattype> b(2.0f*M_PI, 2.0f*M_PI); // Upper limit
27 cout << "Integration path: b-a = " << b-a << '\n';
28 Inttype n_eqs = 2; // Number of equations in ODE system
29 vector<complex<Floattype> > y_start(n_eqs);
30 complex<Floattype> y0(0.0f, 0.0f);
31 complex<Floattype> y1(1.0f, 1.0f);
32 y_start[0] = y0;
33 y_start[1] = y1;
34 Floattype h_start = 0.01f;
35
36 vector<Floattype> precs; // Vector containing precision values
37 vector<Inttype> steps; // Vector containing number of steps required
38
39 for (Floattype prec=eps_machine*10.0f; prec<0.1f; prec *= 10.0f) {
40 ODE ode(func1, // ODE system
41 y_start, // Initial values
42 a, // Lower limit
43 b, // Upper limit
44 h_start, // Start value of step size
45 100000, // Max. number of steps
46 prec, // Absolute precision
47 prec); // Relative precision
48 precs.push_back(prec); // Save precision
49 steps.push_back(ode.steps()); // Save number of steps taken
50 }
51
52 // Save results to text file
53 std::ofstream ost; // Out stream object
54 ost.open("funcD.dat"); // Open outfile for write
55 if (!ost) {
56 std::cerr << "Error, can't open output file!\n";
57 return 1;
58 }
59 // Write output values
60 for (Inttype i=0; i<precs.size(); ++i)
61 ost << precs[i] << '\t' << steps[i] << '\n';
62
63 // Close file
64 ost.close();
65
66 // Return successfully
67 return 0;
68 }