C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
Simple programming examples

On this page we will try to show you, how to use the C-XSC class library with some short and simple examples.

Example 1 - An Introduction

We will start with the following simple program showing you how to use intervals, compute one operation and print out the result:

1 #include "interval.hpp" // predefined interval arithmetic
2 #include <iostream>
3 using namespace cxsc;
4 using namespace std;
5 
6 int main()
7 {
8  interval a, b; // Standard intervals
9  a = 1.0; // a = [1.0,1.0]
10  b = 3.0; // b = [3.0,3.0]
11  cout << "a/b = " << a/b << endl;
12  return 0;
13 }
14 
15 /* --------------------------- Output ------------------------------
16 a/b = [ 0.333333, 0.333334]
17 ------------------------------------------------------------------*/

Let's start investigating the interesting lines. With the line

0 #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.

1 using namespace cxsc;

Then you declare your variables and assign adequate values in the following lines.

3  interval a, b; // Standard intervals
8  a = 1.0; // a = [1.0,1.0]
9  b = 3.0; // b = [3.0,3.0]

Finally you want to print out the result for your desired computation.

10  cout << "a/b = " << a/b << endl;

So it's just that easy to use C-XSC.

Example 2 - Input / Output

1 #include <iostream>
2 #include "interval.hpp"
3 using namespace cxsc;
4 using namespace std;
5 
6 int main()
7 {
8  real a, b;
9 
10  cout << "Please enter real a: ";
11  cout << RndDown; // set rounding mode
12  cin >> a; // read a rounded downwards
13  cout << SetPrecision(7,4); // set field width and number of digits
14  // for output
15  cout << a << endl; // print a rounded downwards to 4 digits
16  cout << RndUp;
17  cout << a << endl; // print a rounded upwards to 4 digits
18 
19  "0.3" >> b; // convert the string "0.3" to a floating
20  // point value b using rounding down
21 
22  cout << SetPrecision(18,15); // from now on print 15 digits
23 
24  cout << b << endl; // print b rounded upwards to 15 digits
25  cout << RndDown;
26  cout << b << endl; // print b rounded downwards to 15 digits
27  interval x;
28  "[1.5, 2]" >> x; // string to interval conversion
29  cout << x << endl; // print interval x using default setting
30 
31  cout << SaveOpt; // push I/O parameters to internal stack
32  cout << SetPrecision(10,7); // modifies output field width and
33  // number of digits to be print
34  cout << x << endl; // print x in the modified format
35  cout << RestoreOpt; // pop parameters from internal stack
36  cout << x << endl; // again, print x using the former format
37  return 0;
38 }
39 
40 /* --------------------------- Output ------------------------------
41 Please enter real a: 0.3
42  0.2999
43  0.3000
44  0.300000000000001
45  0.300000000000000
46 [ 1.500000000000000, 2.000000000000000]
47 [ 1.5000000, 2.0000000]
48 [ 1.500000000000000, 2.000000000000000]
49 ------------------------------------------------------------------*/

Example 3 - Compute all zeros of a function

1 // Compute all zeros of the function
2 //
3 // (x-1)*(exp(-3*x) - power(sin(x), 3))
4 //
5 // This example needs the CToolbox sources additionally!
6 
7 #include "nlfzero.hpp" // Nonlinear equations module
8 #include "stacksz.hpp" // To increase stack size for some
9  // special C++ compiler
10 using namespace cxsc;
11 using namespace std;
12 
13 
14 DerivType f ( const DerivType& x ) // Sample function
15 {
16  return (x-1)*( exp(-3*x) - power(sin(x),3) );
17 }
18 
19 // The class DerivType allows the computation of f, f', and f'' using
20 // automatic differentiation; see "C++ Toolbox for Verified Computing"
21 
22 int main()
23 {
24  interval SearchInterval;
25  real Tolerance;
26  ivector Zero;
27  intvector Unique;
28  int NumberOfZeros, i, Error;
29 
30  cout << SetPrecision(23,15) << Scientific; // Output format
31 
32  cout << "Search interval : ";
33  cin >> SearchInterval;
34  cout << "Tolerance (relative): ";
35  cin >> Tolerance;
36  cout << endl;
37 
38  // Call the function 'AllZeros()' from the C++ Toolbox
39  AllZeros(f,SearchInterval,Tolerance,Zero,Unique,NumberOfZeros,Error);
40 
41  for ( i = 1; i <= NumberOfZeros; i++) {
42  cout << Zero[i] << endl;
43  if (Unique[i])
44  cout << "encloses a locally unique zero!" << endl;
45  else
46  cout << "may contain a zero (not verified unique)!" << endl;
47  }
48  cout << endl << NumberOfZeros << " interval enclosure(s)" << endl;
49  if (Error) cout << endl << AllZerosErrMsg(Error) << endl;
50  return 0;
51 }
52 
53 /* --------------------------- Output ------------------------------
54 Search interval : [-1,15]
55 Tolerance (relative) : 1e-10
56 
57 [ 5.885327439818601E-001, 5.885327439818619E-001]
58 encloses a locally unique zero!
59 [ 9.999999999999998E-001, 1.000000000000001E+000]
60 encloses a locally unique zero!
61 [ 3.096363932404308E+000, 3.096363932416931E+000]
62 encloses a locally unique zero!
63 [ 6.285049273371415E+000, 6.285049273396501E+000]
64 encloses a locally unique zero!
65 [ 9.424697254738511E+000, 9.424697254738533E+000]
66 encloses a locally unique zero!
67 [ 1.256637410119757E+001, 1.256637410231546E+001]
68 encloses a locally unique zero!
69 
70 6 interval enclosure(s)
71 ------------------------------------------------------------------*/

