Lesson N.8 of July 23, 2014

wrap a problem into a class

File example11.cc

:open:
 1/*
 2
 3  wrap a problem into a class
 4
 5*/
 6
 7#include <iostream>
 8#include <iomanip>
 9#include <vector>
10#include <cmath>
11
12class Problem {
13  double a ;
14public:
15  Problem( ) : a(0) { std::cout << "call constuctor 1\n" ; }
16  Problem( double _a ) : a(_a) { std::cout << "call constuctor 2\n" ; }
17  
18  double f( double x ) const { return x*x-a ; }
19  double df( double x ) const { return 2*x ; }
20
21} ;
22
23// Newton--Raphson
24bool
25Newton( Problem & pb,
26        double  & x, // x is the guess and at exit the approximate solution
27        int       max_iter,
28        double    tol ) {
29  for ( int i = 0 ; i < max_iter ; ++i ) {
30    double h = pb.f(x)/pb.df(x) ;
31    x -= h ;
32    if ( std::abs(h) < tol ) return true ;
33  }
34  return false ;
35}
36
37int
38main() {
39  Problem pb1(2) ;
40  double x = 2 ;
41  bool ok = Newton( pb1, x, 10, 1E-10 ) ;
42  if ( ok ) {
43    std::cout << "Newton converged, x = " << x << '\n' ;
44  } else {
45    std::cout << "Newton failed, x = " << x << '\n' ;
46  }
47  return 0 ;
48}

use template for the solver

File example12.cc

:open:
 1/*
 2
 3 Use template for the solver
 4
 5*/
 6
 7#include <iostream>
 8#include <iomanip>
 9#include <vector>
10#include <cmath>
11
12class Problem1 {
13  double a ;
14public:
15  Problem1( ) : a(0) { std::cout << "call constuctor 1\n" ; }
16  Problem1( double _a ) : a(_a) { std::cout << "call constuctor 2\n" ; }
17  
18  double f( double x ) const { return x*x-a ; }
19  double df( double x ) const { return 2*x ; }
20
21} ;
22
23class Problem2 {
24public:
25  Problem2( ) { std::cout << "call constuctor 1 (Problem2) \n" ; }
26  double f( double x ) const { return exp(x)-3 ; }
27  double df( double x ) const { return exp(x) ; }
28
29} ;
30
31// Newton--Raphson
32template <typename PROBLEM>
33inline
34bool
35Newton( PROBLEM & pb,
36        double  & x, // x is the guess and at exit the approximate solution
37        int       max_iter,
38        double    tol ) {
39  for ( int i = 0 ; i < max_iter ; ++i ) {
40    double h = pb.f(x)/pb.df(x) ;
41    x -= h ;
42    if ( std::abs(h) < tol ) return true ;
43  }
44  return false ;
45}
46
47int
48main() {
49  Problem1 pb1(2) ;
50  Problem2 pb2 ;
51  double x = 2 ;
52  bool ok = Newton( pb1, x, 10, 1E-10 ) ;
53  if ( ok ) {
54    std::cout << "Newton converged, x = " << x << '\n' ;
55  } else {
56    std::cout << "Newton failed, x = " << x << '\n' ;
57  }
58  x = 2 ;
59  ok = Newton( pb2, x, 10, 1E-10 ) ;
60  if ( ok ) {
61    std::cout << "Newton converged, x = " << x << '\n' ;
62  } else {
63    std::cout << "Newton failed, x = " << x << '\n' ;
64  }
65  return 0 ;
66}

Use inheritance

File example13.cc

:open:
 1/*
 2
 3  Inheritance and usage of it
 4
 5*/
 6
 7#include <iostream>
 8#include <iomanip>
 9#include <vector>
