DERIVATIVE

template <class T>
class Difference{
    T a, b, x;
    T fun(T x){
        return cos(x/2);
    }
    public:
    Difference(T a, T b, T x){
        try{
            if(x<a || x>b){
                std::string boo {"Podane x nie nalezy do siatki.\n"};
                throw boo;
            }
            this->a = a;
            this->b = b;
            this->x = x;
        }catch(std::string &ex){
            std::cout<<ex;
            exit(0);
        }
    }

    //For derivatives where x is in inside node.
    T forawardDifference(T h){
        if(h!=0 && x>a && x<b){
            return (fun(x+h) - fun(x))/h;
        }
        return -100;
    }
    T backwardDifference(T h){
        if(h!=0 && x>a && x<b){
            return (fun(x)-fun(x-h))/h;
        }
        return -100;
    }
    T centralDifference(T h){
        if(h!=0 && x>a && x<b){
            return (fun(x+h) - fun(x-h))/(2*h);
        }
        return -100;
    }
    //To calculate derivative where x is in first or end node
    T firstOrEndNodeDerivativeTwoPoint(T h){
        if(h!=0 && x==a){
            return (fun(a+h)-fun(a))/h;
        }
        else if(h!=0 && x==b){
            return (fun(b)-fun(b-h))/h;
        }
        return -100;
    }
    T firstOrEndNodeDerivativeThreePoint(T h){
        if(h!=0 && x==a){
            return (-3.0/2.0*fun(a)+2*fun(a+h)-0.5*fun(a+2*h))/h;
        }
        else if(h!=0 && x==b){
            return (0.5*fun(b-2*h)-2*fun(b-h)+3.0/2.0*fun(b))/h;
        }
        return -100;
    }

    void changeX(T x){
        try{
            if(x<a || x>b){
                std::string x{"Podane x wychodzi poza siatke.\nX nie zostalo zmienione.\n"};
                throw x;
            }
            this->x = x;
        }catch(std::string &ex){
            std::cout<<ex;
        }
    }
};


MAIN

#include <iostream>
#include <cmath>
#include <vector>
#include <fstream>
#include <cmath>
#include "derivative.h"

#define M_PI 3.14159265359

template <typename T>
T difference(T x){
    return -sin(x/2.0)*1.0/2.0;
}

int main(){
    std::cout<<difference(0)<<" "<<difference(M_PI)<<" "<<difference(M_PI/2.0)<<std::endl;
    const double stepD = 0.0000000001;
    const double difABD = difference(M_PI/2.0);
    Difference<double> diffD(0.0, M_PI, (M_PI)/2);
    std::vector<double> forwDiffD, backwDiffD, centralDiffD, firstDiffTwoPointD,
                  firstDiffThreePointD, lastDiffTwoPointD, lastDiffThreePointD;
    for(double h = stepD; h<1.4; h+=stepD*100000){
        forwDiffD.push_back(fabs(difABD-diffD.forawardDifference(h)));
        backwDiffD.push_back(fabs(difABD-diffD.backwardDifference(h)));
        centralDiffD.push_back(fabs(difABD-diffD.centralDifference(h)));
    }
    const double difAD = difference(0);
    diffD.changeX(0.0);
    for(double h = stepD; h<1.4; h+=stepD*100000){
        firstDiffTwoPointD.push_back((difAD-diffD.firstOrEndNodeDerivativeTwoPoint(h)));
        firstDiffThreePointD.push_back(fabs(difAD-diffD.firstOrEndNodeDerivativeThreePoint(h)));
    }
    const double difBD = difference(M_PI);
    diffD.changeX(M_PI);
    for(double h = stepD; h<1.4; h+=stepD*100000){
        lastDiffTwoPointD.push_back(fabs(difBD-diffD.firstOrEndNodeDerivativeTwoPoint(h)));
        lastDiffThreePointD.push_back(fabs(difBD-diffD.firstOrEndNodeDerivativeThreePoint(h)));
    }
    std::ofstream file1;
    file1.open("dataDouble.dat");
    file1<<"H"<<" "<<"ErrorForw"<<" "<<"ErrorBack"<<" "<<
           "ErrorCent"<<" "<<"ErrorFirstTw"<<" "<<"ErrorFirstTh"<<
           " "<<"ErrorEndTw"<<" "<<"ErrorEndTh"<<std::endl;
    int i = 0;
    for(double h = stepD; h<1.4; h+=stepD*100000){
        file1 << h << " "<<forwDiffD.at(i);
        file1 <<" "<<backwDiffD.at(i);
        file1 <<" "<<centralDiffD.at(i);
        file1 <<" "<<firstDiffTwoPointD.at(i);
        file1 <<" "<<firstDiffThreePointD.at(i);
        file1 <<" "<<lastDiffTwoPointD.at(i);
        file1 <<" "<<lastDiffThreePointD.at(i)<<std::endl;
        i++;
    }
    file1.close();
    std::cout<<difference(0)<<" "<<difference(M_PI)<<" "<<difference(M_PI/2.0)<<std::endl;
    const float step = 0.0000000001;
    const float DIFAB = difference(M_PI/2.0);
    Difference<float> diff(0.0, M_PI, (M_PI)/2);
    std::vector<float> forwDiff, backwDiff, centralDiff, firstDiffTwoPoint,
                  firstDiffThreePoint, lastDiffTwoPoint, lastDiffThreePoint;
    for(float h = step; h<1.4; h+=step*100000){
        forwDiff.push_back(fabs(DIFAB-diff.forawardDifference(h)));
        backwDiff.push_back(fabs(DIFAB-diff.backwardDifference(h)));
        centralDiff.push_back(fabs(DIFAB-diff.centralDifference(h)));
    }
    const float DIFA = difference(0);
    diff.changeX(0.0);
    for(float h = step; h<1.4; h+=step*100000){
        firstDiffTwoPoint.push_back((DIFA-diff.firstOrEndNodeDerivativeTwoPoint(h)));
        firstDiffThreePoint.push_back(fabs(DIFA-diff.firstOrEndNodeDerivativeThreePoint(h)));
    }
    const float DIFB = difference(M_PI);
    diff.changeX(M_PI);
    for(float h = step; h<1.4; h+=step*100000){
        lastDiffTwoPoint.push_back(fabs(DIFB-diff.firstOrEndNodeDerivativeTwoPoint(h)));
        lastDiffThreePoint.push_back(fabs(DIFB-diff.firstOrEndNodeDerivativeThreePoint(h)));
    }
    file1.open("dataFloat.dat");
    file1<<"H"<<" "<<"ErrorForw"<<" "<<"ErrorBack"<<" "<<
           "ErrorCent"<<" "<<"ErrorFirstTw"<<" "<<"ErrorFirstTh"<<
           " "<<"ErrorEndTw"<<" "<<"ErrorEndTh"<<std::endl;
    i = 0;
    for(float h = step; h<1.4; h+=step*100000){
        file1 << h << " "<<forwDiff.at(i);
        file1 <<" "<<backwDiff.at(i);
        file1 <<" "<<centralDiff.at(i);
        file1 <<" "<<firstDiffTwoPoint.at(i);
        file1 <<" "<<firstDiffThreePoint.at(i);
        file1 <<" "<<lastDiffTwoPoint.at(i);
        file1 <<" "<<lastDiffThreePoint.at(i)<<std::endl;
        i++;
    }
    file1.close();
    return 0;
}