/* Numerical.h
A C++ class for basic numerical calculations.
------------
Example usage for a function f(x).
Numrical N;
N.Diff(f, 5, 1) -> first derivative of f(x) at x=5.
N.Ints(f, 1, 4) -> integral of f(x) at between 1 and 4.
N.Root(f, 1) -> search for the root of f(x) around x=1
N.Opti(f, 1) -> search for the optimum value of f(x) around x=1
------------
Contact : bingul(at)gantep.edu.tr
ChangeLog
07.03.2014 The class is first written.
07.03.2014 Function-base methods are added.
11.03.2014 Array-base and vector-base methods are added.
*/
#ifndef NUMERICAL_H
#define NUMERICAL_H 1
#include <cmath>
#include <string>
#include <vector>
#include <iostream>
using namespace std;
/* **************************************************** */
class Numerical{
private:
double tolerance, h;
int maxiter;
public:
Numerical();
~Numerical();
void Mess(string s);
// Function-base methods
double Intt(double (*fun)(double x), double a=0.0, double b=1.0, int n=1000);
double Ints(double (*fun)(double x), double a=0.0, double b=1.0, int n=1000);
double Diff(double (*fun)(double x), double x, int n=1);
double Root(double (*fun)(double x), double x=1.0);
double Opti(double (*fun)(double x), double x=1.0);
// Array-base or vector-base methods
double Diff(double x[], double y[], int size, int index, int n=1);
double Intt(double x[], double y[], int size);
double Intt(vector<double> x, vector<double> y);
};
/* **************************************************** */
// -------------------------------------------------------------------------------
// The constructor function
Numerical::Numerical(){
tolerance = 1.0e-6;
h = 1.0e-3;
maxiter = 50;
}
// -------------------------------------------------------------------------------
Numerical::~Numerical(){
// empty
}
// -------------------------------------------------------------------------------
// Prints a message
void Numerical::Mess(string s){
cout << endl << "Numerical::Mess INFO \a" << s << endl;
}
/* Function-base methods */
// -------------------------------------------------------------------------------
// Retuns root of the function fun. x is the initial estimate.
// Numerical method: Newton-Rapson
double Numerical::Root(double (*fun)(double x), double x){
double err, xr = x;
int iter=0;
while(1){
err = fun(xr) / Diff(fun,xr,1);
xr = xr - err;
if(fabs(err)<tolerance) break;
iter++;
if(iter>maxiter) {
Mess("Loop does not converge. You may change initial estimate.");
break;
}
}
return xr;
}
// -------------------------------------------------------------------------------
// Retunrs optimum value of the function fun. x is the initial estimate
// Numerical method: Newton-Rapson
double Numerical::Opti(double (*fun)(double x), double x) {
double err, xo = x;
int iter=0;
while(1){
err = Diff(fun,xo,1) / Diff(fun,xo,2);
xo = xo - err;
if(fabs(err)<tolerance) break;
iter++;
if(iter>maxiter) {
Mess("Loop does not converge. You may change initial estimate.");
break;
}
}
return xo;
}
// -------------------------------------------------------------------------------
// Returns first or second derivative of the function at point x.
// Numerical method: Central Difference Appriximation
double Numerical::Diff(double (*fun)(double x), double x, int n){
if(n<1 || n>2) {
Mess("Order can be 1 or 2.");
return 0.0;
}
if(n==1) return (fun(x+h)-fun(x-h))/(2.0*h);
else return (fun(x+h)+fun(x-h)-2*fun(x))/(h*h);
}
// -------------------------------------------------------------------------------
// Returns integral of the function between a and b for n+1 equally-spaced points.
// Numerical method: Trapezoidal method
double Numerical::Intt(double (*fun)(double x), double a, double b, int n){
if(n<1) {
Mess("n must be at least 1.");
return 0.0;
}
double dx = (b-a)/n;
double etf = (fun(a)+fun(b))/2;
for (int i=1; i<n; i++){
etf += fun(a+i*dx);
}
return etf*dx;
}
// -------------------------------------------------------------------------------
// Returns integral of the function between a and b for n+1 equally-spaced points.
// Numerical method: Simpson's method
double Numerical::Ints(double (*fun)(double x), double a, double b, int n){
if(n<1 || n%2==1) {
Mess("n must be at least 1 and even integer.");
return 0.0;
}
double dx = (b-a)/n;
double esf = fun(a)+fun(b);
for (int i=1; i<=n-1; i++){
esf += 2.0 * pow(2.0,double(i%2)) * fun(a+i*dx);
}
return esf*dx/3.0;
}
/* Array-base or vector-base methods */
// -------------------------------------------------------------------------------
// Returns first or second derivative of array y of given size at given index.
// Numerical method: Central Difference Appriximation
double Numerical::Diff(double x[], double y[], int size, int index, int n){
if(n<1 || n>2) {
Mess("Order can be 1 or 2.");
return 0.0;
}
if(index<1 || index>size-2) {
Mess("Index value is out of range.");
return 0.0;
}
double dx = (x[index+1]-x[index-1])/2.0;
if(n==1) return (y[index+1]-y[index-1])/(2.0*dx);
else return (y[index+1]+y[index-1]-2*y[index])/(dx*dx);
}
// -------------------------------------------------------------------------------
// Returns integral of the array y of given size.
// Numerical method: Trapezoidal method
double Numerical::Intt(double x[], double y[], int size){
double etf = 0.0;
for(int i=0; i<size-1; i++){
etf += 0.5 * (y[i+1]+y[i]) * (x[i+1]-x[i]);
}
return etf;
}
// -------------------------------------------------------------------------------
// Returns integral of the vector v.
// Numerical method: Trapezoidal method
double Numerical::Intt(vector<double> x, vector<double> y){
if(x.size()!=y.size()){
Mess("Size of vectors must be the same.");
return 0.0;
}
double etf = 0.0;
for(unsigned int i=0; i<y.size()-1; i++){
etf += 0.5 * (y[i+1]+y[i]) * (x[i+1]-x[i]);
}
return etf;
}
#endif