Метод RKF45? - C++ - Ответ 9055556

#include <stdio.h>
#include <iostream>
#include <math.h>
using namespace std;
 
double const tol = 0.001;
 
double f1(double x1, double x2, double t)// 
{
    return(-52*x1-100*x2+exp(-t));               
}
 
double f2(double x1, double t)//
{
    return(x1+sin(t));               
} 
double s(double h, double x, double y)//
{
   return pow((tol*h)/(2*(fabs(y-x))),1.0/4.0);
}
 
int main()
{
    setlocale(LC_ALL,"rus");
   
    double t, s1,s2, h,tol, eps, E1, E2, k1, k2, k3, k4, k5, k6, y1, y2, m1, m2, m3, m4, m5, m6, x1, x2, x11,x22,y11,y22;
    //cout << "  "<<endl;
    //cout<<"h = ";//
    //cin >> h;
    
    t = 0; // 
    x1 = 1; // 
 
    x2 = 0; // 
    h = 0.001;   //2 or 3
    eps = 0.0001;//5 or 4
    int i = 0;
    y1 = x1;
    y2 = x2; 
    //s1 = s(h,x1,y1);
    //s2 = s(h,x2,y2);
    //while(E1 > eps || E2 > eps)
    while(true)
    {
 
        k1 = h*f1(x1, x2, t);
        m1 = h*f2(x1, t);
        
        k2 = h*f1(x1 + k1/4.0, x2 + m1/4.0, t + h/4.0);
        m2 = h*f2(x1 + k1/4.0, t + h/4.0);
 
        k3 = h*f1(x1 + (3.0*k1)/32.0 + (9.0*k2)/32.0, x2 + (3.0*m1)/32.0 + (9.0*m2)/32.0, t + (3.0*h)/8.0);
        m3 = h*f2(x1 + (3.0*k1)/32.0 + (9.0*k2)/32.0, t + (3.0*h)/8.0);
 
        k4 = h*f1(x1 + (1932.0*k1)/2197.0 - (7200.0*k2)/2197.0 + (7296.0*k3)/2197.0,
                  x2 + (1932.0*m1)/2197.0 - (7200.0*m2)/2197.0 + (7296.0*m3)/2197.0, t + (12.0*h)/13.0);
        m4 = h*f2(x1 + (1932.0*k1)/2197.0 - (7200.0*k2)/2197.0 + (7296.0*k3)/2197.0, t + (12.0*h)/13.0);
        
        k5 = h*f1(x1 + (439.0*k1)/216.0 - 8.0*k2 + (3680.0*k3)/513.0 - (815*k4)/4104.0,
                  x2 + (439.0*m1)/216.0 - 8.0*m2 + (3680.0*m3)/513.0 - (815*m4)/4104.0, t + h);
        m5 = h*f2(x1 + (439.0*k1)/216.0 - 8.0*k2 + (3680.0*k3)/513.0 - (815*k4)/4104.0, t + h);
 
        k6 = h*f1(x1 - (8.0*k1)/27.0 + 2*k2 - (3544.0*k3)/2565.0 + (1895.0*k4)/4104 - (11.0*k5)/40.0,
                  x2 - (8.0*m1)/27.0 + 2*m2 - (3544.0*m3)/2565.0 + (1895.0*m4)/4104 - (11.0*m5)/40.0,
                  t + h/2.0);
        m6 = h*f2(x1 - (8.0*k1)/27.0 + 2*k2 - (3544.0*k3)/2565.0 + (1895.0*k4)/4104 - (11.0*k5)/40.0, t + h/2.0);
 
        x11 = x1;
        x22 = x2;
        
        //до 5 
        x1 = x1 + (25.0*k1)/216.0 + (1408.0*k3)/2565.0 + (2197.0*k4)/4101.0 - k5/7.0;
        x2 = x2 + (25.0*m1)/216.0 + (1408.0*m3)/2565.0 + (2197.0*m4)/4101.0 - m5/7.0; 
        
        //cout<<"t = "<< t << " x1 = "<< x1 << " x2 = " << x2 << endl;
        //до 6
        y11 = x11;
        y22 = x22;
 
        y1 = y1 + (16.0*k1)/135.0 + (6656.0*k3)/12825.0 + (28561*k4)/56430.0 - (9*k5)/50.0 + (2.0*k6)/55.0; //???
        y2 = y2 + (16.0*m1)/135.0 + (6656.0*m3)/12825.0 + (28561*m4)/56430.0 - (9*m5)/50.0 + (2.0*m6)/55.0; //???
        //if(fabs(x1-x11)< 0.0001||fabs(x2-x22)< 0.0001)
        //{
        //  break;
        //}
        if((eps<=fabs(y1 - x1) && fabs(y1 - x1) < eps*10.0) && (eps<=fabs(y1 - x1)&&fabs(y1 - x1) < eps*10.0))
        {
            cout<<"Ответ:"<<endl;
            cout<<"t = "<< t << " x1 = "<< x1<< " y1 = "<< y1<< " x2 = "<< x2<< " y2 = "<< y2<<endl;
            break;
        }
        else if(fabs(y1 - x1) == eps/10.0||fabs(y1 - x1) < eps/10.0)
        {
            //h++;
            h = h*s(h,x1,y1);
            t = t + h;
 
        }
        else
        {
            //h--;
            h = h*s(h,x1,y1);
            //t = t + h;
            x1 = x11;
            x2 = x22;
            y1 = y11;
            y2 = y22;
        }
        //cмотрим шаг, если погрешность не устраивает . то меняем шаг и тогда удаляем  последнюю итерацию(вроде)
        //если погрешность устраивает то шаг оставляем и продолжаем.
        //E1 = ();
        //cout<<"E1 = "<<E1<<"  E2 = "<<E2;
        
    }
    
    //cout<<"t = "<< t << " x1 = "<< x1 << " x2 = " << x2 << endl;
    
    system("pause");
    return 0;
}