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:

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.

Example 2 - Input / Output

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 ------------------------------------------------------------------*/

Example 3 - Compute all zeros of a function

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 ------------------------------------------------------------------*/

Example 4 - Interval Newton method

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 */

Example 5 -

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 ------------------------------------------------------------------*/

Example 6 -

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 

Example 7 -

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 ------------------------------------------------------------------*/

Example 8 -

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 ------------------------------------------------------------------*/