diff_eq_order_1.cpp

#include <iostream>
#include <string>
#include <vector>
#include <math.h>
#include <tuple>
#include "matplotlibcpp.h"

namespace plt = matplotlibcpp;

//-------------------------------------------------------------------

void plotNumericalApproach(std::vector<double> ii1, std::vector<double> yy1, std::vector<double> ii2, std::vector<double> yy2,std::vector<double> ii3, std::vector<double> yy3)
{
    plt::figure_size(1200, 780);
    plt::title("Differential equation. Solution");
    plt::plot(ii1, yy1);
    plt::plot(ii2, yy2);
    plt::plot(ii3, yy3);
    plt::show();
}
//-------------------------------------------------------------------
//used for extimation initial value
double solution(double t)
{

    int C = 1;
    return (C * std::exp(t) - t - 1);
}

//-------------------------------------------------------------------
//the function we approximate
double diffFunction(double t, double x)
{

    return x + t;
}

//-------------------------------------------------------------------

std::tuple<std::vector<double>, std::vector<double>> methodEuler(double (*function)(double t, double x), double runSize)
{

    std::vector<double> diffEq;
    std::vector<double> time;
    double xi;
    float dt = 0.5;
    float t = 0;
    double x0 = solution(0);

    diffEq.push_back(x0);
    time.push_back(0);

    for (int i = 1; i < runSize; i++)
    {
        xi = diffEq[i - 1] + dt * function(t, diffEq[i - 1]);
        diffEq.push_back(xi);
        time.push_back(i);
        t += dt;
    }

    return std::make_tuple(time, diffEq);
}

//-------------------------------------------------------------------

std::tuple<std::vector<double>, std::vector<double>> methodMidPoint(double (*function)(double t, double x), double runSize)
{

    std::vector<double> diffEq;
    std::vector<double> time;

    double xi;
    float dt = 0.5;
    float t = 0;
    double x0 = solution(0);

    diffEq.push_back(x0);
    time.push_back(0);
    for (int i = 1; i < runSize; i++)
    {
        xi = diffEq[i - 1] + dt * function(t + dt * 0.5, diffEq[i - 1] + 0.5 * dt * function(t, diffEq[i - 1]));
        diffEq.push_back(xi);
        time.push_back(i);
        t += dt;
    }

    return std::make_tuple(time, diffEq);
}

//-----------------------------------------------------------

std::tuple<std::vector<double>, std::vector<double>> methodRuneKutt(double (*function)(double t, double x), double runSize)
{

    std::vector<double> diffEq;
    std::vector<double> time;
    double xi, k1, k2, k3, k4, k5, k6;
    float dt = 0.5;
    float t = 0;
    double x0 = solution(0);

    time.push_back(0);
    diffEq.push_back(x0);

    for (int i = 1; i < runSize; i++)
    {
        k1 = dt * function(t, diffEq[i - 1]);
        k2 = dt * function(t + 0.5 * dt, diffEq[i - 1] + 0.5 * k1);
        k3 = dt * function(t + 0.5 * dt, diffEq[i - 1] + 0.5 * k2);
        k4 = dt * function(t + dt, diffEq[i - 1] + k3);
        xi = diffEq[i-1] + 1.0/6.0*(k1 + 2*k2 + 2*k3 + k4 );

        diffEq.push_back(xi);
        time.push_back(i);
        t += dt;
    }

    return std::make_tuple (time, diffEq);
}
//-------------------------------------------------------------------

int main()

{

    std::tuple<std::vector<double>, std::vector<double>>  diffSolutionEuler = methodEuler(diffFunction, 6);
    std::vector<double> timeEuler = std::get<0>(diffSolutionEuler);
    std::vector<double> solutionEuler = std::get<1>(diffSolutionEuler);
 
    for (auto &i : std::get<1>(diffSolutionEuler))
    {

        std::cout << i << std::endl;
    }

    std::cout << " -------------------------------------------- " << std::endl;

    std::tuple<std::vector<double>, std::vector<double>> diffSolutionMidPoint = methodMidPoint(diffFunction, 6);
    std::vector<double> timeMP = std::get<0>(diffSolutionMidPoint);
    std::vector<double> solutionMidPoint = std::get<1>(diffSolutionMidPoint);
 
    for (auto &i : std::get<1>(diffSolutionMidPoint))
    {

        std::cout << i << std::endl;

    }

    std::cout << " -------------------------------------------- " << std::endl;

    std::tuple<std::vector<double>, std::vector<double>>  diffRuneKutt = methodRuneKutt(diffFunction, 6);
    

    std::vector<double> time = std::get<0>(diffRuneKutt);
    std::vector<double> solutionRuneKutt = std::get<1>(diffRuneKutt);
    plotNumericalApproach(timeEuler, solutionEuler, timeMP, solutionMidPoint, time, solutionRuneKutt);


    for (auto &i : std::get<1>(diffRuneKutt))
    {

        std::cout << i <<  std::endl;

    }

}