C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
On this page we will try to show you, how to use the C-XSC class library with some short and simple examples.
We will start with the following simple program showing you how to use intervals, compute one operation and print out the result:
00001 #include "interval.hpp" // predefined interval arithmetic 00002 #include <iostream> 00003 using namespace cxsc; 00004 using namespace std; 00005 00006 int main() 00007 { 00008 interval a, b; // Standard intervals 00009 a = 1.0; // a = [1.0,1.0] 00010 b = 3.0; // b = [3.0,3.0] 00011 cout << "a/b = " << a/b << endl; 00012 return 0; 00013 } 00014 00015 /* --------------------------- Output ------------------------------ 00016 a/b = [ 0.333333, 0.333334] 00017 ------------------------------------------------------------------*/
Let's start investigating the interesting lines. With the line
#include "interval.hpp" // predefined interval arithmetic
you include the functionality of the interval class of C-XSC in your program. After that you have to inform the compiler about the namespace cxsc - where all classes and methods of C-XSC are stored in - to use C-XSC without fully qualified identifiers.
using namespace cxsc;
Then you declare your variables and assign adequate values in the following lines.
interval a, b; // Standard intervals a = 1.0; // a = [1.0,1.0] b = 3.0; // b = [3.0,3.0]
Finally you want to print out the result for your desired computation.
cout << "a/b = " << a/b << endl;
So it's just that easy to use C-XSC.
00001 #include <iostream> 00002 #include "interval.hpp" 00003 using namespace cxsc; 00004 using namespace std; 00005 00006 int main() 00007 { 00008 real a, b; 00009 00010 cout << "Please enter real a: "; 00011 cout << RndDown; // set rounding mode 00012 cin >> a; // read a rounded downwards 00013 cout << SetPrecision(7,4); // set field width and number of digits 00014 // for output 00015 cout << a << endl; // print a rounded downwards to 4 digits 00016 cout << RndUp; 00017 cout << a << endl; // print a rounded upwards to 4 digits 00018 00019 "0.3" >> b; // convert the string "0.3" to a floating 00020 // point value b using rounding down 00021 00022 cout << SetPrecision(18,15); // from now on print 15 digits 00023 00024 cout << b << endl; // print b rounded upwards to 15 digits 00025 cout << RndDown; 00026 cout << b << endl; // print b rounded downwards to 15 digits 00027 interval x; 00028 "[1.5, 2]" >> x; // string to interval conversion 00029 cout << x << endl; // print interval x using default setting 00030 00031 cout << SaveOpt; // push I/O parameters to internal stack 00032 cout << SetPrecision(10,7); // modifies output field width and 00033 // number of digits to be print 00034 cout << x << endl; // print x in the modified format 00035 cout << RestoreOpt; // pop parameters from internal stack 00036 cout << x << endl; // again, print x using the former format 00037 return 0; 00038 } 00039 00040 /* --------------------------- Output ------------------------------ 00041 Please enter real a: 0.3 00042 0.2999 00043 0.3000 00044 0.300000000000001 00045 0.300000000000000 00046 [ 1.500000000000000, 2.000000000000000] 00047 [ 1.5000000, 2.0000000] 00048 [ 1.500000000000000, 2.000000000000000] 00049 ------------------------------------------------------------------*/
00001 // Compute all zeros of the function 00002 // 00003 // (x-1)*(exp(-3*x) - power(sin(x), 3)) 00004 // 00005 // This example needs the CToolbox sources additionally! 00006 00007 #include "nlfzero.hpp" // Nonlinear equations module 00008 #include "stacksz.hpp" // To increase stack size for some 00009 // special C++ compiler 00010 using namespace cxsc; 00011 using namespace std; 00012 00013 00014 DerivType f ( const DerivType& x ) // Sample function 00015 { 00016 return (x-1)*( exp(-3*x) - power(sin(x),3) ); 00017 } 00018 00019 // The class DerivType allows the computation of f, f', and f'' using 00020 // automatic differentiation; see "C++ Toolbox for Verified Computing" 00021 00022 int main() 00023 { 00024 interval SearchInterval; 00025 real Tolerance; 00026 ivector Zero; 00027 intvector Unique; 00028 int NumberOfZeros, i, Error; 00029 00030 cout << SetPrecision(23,15) << Scientific; // Output format 00031 00032 cout << "Search interval : "; 00033 cin >> SearchInterval; 00034 cout << "Tolerance (relative): "; 00035 cin >> Tolerance; 00036 cout << endl; 00037 00038 // Call the function 'AllZeros()' from the C++ Toolbox 00039 AllZeros(f,SearchInterval,Tolerance,Zero,Unique,NumberOfZeros,Error); 00040 00041 for ( i = 1; i <= NumberOfZeros; i++) { 00042 cout << Zero[i] << endl; 00043 if (Unique[i]) 00044 cout << "encloses a locally unique zero!" << endl; 00045 else 00046 cout << "may contain a zero (not verified unique)!" << endl; 00047 } 00048 cout << endl << NumberOfZeros << " interval enclosure(s)" << endl; 00049 if (Error) cout << endl << AllZerosErrMsg(Error) << endl; 00050 return 0; 00051 } 00052 00053 /* --------------------------- Output ------------------------------ 00054 Search interval : [-1,15] 00055 Tolerance (relative) : 1e-10 00056 00057 [ 5.885327439818601E-001, 5.885327439818619E-001] 00058 encloses a locally unique zero! 00059 [ 9.999999999999998E-001, 1.000000000000001E+000] 00060 encloses a locally unique zero! 00061 [ 3.096363932404308E+000, 3.096363932416931E+000] 00062 encloses a locally unique zero! 00063 [ 6.285049273371415E+000, 6.285049273396501E+000] 00064 encloses a locally unique zero! 00065 [ 9.424697254738511E+000, 9.424697254738533E+000] 00066 encloses a locally unique zero! 00067 [ 1.256637410119757E+001, 1.256637410231546E+001] 00068 encloses a locally unique zero! 00069 00070 6 interval enclosure(s) 00071 ------------------------------------------------------------------*/
00001 // Interval Newton method using ordinary interval arithmetic 00002 // Verified computation of a zero of the function f() 00003 00004 #include <iostream> 00005 #include "interval.hpp" // Include interval arithmetic package 00006 #include "imath.hpp" // Include interval standard functions 00007 using namespace std; 00008 using namespace cxsc; 00009 00010 interval f(const real x) 00011 { // Function f 00012 interval y(x); // y is a thin interval initialized by x 00013 return sqrt(y) + (y+1)*cos(y); // Interval arithmetic is used 00014 } 00015 00016 interval deriv(const interval& x) 00017 { // Derivative function f' 00018 return 1/(2*sqrt(x)) + cos(x) - (x+1)*sin(x); 00019 } 00020 00021 bool criter(const interval& x) // Computing: f(a)*f(b) < 0 and 00022 { // not 0 in f'([x])? 00023 return Sup( f(Inf(x))*f(Sup(x)) ) < 0.0 && !(0.0 <= deriv(x)); 00024 } // '<=' means 'element of' 00025 00026 int main(void) 00027 { 00028 interval x, xOld; 00029 cout << SetPrecision(20,15); // Number of mantissa digits in I/O 00030 x= interval(2,3); 00031 cout << "Starting interval is [2,3]" << endl; 00032 if (criter(x)) 00033 { // There is exactly one zero of f in the interval x 00034 do { 00035 xOld = x; 00036 cout << "Actual enclosure is " << x << endl; 00037 x = (mid(x)-f(mid(x))/deriv(x)) & x; // Iteration formula 00038 } while (x != xOld); 00039 cout << "Final enclosure of the zero: " << x << endl; 00040 } 00041 else 00042 cout << "Criterion not satisfied!" << endl; 00043 return 0; 00044 } 00045 00046 /* Output: 00047 00048 Starting interval is [2,3] 00049 Actual enclosure is [ 2.000000000000000, 3.000000000000000] 00050 Actual enclosure is [ 2.000000000000000, 2.218137182953809] 00051 Actual enclosure is [ 2.051401462380920, 2.064726329907714] 00052 Actual enclosure is [ 2.059037791936965, 2.059053011233253] 00053 Actual enclosure is [ 2.059045253413831, 2.059045253416460] 00054 Actual enclosure is [ 2.059045253415142, 2.059045253415145] 00055 Final enclosure of the zero: [ 2.059045253415142, 00056 2.059045253415145] 00057 00058 */
00001 #include "l_interval.hpp" // interval staggered arithmetic 00002 #include <iostream> 00003 using namespace cxsc; 00004 using namespace std; 00005 00006 int main() 00007 { 00008 l_interval a, b; // Multiple-precision intervals 00009 a = 1.0; 00010 b = 3.0; 00011 stagprec = 2; // global integer variable 00012 cout << SetDotPrecision(16*stagprec, 16*stagprec-3) << RndNext; 00013 // I/O for variables of type l_interval is done using 00014 // the long accumulator (i.e. a dotprecision variable) 00015 00016 cout << "a/b = " << a/b << endl; 00017 return(0); 00018 } 00019 00020 /* --------------------------- Output ------------------------------ 00021 a/b = [ 0.33333333333333333333333333333, 0.33333333333333333333333333334] 00022 ------------------------------------------------------------------*/
00001 // Interval Newton method using a multi-precision interval data type 00002 // Verified computation of a zero of the function f() to high accuracy 00003 00004 #include <iostream> 00005 #include "l_interval.hpp" // Include multi-precision intervals 00006 #include "l_imath.hpp" // Include multi-precision math functions 00007 #include "l_rmath.hpp" 00008 using namespace std; 00009 using namespace cxsc; 00010 00011 l_interval f(const l_real& x) // Function f 00012 { 00013 l_interval y(x); // y is a thin interval initialized by x 00014 return sqrt(y) + (y+1)*cos(y); // Use multi-precision interval arithmetic 00015 } 00016 00017 l_interval deriv(const l_interval& x) // Derivative function f' 00018 { 00019 return 1/(2*sqrt(x)) + cos(x) - (x+1)*sin(x); 00020 } 00021 00022 bool criter(const l_interval& x) // Computing: f(a)*f(b) < 0 and 00023 { // not 0 in f'([x])? 00024 return Sup( f(Inf(x))*f(Sup(x)) ) < 0.0 && !(0.0 <= deriv(x)); 00025 } // '<=' means 'element of' 00026 00027 int main(void) 00028 { 00029 l_interval x, xOld; 00030 stagprec= 3; // Set precision of the staggered correction arithmetic 00031 x= l_interval(2,3); 00032 cout << "Starting interval is [2,3]" << endl; 00033 cout << SetDotPrecision(16*stagprec, 16*stagprec-3) << RndNext; 00034 // I/O for variables of type l_interval is done using 00035 // the long accumulator (i.e. a dotprecision variable) 00036 00037 if (criter(x)) 00038 { // There is exactly one zero of f in the interval x 00039 do { 00040 xOld = x; 00041 cout << "Diameter of actual enclosure: " << real(diam(x)) << endl; 00042 x = (mid(x)-f(mid(x))/deriv(x)) & x; // Iteration formula 00043 } while (x != xOld); // & means intersection 00044 cout << "Final enclosure of the zero: " << x << endl; 00045 } 00046 else 00047 cout << "Criterion not satisfied!" << endl; 00048 return 0; 00049 } 00050 00051 /* Output: 00052 00053 Starting interval is [2,3]) 00054 Diameter of actual enclosure: 1.000000 00055 Diameter of actual enclosure: 0.218137 00056 Diameter of actual enclosure: 0.013325 00057 Diameter of actual enclosure: 1.521930E-005 00058 Diameter of actual enclosure: 2.625899E-012 00059 Diameter of actual enclosure: 4.708711E-027 00060 Diameter of actual enclosure: 2.138212E-049 00061 Final enclosure of the zero: 00062 [ 2.059045253415143788680636155343254522623083897, 00063 2.059045253415143788680636155343254522623083898 ] 00064 */ 00065 00066
00001 // Runge-Kutta Method 00002 00003 #include <iostream> 00004 #include "rvector.hpp" // Include dynamic arrays (real vectors) 00005 using namespace std; 00006 using namespace cxsc; 00007 00008 rvector F(const real x, const rvector& Y) // Function definition 00009 { 00010 rvector Z(3); // Constructor call 00011 00012 Z[1] = Y[2] * Y[3]; // F is independent of x 00013 Z[2] = -Y[1] * Y[3]; 00014 Z[3] = -0.522 * Y[1] * Y[2]; 00015 return Z; 00016 } 00017 00018 void Init(real& x, real& h, rvector& Y) 00019 { // Initialization 00020 x = 0.0; h = 0.1; 00021 Y[1] = 0.0; Y[2] = 1.0; Y[3] = 1.0; 00022 } 00023 00024 int main(void) 00025 { 00026 real x, h; // Declarations and dynamic 00027 rvector Y(3), K1, K2, K3, K4; // memory allocation 00028 // Automatic resize of Ki below 00029 Init(x, h, Y); 00030 for (int i=1; i<=3; i++) { // 3 Runge-Kutta steps 00031 K1 = h * F(x, Y); // with array result 00032 K2 = h * F(x + h / 2, Y + K1 / 2); 00033 K3 = h * F(x + h / 2, Y + K2 / 2); 00034 K4 = h * F(x + h, Y + K3); 00035 Y = Y + (K1 + 2 * K2 + 2 * K3 + K4) / 6; 00036 x += h; 00037 cout << SetPrecision(10,6) << Dec; // I/O modification 00038 cout << "Step: " << i << ", " 00039 << "x = " << x << endl; 00040 cout << "Y = " << endl << Y << endl; 00041 } 00042 return 0; 00043 } 00044 00045 /* --------------------------- Output ------------------------------ 00046 Step: 1, x = 0.100000 00047 Y = 00048 0.099747 00049 0.995013 00050 0.997400 00051 00052 Step: 2, x = 0.200000 00053 Y = 00054 0.197993 00055 0.980203 00056 0.989716 00057 00058 Step: 3, x = 0.300000 00059 Y = 00060 0.293320 00061 0.956014 00062 0.977286 00063 ------------------------------------------------------------------*/
00001 // Trace of a (complex) matrix product 00002 // Let C denote the matrix product A*B. 00003 // Then the diagonal entries of C are added to get the trace of C. 00004 00005 #include <iostream> 00006 #include "cmatrix.hpp" // Use the complex matrix package 00007 using namespace std; 00008 using namespace cxsc; 00009 00010 int main() 00011 { 00012 int n; 00013 cout << "Please enter the matrix dimension n: "; cin >> n; 00014 cmatrix A(n,n), B(n,n); // Dynamic allocation of A, B 00015 cdotprecision accu; // Allows exact computation of dotproducts 00016 cout << "Please enter the matrix A:" << endl; cin >> A; 00017 cout << "Please enter the matrix B:" << endl; cin >> B; 00018 accu = 0.0; // Clear accumulator 00019 for (int i=1; i<=n; i++) accumulate(accu, A[i], B[Col(i)]); 00020 // A[i] and B[Col(i)] are subarrays of type rvector. 00021 00022 // The exact result stored in the complex dotprecision variable accu 00023 // is rounded to the nearest complex floating point number: 00024 complex result = rnd(accu); 00025 cout << SetPrecision(12,6) << RndNext << Dec; 00026 cout << "Trace of product matrix: " << result << endl; 00027 return 0; 00028 } 00029 00030 /* --------------------------- Output ------------------------------ 00031 Please enter the matrix dimension n: 3 00032 Please enter the matrix A: 00033 (1,0) (2,0) (3,0) 00034 (4,0) (5,0) (6,0) 00035 (7,0) (8,0) (9,0) 00036 Please enter the matrix B: 00037 (10,0) (11,0) (12,0) 00038 (13,0) (14,0) (15,0) 00039 (16,0) (17,0) (18,0) 00040 Trace of product matrix: ( 666.000000, 0.000000) 00041 ------------------------------------------------------------------*/