10#include <cmath>
11
12class Problem {
13public:
14  Problem( ) { std::cout << "call constuctor Problem class\n" ; }
15  virtual double f( double x ) const = 0 ;
16  virtual double df( double x ) const = 0 ;
17} ;
18
19class Problem1 : public Problem {
20  double a ;
21public:
22  Problem1( ) : a(0) { std::cout << "call constuctor Problem1\n" ; }
23  Problem1( double _a ) : a(_a) { std::cout << "call constuctor Problem1\n" ; }
24  
25  virtual double f( double x ) const { return x*x-a ; }
26  virtual double df( double x ) const { return 2*x ; }
27
28} ;
29
30class Problem2 : public Problem {
31public:
32  Problem2( ) { std::cout << "call constuctor 1 Problem2 \n" ; }
33
34  virtual double f( double x ) const { return exp(x)-3 ; }
35  virtual double df( double x ) const { return exp(x) ; }
36
37} ;
38
39// Newton--Raphson
40bool
41Newton( Problem * pb,
42        double  & x, // x is the guess and at exit the approximate solution
43        int       max_iter,
44        double    tol ) {
45  for ( int i = 0 ; i < max_iter ; ++i ) {
46    double h = pb->f(x)/pb->df(x) ;
47    x -= h ;
48    if ( std::abs(h) < tol ) return true ;
49  }
50  return false ;
51}
52
53int
54main() {
55  Problem1 pb1(2) ;
56  Problem2 pb2 ;
57  double x = 2 ;
58  bool ok = Newton( &pb1, x, 10, 1E-10 ) ;
59  if ( ok ) {
60    std::cout << "Newton converged, x = " << x << '\n' ;
61  } else {
62    std::cout << "Newton failed, x = " << x << '\n' ;
63  }
64  x = 2 ;
65  ok = Newton( &pb2, x, 10, 1E-10 ) ;
66  if ( ok ) {
67    std::cout << "Newton converged, x = " << x << '\n' ;
68  } else {
69    std::cout << "Newton failed, x = " << x << '\n' ;
70  }
71  return 0 ;
72}

A BVP solver

File bvp.hh

:open:
 1#ifndef BVP_HH
 2#define BVP_HH
 3
 4#include <iostream>
 5#include <iomanip>
 6#include <vector>
 7#include <cmath>
 8
 9/*
10
11 
12
13 */
14
15class BVP {
16public:
17  typedef double value_type ;
18
19private:
20  std::vector<value_type> x, y ;
21
22public:
23
24  BVP( ) { std::cout << "call constuctor BVP class\n" ; }
25
26  virtual ~BVP( ) { std::cout << "call descructor BVP class\n" ; }
27
28  virtual value_type p( value_type x ) const = 0 ;
29  virtual value_type q( value_type x ) const = 0 ;
30  virtual value_type r( value_type x ) const = 0 ;
31
32  void
33  solve( value_type a, value_type Ya,
34         value_type b, value_type Yb,
35         int        N ) ;
36
37  std::vector<value_type> const & get_x() const { return x ; }
38  std::vector<value_type> const & get_y() const { return y ; }
39
40} ;
41
42#endif

File bvp.cc

