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
• Example 2 - Input / Output
• Example 3 - Compute all zeros of a function
• Example 4 - Interval Newton method
• Example 5 -
• Example 6 -
• Example 7 -
• Example 8 -

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:

#include "interval.hpp"  // predefined interval arithmetic
#include <iostream>
using namespace cxsc;
using namespace std;

int main()
{
  interval a, b;            // Standard intervals     
  a = 1.0;                  // a   = [1.0,1.0]       
  b = 3.0;                  // b   = [3.0,3.0]        
  cout << "a/b = " << a/b << endl;
  return 0;
}

/* --------------------------- Output ------------------------------
a/b = [  0.333333,  0.333334]
------------------------------------------------------------------*/

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

#include <iostream>
#include "interval.hpp"
using namespace cxsc;
using namespace std;

int main()
{
  real a, b;

  cout << "Please enter real a: ";
  cout << RndDown;             // set rounding mode 
  cin  >> a;                   // read a rounded downwards                  
  cout << SetPrecision(7,4);   // set field width and number of digits 
                               // for output
  cout << a << endl;           // print a rounded downwards to 4 digits  
  cout << RndUp;              
  cout << a << endl;           // print a rounded upwards to 4 digits 
    
  "0.3" >> b;                  // convert the string "0.3" to a floating 
                               // point value b using rounding down    

  cout << SetPrecision(18,15); // from now on print 15 digits

  cout << b << endl;           // print b rounded upwards to 15 digits 
  cout << RndDown;
  cout << b << endl;           // print b rounded downwards to 15 digits 
  interval x;
  "[1.5, 2]" >> x;             // string to interval conversion 
  cout << x << endl;           // print interval x using default setting 

  cout << SaveOpt;             // push I/O parameters to internal stack 
  cout << SetPrecision(10,7);  // modifies output field width and 
                               // number of digits to be print 
  cout << x << endl;           // print x in the modified format 
  cout << RestoreOpt;          // pop parameters from internal stack   
  cout << x << endl;           // again, print x using the former format 
  return 0;
}

/* --------------------------- Output ------------------------------
Please enter real a: 0.3
 0.2999
 0.3000
 0.300000000000001
 0.300000000000000
[ 1.500000000000000, 2.000000000000000]
[ 1.5000000, 2.0000000]
[ 1.500000000000000, 2.000000000000000]
------------------------------------------------------------------*/

Example 3 - Compute all zeros of a function

// Compute all zeros of the function 
//
//              (x-1)*(exp(-3*x) - power(sin(x), 3))
//
// This example needs the CToolbox sources additionally!

#include "nlfzero.hpp"     // Nonlinear equations module
#include "stacksz.hpp"     // To increase stack size for some
                           // special C++ compiler
using namespace cxsc;
using namespace std;


DerivType f ( const DerivType& x )                    // Sample function
{ 
   return (x-1)*( exp(-3*x) - power(sin(x),3) ); 
}

// The class DerivType allows the computation of f, f', and f'' using
// automatic differentiation; see "C++ Toolbox for Verified Computing"

int main()
{
  interval  SearchInterval;
  real      Tolerance;
  ivector   Zero; 
  intvector Unique;
  int       NumberOfZeros, i, Error;

  cout << SetPrecision(23,15) << Scientific;   // Output format

  cout << "Search interval     : ";
  cin  >> SearchInterval;
  cout << "Tolerance (relative): ";
  cin  >> Tolerance;
  cout << endl;

  // Call the function 'AllZeros()' from the C++ Toolbox
  AllZeros(f,SearchInterval,Tolerance,Zero,Unique,NumberOfZeros,Error);

  for ( i = 1; i <= NumberOfZeros; i++) {
    cout << Zero[i] << endl;
    if (Unique[i])
      cout << "encloses a locally unique zero!" << endl;
    else
      cout << "may contain a zero (not verified unique)!" << endl;
  }
  cout << endl << NumberOfZeros << " interval enclosure(s)" << endl;
  if (Error) cout << endl << AllZerosErrMsg(Error) << endl;
  return 0;
}

/* --------------------------- Output ------------------------------
Search interval      : [-1,15]
Tolerance (relative) : 1e-10

[ 5.885327439818601E-001, 5.885327439818619E-001]
encloses a locally unique zero!
[ 9.999999999999998E-001, 1.000000000000001E+000]
encloses a locally unique zero!
[ 3.096363932404308E+000, 3.096363932416931E+000]
encloses a locally unique zero!
[ 6.285049273371415E+000, 6.285049273396501E+000]
encloses a locally unique zero!
[ 9.424697254738511E+000, 9.424697254738533E+000]
encloses a locally unique zero!
[ 1.256637410119757E+001, 1.256637410231546E+001]
encloses a locally unique zero!

6 interval enclosure(s)
------------------------------------------------------------------*/