Example 4 - Interval Newton method

1 // Interval Newton method using ordinary interval arithmetic
2 // Verified computation of a zero of the function f()
3 
4 #include <iostream>
5 #include "interval.hpp" // Include interval arithmetic package
6 #include "imath.hpp" // Include interval standard functions
7 using namespace std;
8 using namespace cxsc;
9 
10 interval f(const real x)
11 { // Function f
12  interval y(x); // y is a thin interval initialized by x
13  return sqrt(y) + (y+1)*cos(y); // Interval arithmetic is used
14 }
15 
16 interval deriv(const interval& x)
17 { // Derivative function f'
18  return 1/(2*sqrt(x)) + cos(x) - (x+1)*sin(x);
19 }
20 
21 bool criter(const interval& x) // Computing: f(a)*f(b) < 0 and
22 { // not 0 in f'([x])?
23  return Sup( f(Inf(x))*f(Sup(x)) ) < 0.0 && !(0.0 <= deriv(x));
24 } // '<=' means 'element of'
25 
26 int main(void)
27 {
28  interval x, xOld;
29  cout << SetPrecision(20,15); // Number of mantissa digits in I/O
30  x= interval(2,3);
31  cout << "Starting interval is [2,3]" << endl;
32  if (criter(x))
33  { // There is exactly one zero of f in the interval x
34  do {
35  xOld = x;
36  cout << "Actual enclosure is " << x << endl;
37  x = (mid(x)-f(mid(x))/deriv(x)) & x; // Iteration formula
38  } while (x != xOld);
39  cout << "Final enclosure of the zero: " << x << endl;
40  }
41  else
42  cout << "Criterion not satisfied!" << endl;
43  return 0;
44 }
45 
46 /* Output:
47 
48 Starting interval is [2,3]
49 Actual enclosure is [ 2.000000000000000, 3.000000000000000]
50 Actual enclosure is [ 2.000000000000000, 2.218137182953809]
51 Actual enclosure is [ 2.051401462380920, 2.064726329907714]
52 Actual enclosure is [ 2.059037791936965, 2.059053011233253]
53 Actual enclosure is [ 2.059045253413831, 2.059045253416460]
54 Actual enclosure is [ 2.059045253415142, 2.059045253415145]
55 Final enclosure of the zero: [ 2.059045253415142,
56  2.059045253415145]
57 
58 */

Example 5 -

1 #include "l_interval.hpp" // interval staggered arithmetic
2 #include <iostream>
3 using namespace cxsc;
4 using namespace std;
5 
6 int main()
7 {
8  l_interval a, b; // Multiple-precision intervals
9  a = 1.0;
10  b = 3.0;
11  stagprec = 2; // global integer variable
12  cout << SetDotPrecision(16*stagprec, 16*stagprec-3) << RndNext;
13  // I/O for variables of type l_interval is done using
14  // the long accumulator (i.e. a dotprecision variable)
15 
16  cout << "a/b = " << a/b << endl;
17  return(0);
18 }
19 
20 /* --------------------------- Output ------------------------------
21 a/b = [ 0.33333333333333333333333333333, 0.33333333333333333333333333334]
22 ------------------------------------------------------------------*/

