Взаимная ортогональность волновых функций - броненосец - PullRequest
Новая
       46

Взаимная ортогональность волновых функций - броненосец

0 голосов
/ 05 ноября 2018

Не уверен, что это лучшее место, чтобы спросить это.

Итак, я пытаюсь изучить ортогональность решений волновой функции, вычисляя интеграл произведения двух решений разных порядков m и n. Теперь я перехожу к той части, где мне нужно сделать произведение из двух матриц Эрмита разных размеров, которые я не могу математически выполнить, одна из которых 3х20, а другая 4х20. Есть ли способ обойти это? 

   arma::mat Orthonormality::gaussHermiteG(int n, int m, arma::mat Z)
{
    Miscellaneous misc;
    Calcul *caln = new Calcul(n,Z);
    Calcul *calm = new Calcul(m,Z);
    double f1;

    arma::mat Hnm;
    arma::mat res;
    f1 = (1 / std::sqrt(std::exp(n * std::log(2)) * misc.factorial(n))) * (1 / std::sqrt(std::exp(m * std::log(2)) * misc.factorial(m))) * std::sqrt(1 / M_PI);
    Hnm = caln->calculPolynomeHermite() % calm->calculPolynomeHermite();
    res = f1 * Hnm;

    return res;

}

Вот моя функция для получения квадратуры. Это способ доказать ортогональность, или я делаю это неправильно? 

long double Orthonormality::quadrature(int n, int m)
{
    arma::mat gx;   

    arma::mat gauss_point = {{
            -2.453407083009012499038365306336166239661e-1,
            2.453407083009012499038365306336166239661e-1,
            -7.374737285453943587056051442521042290772e-1,
            7.374737285453943587056051442521042290772e-1,
            1.234076215395323007885818346959410229585,
            -1.234076215395323007885818346959410229584,
            -1.738537712116586206780865662136406442958,
            1.738537712116586206780865662136406442953,
            2.254974002089275523082333344734565128082,
            -2.254974002089275523082333344734565128065,
            -2.788806058428130480525033756403185410695,
            2.788806058428130480525033756403185410655,
            3.347854567383216326914924522996463698566,
            -3.347854567383216326914924522996463698495,
            -3.94476404011562521037562880052441180715,
            3.944764040115625210375628800524411807067,
            4.603682449550744273077675248978347585171,
            -4.603682449550744273077675248978347585109,
            5.387480890011232862016900410681120753981,
            -5.387480890011232862016900410681120754003,
        }
    };

    arma::mat gauss_point_weight = {{
            4.622436696006100896503286398612081142142e-1,
            4.622436696006100896503286398612081142142e-1,
            2.866755053628341297196597062280879168236e-1,
            2.866755053628341297196597062280879168236e-1,
            1.090172060200233200137550335354255770852e-1,
            1.090172060200233200137550335354255770846e-1,
            2.481052088746361088216495255894039439922e-2,
            2.481052088746361088216495255894039440028e-2,
            3.24377334223786183218324713235370544232e-3,
            3.243773342237861832183247132353705443042e-3,
            2.283386360163539672571459179634955394906e-4,
            2.283386360163539672571459179634955393512e-4,
            7.802556478532063694145991999647569104495e-6,
            7.802556478532063694145991999647569095955e-6,
            1.086069370769281693999524563447163430255e-7,
            1.086069370769281693999524563447163432688e-7,
            4.399340992273180553628851455467928211995e-10,
            4.399340992273180553628851455467928212879e-10,
            2.229393645534151292522500616029095785758e-13,
            2.22939364553415129252250061602909578525e-13,
        }
    };


    gx = Orthonormality::gaussHermiteG(n, m, gauss_point);

    arma::mat res;
    res = gx * gauss_point_weight.t();
    long double resDouble = res(0, 0);
    return resDouble;
}

Вот функция полинома Эрмита и ее выход для 3 и 4 режимов: 

mat Calcul::calculPolynomeHermite(int n_max, mat z)
{   
    mat H(n_max, z.n_elem);

    if (n_max == 0)
    {
        H = z.ones(size(z));
    }

    else
    {
        if (n_max == 1)
        {
            return z.for_each([](arma::mat::elem_type& val)
            {
                val = 2 * val;
            });
        }
        else {
            for(int i = 0; i < z.n_elem; ++i)
            {
                H(0, i) = 1;
            }
            rowvec h2 = rowvec(z.n_elem);
            h2 = 2 * z;

            H.row(1) = h2;

            for(int i = 2; i < n_max; i++)
            {
                rowvec hn = rowvec(z.n_elem);
                hn = h2 % H.row(i - 1) - (2 * i) * H.row(i - 2);
                H.row(i) = hn;
            }
        }       
    }
    return H;
}

вывод: 

H(3,z):
    1.0000    1.0000    1.0000    1.0000    1.0000
   -4.0000   -2.0000         0    2.0000    4.0000
   12.0000         0   -4.0000         0   12.0000

H(4,z):
    1.0000    1.0000    1.0000    1.0000    1.0000
   -4.0000   -2.0000         0    2.0000    4.0000
   12.0000         0   -4.0000         0   12.0000
  -24.0000   12.0000         0  -12.0000   24.0000
...