ObviousNistrem4(main), Implicit Runge Kutta 1(auxiliary) on

这是我所做的，它通过使用一步 RK 找到 3 个第一个值，但随后我有一些问题如何在这里使用 Nistrema4th。

在第一部分中，我得到eps，然后使用一步 Runge Kutta 并将这 3 个值作为具有最佳h（这是最小步长）的第一、第二和第三个值放入数组中。然后我需要从第 4 个开始并使用 Nistrema 找到其他人。

double fx(double x, double y, double t)
{
    return (2*x)-y+(t*t)-2*(sin(t)+1)+cos(t);
}
double fy(double x, double y, double t)
{
    return x+(2*y)-sin(t)-(2*(t*t))+(2*t)-1;
}
double precise_x(double t)
{
    return sin(t)+1;
}
double precise_y(double t)
{
    return (t*t);
}
int main(int argc, char* argv[])
{
    double h=0.025,h_opt, x[1000],y[1000],x_pre[10000],y_pre[10000],eps,current_x,current_y,x_prev,y_prev;
    double a=0.,b=1.,t=a, x_prev1,y_prev1,current_x1,current_y1, h_temp[3],x1,y1;
    int k=20,i;
    cout << "Input value of the epsilon, please: ";
    cin >> eps;
    a = 0;  b = 1;  t = a; x[0] = 1; y[0] = 0;
    current_x = x[0];
    current_y = y[0];
    int flag = 0;
    current_x = x[0];
    current_y = y[0];
    while(!flag)
    {
        for(i=0;i<k;i++)
        {
            x_prev = current_x;
            y_prev = current_y;
            current_x = x_prev + (h*fx(x_prev,y_prev,t+h));
            current_y = y_prev + (h*fy(x_prev,y_prev,t+h));
            x[0]=x[1];
            x[1]=x[2];
            x[2]=x[3];
            x[3]=current_x;
            y[0]=y[1];
            y[1]=y[2];
            y[2]=y[3];
            y[3]=current_y;
            t+=h;
        }
        if((fabs(x_prev-current_x)+fabs(y_prev-current_y)<eps))
        {
            flag = 1;
        }
        else
        {
            h/=2;
            t=0;
            current_x=1.;
            current_y=0.;
        }
        h_opt=h;
    }
 //===============???????????????????????????????????????????===
    //use formulas of Nisterm (4th order ) to get the results
    double n = 1./h_opt;
    for(i=0;i<=n;i++)
    {
        x[i+4] = x[i+2] + (h_opt/3.0)*(8*fx(x[i+3],y[i+3],(i+3)*h)
                           -(5*fx(x[i+2],y[i+2],(i+2)*h))
                           +4*fx(x[i+1],y[i+1],(i+1)*h)
                           -fx(x[i],y[i],i*h));
        y[i+4] = y[i+2] + (h_opt/3.0)*(8*fy(x[i+3],y[i+3],(i+3)*h)
                           -(5*fy(x[i+2],y[i+2],(i+2)*h))
                           +4*fy(x[i+1],y[i+1],(i+1)*h)
                           -fy(x[i],y[i],i*h));
    }
    t = a;
    for(i=0;i<=n;i++)
    {
        cout<<"n j "<<j;
        x_pre[i] = precise_x(t);
        y_pre[i] = precise_y(t);
        t += h_opt;
    }
    cout << "xtyttprecise_x precise_y" << endl;
    for(int i = 0; i <=k; i++)
    {
        cout  << x[i] << "t" << y[i]<<"t"<< x_pre[i] <<"t"<<y_pre[i]<<endl;
    }    
    system("PAUSE");
    return 0;
}