Example 4 - Interval Newton method

// Interval Newton method using ordinary interval arithmetic
// Verified computation of a zero of the function f() 

#include <iostream>
#include "interval.hpp"          // Include interval arithmetic package
#include "imath.hpp"             // Include interval standard functions
using namespace std;
using namespace cxsc;

interval f(const real x)
{                                // Function f
  interval y(x);                 // y is a thin interval initialized by x
  return sqrt(y) + (y+1)*cos(y); // Interval arithmetic is used
}

interval deriv(const interval& x)
{                                // Derivative function f'
   return 1/(2*sqrt(x)) + cos(x) - (x+1)*sin(x);
}

bool criter(const interval& x)    // Computing: f(a)*f(b) < 0   and
{                                 //            not 0 in f'([x])?
   return Sup( f(Inf(x))*f(Sup(x)) ) < 0.0  &&  !(0.0 <= deriv(x)); 
}                                 // '<=' means 'element of'
 
int main(void)
{
   interval x, xOld;
   cout << SetPrecision(20,15); // Number of mantissa digits in I/O
   x= interval(2,3);
   cout << "Starting interval is [2,3]" << endl;
   if (criter(x))
   {  // There is exactly one zero of f in the interval x
      do {
         xOld = x;
         cout << "Actual enclosure is " << x << endl;
         x = (mid(x)-f(mid(x))/deriv(x)) & x; // Iteration formula
      } while (x != xOld);   
      cout << "Final enclosure of the zero: " << x << endl;
   } 
   else
      cout << "Criterion not satisfied!" << endl;
   return 0;
}

/* Output:

Starting interval is [2,3]
Actual enclosure is [   2.000000000000000,   3.000000000000000]
Actual enclosure is [   2.000000000000000,   2.218137182953809]
Actual enclosure is [   2.051401462380920,   2.064726329907714]
Actual enclosure is [   2.059037791936965,   2.059053011233253]
Actual enclosure is [   2.059045253413831,   2.059045253416460]
Actual enclosure is [   2.059045253415142,   2.059045253415145]
Final enclosure of the zero: [ 2.059045253415142, 
                               2.059045253415145]

*/

Example 5 -

#include "l_interval.hpp"  // interval staggered arithmetic
#include <iostream>
using namespace cxsc;
using namespace std;

int main() 
{
  l_interval a, b;         // Multiple-precision intervals 
  a = 1.0;
  b = 3.0;
  stagprec = 2;            // global integer variable      
  cout << SetDotPrecision(16*stagprec, 16*stagprec-3) << RndNext;
  // I/O for variables of type l_interval is done using
  // the long accumulator (i.e. a dotprecision variable)   

  cout << "a/b = " << a/b << endl;  
  return(0); 
}

/* --------------------------- Output ------------------------------
a/b = [ 0.33333333333333333333333333333, 0.33333333333333333333333333334]
------------------------------------------------------------------*/

Example 6 -

// Interval Newton method using a multi-precision interval data type
// Verified computation of a zero of the function f() to high accuracy

#include <iostream>             
#include "l_interval.hpp"         // Include multi-precision intervals
#include "l_imath.hpp"            // Include multi-precision math functions
#include "l_rmath.hpp"
using namespace std;
using namespace cxsc;

l_interval f(const l_real& x)     // Function f
{ 
   l_interval y(x);               // y is a thin interval initialized by x  
   return sqrt(y) + (y+1)*cos(y); // Use multi-precision interval arithmetic
}

l_interval deriv(const l_interval& x) // Derivative function f'
{                                 
   return 1/(2*sqrt(x)) + cos(x) - (x+1)*sin(x);
}

bool criter(const l_interval& x)  // Computing: f(a)*f(b) < 0   and
{                                 //            not 0 in f'([x])?
   return Sup( f(Inf(x))*f(Sup(x)) ) < 0.0  &&  !(0.0 <= deriv(x)); 
}                                 // '<=' means 'element of'

int main(void)
{
   l_interval x, xOld;
   stagprec= 3; // Set precision of the staggered correction arithmetic 
   x= l_interval(2,3);
   cout << "Starting interval is [2,3]" << endl;
   cout << SetDotPrecision(16*stagprec, 16*stagprec-3) << RndNext;
   // I/O for variables of type l_interval is done using
   // the long accumulator (i.e. a dotprecision variable)   

   if (criter(x)) 
   {  // There is exactly one zero of f in the interval x
      do {
         xOld = x;
         cout << "Diameter of actual enclosure: " << real(diam(x)) << endl;
         x = (mid(x)-f(mid(x))/deriv(x)) & x; // Iteration formula
      } while (x != xOld);                    // &  means  intersection
      cout << "Final enclosure of the zero: " << x << endl;
   }
   else    
      cout << "Criterion not satisfied!" << endl;
   return 0;
}