Example 6 -

1 // Interval Newton method using a multi-precision interval data type
2 // Verified computation of a zero of the function f() to high accuracy
3 
4 #include <iostream>
5 #include "l_interval.hpp" // Include multi-precision intervals
6 #include "l_imath.hpp" // Include multi-precision math functions
7 #include "l_rmath.hpp"
8 using namespace std;
9 using namespace cxsc;
10 
11 l_interval f(const l_real& x) // Function f
12 {
13  l_interval y(x); // y is a thin interval initialized by x
14  return sqrt(y) + (y+1)*cos(y); // Use multi-precision interval arithmetic
15 }
16 
17 l_interval deriv(const l_interval& x) // Derivative function f'
18 {
19  return 1/(2*sqrt(x)) + cos(x) - (x+1)*sin(x);
20 }
21 
22 bool criter(const l_interval& x) // Computing: f(a)*f(b) < 0 and
23 { // not 0 in f'([x])?
24  return Sup( f(Inf(x))*f(Sup(x)) ) < 0.0 && !(0.0 <= deriv(x));
25 } // '<=' means 'element of'
26 
27 int main(void)
28 {
29  l_interval x, xOld;
30  stagprec= 3; // Set precision of the staggered correction arithmetic
31  x= l_interval(2,3);
32  cout << "Starting interval is [2,3]" << endl;
33  cout << SetDotPrecision(16*stagprec, 16*stagprec-3) << RndNext;
34  // I/O for variables of type l_interval is done using
35  // the long accumulator (i.e. a dotprecision variable)
36 
37  if (criter(x))
38  { // There is exactly one zero of f in the interval x
39  do {
40  xOld = x;
41  cout << "Diameter of actual enclosure: " << real(diam(x)) << endl;
42  x = (mid(x)-f(mid(x))/deriv(x)) & x; // Iteration formula
43  } while (x != xOld); // & means intersection
44  cout << "Final enclosure of the zero: " << x << endl;
45  }
46  else
47  cout << "Criterion not satisfied!" << endl;
48  return 0;
49 }
50 
51 /* Output:
52 
53 Starting interval is [2,3])
54 Diameter of actual enclosure: 1.000000
55 Diameter of actual enclosure: 0.218137
56 Diameter of actual enclosure: 0.013325
57 Diameter of actual enclosure: 1.521930E-005
58 Diameter of actual enclosure: 2.625899E-012
59 Diameter of actual enclosure: 4.708711E-027
60 Diameter of actual enclosure: 2.138212E-049
61 Final enclosure of the zero:
62  [ 2.059045253415143788680636155343254522623083897,
63  2.059045253415143788680636155343254522623083898 ]
64 */
65 
66 

Example 7 -

1 // Runge-Kutta Method
2 
3 #include <iostream>
4 #include "rvector.hpp" // Include dynamic arrays (real vectors)
5 using namespace std;
6 using namespace cxsc;
7 
8 rvector F(const real x, const rvector& Y) // Function definition
9 {
10  rvector Z(3); // Constructor call
11 
12  Z[1] = Y[2] * Y[3]; // F is independent of x
13  Z[2] = -Y[1] * Y[3];
14  Z[3] = -0.522 * Y[1] * Y[2];
15  return Z;
16 }
17 
18 void Init(real& x, real& h, rvector& Y)
19 { // Initialization
20  x = 0.0; h = 0.1;
21  Y[1] = 0.0; Y[2] = 1.0; Y[3] = 1.0;
22 }
23 
24 int main(void)
25 {
26  real x, h; // Declarations and dynamic
27  rvector Y(3), K1, K2, K3, K4; // memory allocation
28  // Automatic resize of Ki below
29  Init(x, h, Y);
30  for (int i=1; i<=3; i++) { // 3 Runge-Kutta steps
31  K1 = h * F(x, Y); // with array result
32  K2 = h * F(x + h / 2, Y + K1 / 2);
33  K3 = h * F(x + h / 2, Y + K2 / 2);
34  K4 = h * F(x + h, Y + K3);
35  Y = Y + (K1 + 2 * K2 + 2 * K3 + K4) / 6;
36  x += h;
37  cout << SetPrecision(10,6) << Dec; // I/O modification
38  cout << "Step: " << i << ", "
39  << "x = " << x << endl;
40  cout << "Y = " << endl << Y << endl;
41  }
42  return 0;
43 }
44 
45 /* --------------------------- Output ------------------------------
46 Step: 1, x = 0.100000
47 Y =
48  0.099747
49  0.995013
50  0.997400
51 
52 Step: 2, x = 0.200000
53 Y =
54  0.197993
55  0.980203
56  0.989716
57 
58 Step: 3, x = 0.300000
59 Y =
60  0.293320
61  0.956014
62  0.977286
63 ------------------------------------------------------------------*/

