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>
3using namespace cxsc;
4using namespace std;
5
6int 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 ------------------------------
16a/b = [ 0.333333, 0.333334]
17------------------------------------------------------------------*/
The Scalar Type interval.
Definition interval.hpp:55
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition cdot.cpp:29

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.

1using 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"
3using namespace cxsc;
4using namespace std;
5
6int 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 ------------------------------
41Please 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------------------------------------------------------------------*/
The Scalar Type real.
Definition real.hpp:114

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
10using namespace cxsc;
11using namespace std;
12
13
14DerivType 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
22int 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 ------------------------------
54Search interval : [-1,15]
55Tolerance (relative) : 1e-10
56
57[ 5.885327439818601E-001, 5.885327439818619E-001]
58encloses a locally unique zero!
59[ 9.999999999999998E-001, 1.000000000000001E+000]
60encloses a locally unique zero!
61[ 3.096363932404308E+000, 3.096363932416931E+000]
62encloses a locally unique zero!
63[ 6.285049273371415E+000, 6.285049273396501E+000]
64encloses a locally unique zero!
65[ 9.424697254738511E+000, 9.424697254738533E+000]
66encloses a locally unique zero!
67[ 1.256637410119757E+001, 1.256637410231546E+001]
68encloses a locally unique zero!
69
706 interval enclosure(s)
71------------------------------------------------------------------*/
The Data Type intvector.
Definition intvector.hpp:52
The Data Type ivector.
Definition ivector.hpp:55

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
7using namespace std;
8using namespace cxsc;
9
10interval 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
16interval deriv(const interval& x)
17{ // Derivative function f'
18 return 1/(2*sqrt(x)) + cos(x) - (x+1)*sin(x);
19}
20
21bool 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
26int 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
48Starting interval is [2,3]
49Actual enclosure is [ 2.000000000000000, 3.000000000000000]
50Actual enclosure is [ 2.000000000000000, 2.218137182953809]
51Actual enclosure is [ 2.051401462380920, 2.064726329907714]
52Actual enclosure is [ 2.059037791936965, 2.059053011233253]
53Actual enclosure is [ 2.059045253413831, 2.059045253416460]
54Actual enclosure is [ 2.059045253415142, 2.059045253415145]
55Final enclosure of the zero: [ 2.059045253415142,
56 2.059045253415145]
57
58*/
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition cimatrix.inl:739

Example 5 -

1#include "l_interval.hpp" // interval staggered arithmetic
2#include <iostream>
3using namespace cxsc;
4using namespace std;
5
6int 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 ------------------------------
21a/b = [ 0.33333333333333333333333333333, 0.33333333333333333333333333334]
22------------------------------------------------------------------*/
The Multiple-Precision Data Type l_interval.

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"
8using namespace std;
9using namespace cxsc;
10
11l_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
17l_interval deriv(const l_interval& x) // Derivative function f'
18{
19 return 1/(2*sqrt(x)) + cos(x) - (x+1)*sin(x);
20}
21
22bool 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
27int 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
53Starting interval is [2,3])
54Diameter of actual enclosure: 1.000000
55Diameter of actual enclosure: 0.218137
56Diameter of actual enclosure: 0.013325
57Diameter of actual enclosure: 1.521930E-005
58Diameter of actual enclosure: 2.625899E-012
59Diameter of actual enclosure: 4.708711E-027
60Diameter of actual enclosure: 2.138212E-049
61Final enclosure of the zero:
62 [ 2.059045253415143788680636155343254522623083897,
63 2.059045253415143788680636155343254522623083898 ]
64*/
65
66
The Multiple-Precision Data Type l_real.
Definition l_real.hpp:78
cvector diam(const cimatrix_subv &mv) noexcept
Returns the diameter of the matrix.
Definition cimatrix.inl:738

Example 7 -

1// Runge-Kutta Method
2
3#include <iostream>
4#include "rvector.hpp" // Include dynamic arrays (real vectors)
5using namespace std;
6using namespace cxsc;
7
8rvector 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
18void 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
24int 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 ------------------------------
46Step: 1, x = 0.100000
47Y =
48 0.099747
49 0.995013
50 0.997400
51
52Step: 2, x = 0.200000
53Y =
54 0.197993
55 0.980203
56 0.989716
57
58Step: 3, x = 0.300000
59Y =
60 0.293320
61 0.956014
62 0.977286
63------------------------------------------------------------------*/
The Data Type rvector.
Definition rvector.hpp:58

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
7using namespace std;
8using namespace cxsc;
9
10int 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 ------------------------------
31Please enter the matrix dimension n: 3
32Please 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)
36Please 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)
40Trace of product matrix: ( 666.000000, 0.000000)
41------------------------------------------------------------------*/
The Data Type cdotprecision.
Definition cdot.hpp:61
The Data Type cmatrix.
Definition cmatrix.hpp:514
The Scalar Type complex.
Definition complex.hpp:50