Скорость сходимости не достигает ожидаемой (по решениям уравнения Бюргерса)

Question

Мне нужна помощь, чтобы проверить мой код на C ++, разработанный для решения уравнения Бюргерса с использованием схемы взвешенного ENO третьего порядка. Ожидаемая скорость конвергенции, конечно, 3 (три). Но мой код не соответствует этому, результат достигает только 1 (одного) коэффициента конвергенции. Код использует библиотеку Eigen для векторных операций. Скорость конвергенции вычисляется с использованием этого уравнения: Log10 (Ошибка (0) / Ошибка (1)) / Log10 (N (1) / N (0)).

#include <weno_bur.h>

void Test(const uint& n)
{
    uint k = 3;
    vxd X = vxd::LinSpaced(n+1,(double)-1,(double)1), U, U_Exact;
    U = ICs(X,k);

    double t = (double)0, tmax = 0.2, CFL = 0.6;
    double dx = (double)2/(n+1), dt;

    U_Exact = U_exact(U,X,tmax);

    while(t<tmax)
    {
        dt = DT(U,dx,CFL);

        if(t+dt > tmax){
            dt = tmax - t;
        }

        U = TimeStepping(U,dx,dt,k);

        t += dt;
    }

    std::cout << "Numerical errors (n="<< n <<"): " << (U-U_Exact).array().abs().sum()/n << std::endl;
}

int main()
{

    Test(32);
    Test(64);
    Test(128);

    return 0;
}

Файл weno_bur.h:

/*
 * weno_bur.h
 *
 *  Created on: Nov 8, 2019
 *      Author: rodbround
 */

#ifndef WENO_BUR_H_
#define WENO_BUR_H_

#include <iostream>
#include <fstream>
#include <iomanip>
#include <Eigen/Dense>

typedef unsigned int uint;
typedef Eigen::VectorXd vxd;

using namespace Eigen;

double pi(){
    return 3.141592653589793238462643383279502884;
}

vxd ICs(const vxd& x, const uint& k)
{
    uint n = x.size();
    vxd u = vxd::Zero(n);

    for(uint i = 0; i<n; i++)
    {
        u(i) = sin(pi()*x(i));
    }

    return u;
}

vxd Extend(const vxd& u, const uint& k)
{
    uint n = u.size();
    vxd U(n+k+k);
    U.segment(k,n) = u;
    U.segment(0,k) = u.segment(n-k-1,k).rowwise().reverse();
    U.segment(n+k,k) = u.segment(1,k).rowwise().reverse();

    return U;
}

vxd WENO3(const vxd& V)
{
    double Eps = 1e-15, Asum;
    vxd d(2),b(2),p(2),A(2),w(2),Vc(2);

    d(0) = (double)2/(double)3; d(1) = (double)1/(double)3;
    b(0) = std::pow(V(2) - V(1),2);
    b(1) = std::pow(V(1) - V(0),2);

    A(0) = d(0)/std::pow(b(0) + Eps,2);
    A(1) = d(1)/std::pow(b(1) + Eps,2);
    Asum = A.sum();

    w = A/Asum;

    p(0) = (V(1) + V(2))/(double)2;
    p(1) = (-V(0) + (double)3*V(1))/(double)2;

    Vc(0) = p.dot(w);

    A(0) = d(1)/std::pow(b(0) + Eps,2);
    A(1) = d(0)/std::pow(b(1) + Eps,2);
    Asum = A.sum();

    w = A/Asum;

    p(0) = ((double)3*V(1) - V(2))/(double)2;
    p(1) = (V(0) + V(1))/(double)2;

    Vc(1) = p.dot(w);

    return Vc;
}

vxd LF(const vxd& UL, const vxd& UR, const double& S)
{
    return ((UL.array().pow(2) + UR.array().pow(2))/(double)2 - \
            S*(UR.array() - UL.array()))/(double)2;
}

vxd RHS(const vxd& u, const double& dx, const uint& k)
{
    vxd U = Extend(u,k);
    uint n = U.size(), s = (k-1)/2, r = (k+1)/2, nf = n-s-s;
    vxd UC, UR(nf), UL(nf), FLUX;

    for(uint i=0; i<nf; i++)
    {
        UC = WENO3(U.segment(i,k));
        UL(i) = UC(0); UR(i) = UC(1);
    }

    FLUX = LF(UL.segment(0,n-k),UR.segment(1,n-k),U.maxCoeff());

    return (-(FLUX.segment(r,n-k-k) - FLUX.segment(r-1,n-k-k)))/dx;
}

double DT(const vxd& U, const double& dx, const double& CFL)
{
    return CFL*dx/U.maxCoeff();
}

vxd TimeStepping(const vxd& U, const double& dx, const double& dt, const uint& k)
{
    vxd U0 = U, U1, U2;
    U1 = U0 + RHS(U,dx,k)*dt;
    U2 = ((double)3*U0 + U1 + RHS(U1,dx,k)*dt)/(double)4;
    U0 = (U0 + (double)2*(U2 + RHS(U2,dx,k)*dt))/(double)3;

    return U0;
}

vxd U_exact(const vxd& U, const vxd& X, const double& tmax)
{
    double uo, un, err;
    vxd U__(U.size());

    for(unsigned int i = 0; i < U.size(); i++){
        uo = U(i);
        err = 1e100;
        while(err >= 1e-15){
            un = sin(pi()*(X(i) - uo*tmax));
            err = std::fabs(un - uo);
            uo = un;
        }
        U__(i) = uo;
    }

    return U__;
}

#endif /* WENO_BUR_H_ */

Выходные данные:

Numerical errors (n=32): 0.00854111
Numerical errors (n=64): 0.00393967
Numerical errors (n=128): 0.00179549

Исходя из приведенного выше результата, скорость схождения для n = 64 составляет 1.1191.