Runge-Kutta adaptive method (Метод рунге-Кутты-Мерсона с переменным шагом. Проблема.)

[yuran1987]

Привет всем. Проблема с методом Рунге-Кутты-Мерсона. Не могу понять как правильно реализовать данный метод т.к. количество узлов все время меняется(из-за шага). И решение не сходиться с решением классического метода Рунге-Кутты 4.
Объект trj - это просто структура с одномерными массивами для пяти решений системы дифф.уравнений.
Где ошибка может кто подскажет?!

Код: Выделить всё

class Solver{
public:

std::valarray<double> step_rkm45(const double &t, const std::valarray<double> &u, double h);

protected:
    virtual int dimention()=0;
    //!< \brief Метод определяет размерность системы дифф. уравнений.
    virtual std::valarray<double> f(const double&, const std::valarray<double>&)=0;
    //!< \brief Система уравнений.
private:
    int dim;
    std::valarray<double> k1, k2, k3, k4, k5, err, result;
    //!< Коэффиценты Рунге-Кутты
    double h;
    //!< Шаг равномерной сетки по независимой переменной
    double eps;
    //!< Погрешность на шаге для метода Рунге-Кутта-Мерсона
};

void Solver::setParam(const double &step, const double &err_max)
{
    h = step;
    eps = err_max;
    dim = dimention();
    k1.resize(dim);
    k2.resize(dim);
    k3.resize(dim);
    k4.resize(dim);
    k5.resize(dim);
    err.resize(dim);
    result.resize(dim);
}

std::valarray<double> Solver::step_rkm45(const double &t, const std::valarray<double> &u, double h)
{
    double epsmax = 1E-18;


    #pragma omp parallel for
    for(int i = 0; i < dim; i++)
        k1[i] = f(t,          u)[i];

    #pragma omp parallel for
    for(int i = 0; i < dim; i++)
        k2[i] = f(t + h/3.0,  u + k1*h/3.0)[i];

    #pragma omp parallel for
    for(int i = 0; i < dim; i++)
        k3[i] = f(t + h/3.0,  u + k1*h/6.0  +  k2*h/6.0)[i];

    #pragma omp parallel for
    for(int i = 0; i < dim; i++)
        k4[i] = f(t + h/2.0,  u + k1*h/8.0  +  3.0*k3*h/8.0)[i];


    #pragma omp parallel for
    for(int i = 0; i < dim; i++)
        k5[i] = f(t + h,  u + k1*h/2.0  - 3.0*k2*h/2.0 +2.0*k4*h)[i];

    #pragma omp parallel for
    for(int i = 0; i < dim; i++)
        err[i] = 0.2*fabs((h*(k1[i] + 4*k4[i] + k5[i])/6 - h*(k1[i] -3*k3[i] +4*k4[i])/2));

    if(err.sum()/dim > epsmax){
        h/=2.0;
        step_rkm45(t,u,h);
    }

    #pragma omp parallel for
    for(int i = 0; i < dim; i++)
        result[i] = u[i] + h/6.0 * (k1[i] + 4*k4[i] + k5[i]);

    if(err.sum()/dim < epsmax/64.0)
        h*=2.0;

    return result;
}

Функция f переопределяется в производном классе TrajectPl, вот так вот:

Код: Выделить всё

std::valarray<double> TrajectPl::f(const double &t, const std::valarray<double> &u)
{
    static std::valarray<double> tmp(dimention());
    tmp[0] = u[1];
    tmp[1] = -Z1neg*Uspx(u[0],t)/Vmax + kt*t + k0;
    tmp[2] = 2.*u[3];
    tmp[3] = u[4] - dUxx(u[0],t)*u[2]/Vmax;
    tmp[4] = epsilon*d*Diffus(u[0])/(Vmax*Vmax) - 2*dUxx(u[0],t)*u[3]/Vmax;
    return tmp;
}


int main{
TrajectPl s;

double step = (Tmax - Tmin)/Node;
s.setParam(step,1E-7);
double t= Tmin;

int n=0;
        while(t<Tmax){
            y0 = step_rkm45(trj.T[n],y0, step);
            trj.T[n+1] = trj.T[n] + step;
            trj.X[n+1] = y0[0];
            trj.Y[n+1] = y0[1];
            trj.MA[n+1] = y0[2];
            trj.MB[n+1] = y0[3];
            trj.MC[n+1] = y0[4];
            t+=step;
            //fprintf(stderr,"n = %i \n",n);
            n++;
        }
}

[RasenHerz]

Оформите код как полагается (с помощью тега code и расставления отступов), т.к. выкладывая его в таком виде вы проявляете неуважение к пользователям форума.