:open:
 1#include "bvp.hh"
 2
 3/*
 4
 5
 6 */
 7
 8/*
 9
10  Solve BVP
11  
12     y''(x) + p(x) y'(x) + q(x) y(x) = r(x)
13     y(a) = Ya, y(b) = Yb
14
15  =================
16
17   y[k-1]-2*y[k]+y[k+1]            y[k+1]-y[k-1]
18   --------------------  + p(x[k]) ------------- + q(x[k]) * y[k] = r(x[k]) 
19          h*h                           2*h
20
21   y[0] = Ya
22   y[N] = Yb
23
24  =================
25
26   y[k-1] * ( 1 - h*p(x[k])/2 ) +
27   y[k]   * ( -2 + (h*h)*q(x[k]) ) +
28   y[k+1] * ( 1 + h*p(x[k])/2 ) = (h*h)*r(x[k])
29
30  =================
31
32   alpha[k] * y[k-1] + beta[k] * y[k] + gamma[k] * y[k+1] = omega[k]
33
34   alpha[k] = 1 - h*p(x[k])/2 ;
35   beta[k]  = -2 + (h*h)*q(x[k]) ;
36   gamma[k] = 1 + h*p(x[k])/2 ;
37   omega[k] = (h*h)*r(x[k]) ;
38
39  =================
40
41*/
42
43using namespace std ;
44
45void
46BVP::solve( value_type a, value_type Ya,
47            value_type b, value_type Yb,
48            int        N ) {
49  vector<value_type> alpha(N+1), beta(N+1), gamma(N+1), omega(N+1) ;
50  value_type h = (b-a)/N ;
51  // equation of the first BC
52  alpha[0] = 0 ;
53  beta[0]  = 1 ;
54  gamma[0] = 0 ;
55  omega[0] = Ya ;
56  // equation of the second BC
57  alpha[N] = 0 ;
58  beta[N]  = 1 ;
59  gamma[N] = 0 ;
60  omega[N] = Yb ;
61  
62  // fill the lienear system
63  value_type h2 = h*h ;
64  x.resize(N+1) ;
65  y.resize(N+1) ;
66  x[0] = a ;
67  x[N] = b ;
68
69  for ( int k = 1 ; k < N ; ++k ) {
70    value_type xk = a+h*k  ;
71    value_type bf = h*p(xk)/2 ;
72    x[k]     = xk ;
73    alpha[k] = 1 - bf ;
74    beta[k]  = -2 + h2*q(xk)  ;  
75    gamma[k] = 1 + bf ;
76    omega[k] = h2*r(xk) ;
77  }
78
79  // solve linear system and compute y
80  // eliminate lower diagonal
81  for ( int k = 1 ; k <= N ; ++k ) {
82    value_type bf = alpha[k]/beta[k-1] ;
83    beta[k]  -= bf*gamma[k-1] ;
84    omega[k] -= bf*omega[k-1] ;
85  }
86
87  // solve bidiagonal system
88  y[N] = omega[N] / beta[N] ;
89  for ( int k = N-1 ; k >= 0 ; --k ) y[k] = (omega[k] - gamma[k]*y[k+1]) / beta[k] ;
90
91}

File bvp_test.cc

:open:
 1#include "bvp.hh"
 2#include <iomanip>
 3
 4using namespace std ;
 5
 6/*
 7  y'' + y' - 3*y = 9*exp(3*x) -6
 8
 9  y(0) = 3
10  y(1) = exp(3)+2
11
12  // exact solution
13  y(x) = exp(3*x) + 2 
14
15 */
16
17class my_BVP : public BVP {
18public:
19  virtual value_type p( value_type x ) const { return 1 ; }
20  virtual value_type q( value_type x ) const { return -3 ; }
21  virtual value_type r( value_type x ) const { return 9*exp(3*x) - 6 ; }
22  
23  value_type exact( value_type x ) const { return exp(3*x) + 2 ; }
24
25  void
26  solve( value_type a, value_type b, int N ) 
27  { BVP::solve( a, exact(a), b, exact(b), N ) ; }
28 
29};
30
31int
32main() {
33  my_BVP bvp ;
34  bvp.solve( 0, 1, 100 ) ;
35  vector<BVP::value_type> const & x = bvp.get_x() ;
36  vector<BVP::value_type> const & y = bvp.get_y() ;
37  cout << setw(4) << "index"
38       << " " << setw(12) << "x" 
39       << " " << setw(12) << "y"
40       << " " << setw(12) << "esatta"
41       << '\n' ;
42  for ( int i = 0 ; i <= 100 ; ++i )
43    cout << setw(4) << i
44         << " " << setw(12) << x[i] 
45         << " " << setw(12) << y[i]
46         << " " << setw(12) << bvp.exact(x[i])
47         << '\n' ;
48  return 0 ;
49}

All the files in one zip file

2014_07_23.zip