Example 8 -

1 // Trace of a (complex) matrix product
2 // Let C denote the matrix product A*B.
3 // Then the diagonal entries of C are added to get the trace of C.
4 
5 #include <iostream>
6 #include "cmatrix.hpp" // Use the complex matrix package
7 using namespace std;
8 using namespace cxsc;
9 
10 int main()
11 {
12  int n;
13  cout << "Please enter the matrix dimension n: "; cin >> n;
14  cmatrix A(n,n), B(n,n); // Dynamic allocation of A, B
15  cdotprecision accu; // Allows exact computation of dotproducts
16  cout << "Please enter the matrix A:" << endl; cin >> A;
17  cout << "Please enter the matrix B:" << endl; cin >> B;
18  accu = 0.0; // Clear accumulator
19  for (int i=1; i<=n; i++) accumulate(accu, A[i], B[Col(i)]);
20  // A[i] and B[Col(i)] are subarrays of type rvector.
21 
22  // The exact result stored in the complex dotprecision variable accu
23  // is rounded to the nearest complex floating point number:
24  complex result = rnd(accu);
25  cout << SetPrecision(12,6) << RndNext << Dec;
26  cout << "Trace of product matrix: " << result << endl;
27  return 0;
28 }
29 
30 /* --------------------------- Output ------------------------------
31 Please enter the matrix dimension n: 3
32 Please enter the matrix A:
33 (1,0) (2,0) (3,0)
34 (4,0) (5,0) (6,0)
35 (7,0) (8,0) (9,0)
36 Please enter the matrix B:
37 (10,0) (11,0) (12,0)
38 (13,0) (14,0) (15,0)
39 (16,0) (17,0) (18,0)
40 Trace of product matrix: ( 666.000000, 0.000000)
41 ------------------------------------------------------------------*/
cxsc::cmatrix
The Data Type cmatrix.
Definition: cmatrix.hpp:513
cxsc::mid
cvector mid(const cimatrix_subv &mv)
Returns the middle of the matrix.
Definition: cimatrix.inl:739
cxsc::interval
The Scalar Type interval.
Definition: interval.hpp:54
cxsc::intvector
The Data Type intvector.
Definition: intvector.hpp:51
cxsc::Zero
int Zero(ivector &x)
Checks if vector is zero vector.
Definition: ivector.cpp:80
cxsc::diam
cvector diam(const cimatrix_subv &mv)
Returns the diameter of the matrix.
Definition: cimatrix.inl:738
cxsc::rvector
The Data Type rvector.
Definition: rvector.hpp:57
cxsc::sqrt
cinterval sqrt(const cinterval &z)
Calculates .
Definition: cimath.cpp:1007
cxsc::Col
cimatrix_subv Col(cimatrix &m, const int &i)
Returns one column of the matrix as a vector.
Definition: cimatrix.inl:242
cxsc::l_interval
The Multiple-Precision Data Type l_interval.
Definition: l_interval.hpp:71
cxsc::ivector
The Data Type ivector.
Definition: ivector.hpp:54
cxsc::power
cinterval power(const cinterval &z, int n)
Calculates .
Definition: cimath.cpp:1941
cxsc::cos
cinterval cos(const cinterval &z)
Calculates .
Definition: cimath.cpp:207
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::l_real
The Multiple-Precision Data Type l_real.
Definition: l_real.hpp:77
cxsc::sin
cinterval sin(const cinterval &z)
Calculates .
Definition: cimath.cpp:215
cxsc::cdotprecision
The Data Type cdotprecision.
Definition: cdot.hpp:60
cxsc::exp
cinterval exp(const cinterval &z)
Calculates .
Definition: cimath.cpp:159
cxsc::complex
The Scalar Type complex.
Definition: complex.hpp:49
cxsc::real
The Scalar Type real.
Definition: real.hpp:113