/* Output:

Starting interval is [2,3])
Diameter of actual enclosure:   1.000000
Diameter of actual enclosure:   0.218137
Diameter of actual enclosure:   0.013325
Diameter of actual enclosure: 1.521930E-005
Diameter of actual enclosure: 2.625899E-012
Diameter of actual enclosure: 4.708711E-027
Diameter of actual enclosure: 2.138212E-049 
Final enclosure of the zero: 
   [ 2.059045253415143788680636155343254522623083897,  
     2.059045253415143788680636155343254522623083898 ]
*/

Example 7 -

// Runge-Kutta Method

#include <iostream>
#include "rvector.hpp"             // Include dynamic arrays (real vectors)
using namespace std;
using namespace cxsc;

rvector F(const real x, const rvector& Y) // Function definition
{  
   rvector Z(3);                   // Constructor call
                               
   Z[1] =  Y[2]  * Y[3];           // F is independent of x
   Z[2] = -Y[1]  * Y[3];
   Z[3] = -0.522 * Y[1] * Y[2];
   return Z;
}

void Init(real& x, real& h, rvector& Y)
{                                  // Initialization
   x    = 0.0;   h    = 0.1;
   Y[1] = 0.0;   Y[2] = 1.0;  Y[3] = 1.0;
}

int main(void)
{
   real x, h;                      // Declarations and dynamic
   rvector Y(3), K1, K2, K3, K4;   // memory allocation
                                   // Automatic resize of Ki below 
   Init(x, h, Y);                                          
   for (int i=1; i<=3; i++) {             // 3 Runge-Kutta steps
      K1 = h * F(x, Y);                   // with array result
      K2 = h * F(x + h / 2, Y + K1 / 2);
      K3 = h * F(x + h / 2, Y + K2 / 2);
      K4 = h * F(x + h, Y + K3);
      Y  = Y + (K1 + 2 * K2 + 2 * K3 + K4) / 6;
      x += h;
      cout << SetPrecision(10,6) << Dec;  // I/O modification
      cout << "Step: " << i << ",  "
           << "x =   " << x << endl;
      cout << "Y =   " << endl << Y << endl;
   }
   return 0;
}

/* --------------------------- Output ------------------------------
Step: 1,  x =     0.100000
Y =   
  0.099747
  0.995013
  0.997400

Step: 2,  x =     0.200000
Y =   
  0.197993
  0.980203
  0.989716

Step: 3,  x =     0.300000
Y =   
  0.293320
  0.956014
  0.977286
------------------------------------------------------------------*/

Example 8 -

// Trace of a (complex) matrix product  
// Let C denote the matrix product A*B. 
// Then the diagonal entries of C are added to get the trace of C. 

#include <iostream>
#include "cmatrix.hpp"       // Use the complex matrix package
using namespace std;
using namespace cxsc;
 
int main()
{
   int n;
   cout << "Please enter the matrix dimension n: ";  cin  >> n;
   cmatrix A(n,n), B(n,n);   // Dynamic allocation of A, B
   cdotprecision accu;       // Allows exact computation of dotproducts
   cout << "Please enter the matrix A:" << endl;  cin  >> A;
   cout << "Please enter the matrix B:" << endl;  cin  >> B;
   accu = 0.0;               // Clear accumulator
   for (int i=1; i<=n; i++) accumulate(accu, A[i], B[Col(i)]);
   // A[i] and B[Col(i)] are subarrays of type rvector.

   // The exact result stored in the complex dotprecision variable accu 
   // is rounded to the nearest complex floating point number:
   complex result = rnd(accu);  
   cout << SetPrecision(12,6) << RndNext << Dec;
   cout << "Trace of product matrix: " << result << endl;
   return 0;
}

/* --------------------------- Output ------------------------------
Please enter the matrix dimension n: 3
Please enter the matrix A:
(1,0) (2,0) (3,0)
(4,0) (5,0) (6,0)
(7,0) (8,0) (9,0)
Please enter the matrix B:
(10,0) (11,0) (12,0)
(13,0) (14,0) (15,0)
(16,0) (17,0) (18,0)
Trace of product matrix: (  666.000000,    0.000000) 
------------------------------------------------------------------*/

---