C++ Theta реализация функции

Question

Я пытаюсь реализовать эту функцию:

, но это не работает. Минимальный, проверяемый пример выглядит следующим образом:

#include <iostream> 
#include <cmath> 

int main() 
{ 
    int N {8};  // change this for testing <1..inf> 
    double q {0.1/N}; 
    int countN {static_cast<int>(floor(N/2))}; 
    static const double PI {3.1415926535897932384626433832795}; 

    // Omega[i] = Theta1(u,m)/Theta4(u,m) 
    double Omega[countN]; 
    for (int i=0; i<countN; ++i) 
    { 
    double micro {!(N % 2) * 0.5}; // 0 for odd N, 1/2 for even N 
    double num[countN] {sin(PI/N * (i + 1 - micro))}; 
    double den[countN] {0.5}; 

    for (int m=1; m<4; ++m) 
    { 
     num[i] += pow(-1, m) * pow(q, m*(m + 1)) * sin((2 * m + 1) * PI/N * (i + 1 - micro)); 
     den[i] += pow(-1, m) * pow(q, m*m) * cos(2 * m * PI/N * (i + 1 - micro)); 

    } 
    Omega[i] = fabs(pow(q, 0.25) * num[i]/den[i]); 
    std::cout << " " << Omega[i] << "\n"; 
    } 

    // testing the values, they should be increasing in value 
    for (const auto &elem: Omega) 
    std::cout << elem << " "; 
    std::cout << "\n"; 

    return 0; 
} 

Существует незначительное упрощение по сравнению с оригиналом: Я факторизуюсь 2 как в числителе и знаменателе, и я использовал только q^0.25 вне фракции. Кроме того, countN является r из исходного документа, micro лишь 1/2 для даже N или 0 для нечетных N и i является 0 для индекса массива, но i+1 для вычислений, но это не имеет значения в целом.

Я попытался это с wxMaxima:

Theta[1](x,y):=2*y^0.25*sum((-1)^k*y^(k*(k+1))*sin((2*k+1)*x),k,0,n); 
Theta[4](x,y):=1+2*sum((-1)^k*y^(k^2)*cos(2*k*x),k,1,n); 
n:4$ 
N:8$ 
a:0.05$ 
b(i):=%pi/N*(i-(1-mod(N,2))/2)$ 
for N:8 thru 9 do for i:1 thru N/2 do print(["N=",N,"i=",i],Theta[1](b(i),a)/Theta[4](b(i),a)),numer; 

И результаты, в C++:

(q=0.05; N=8) 
Omega[0]=0.2018370065366672 
Omega[1]=0.06058232646142273 
Omega[2]=0.01205653570636574 
Omega[3]=0.02127667733703158 
(q=0.05; N=9) 
Omega[0]=0.348078726440638 
Omega[1]=0.1178366281313341 
Omega[2]=2.559808325080287e-07 
Omega[3]=0.02178788541277828 

и в wxMaxima:

["N=",8,"i=",1]" "0.2018370065366672" " 
["N=",8,"i=",2]" "0.5439269564954693" " 
["N=",8,"i=",3]" "0.7569342043740249" " 
["N=",8,"i=",4]" "0.850913653939989" " 
["N=",9,"i=",1]" "0.348078726440638" " 
["N=",9,"i=",2]" "0.6165773889432575" " 
["N=",9,"i=",3]" "0.7800391631077094" " 
["N=",9,"i=",4]" "0.8532352152763631 

К моему удивлению, первый член хорошо, для bith N, поэтому я не могу сказать, что в моем коде не так. Может ли кто-нибудь помочь мне обнаружить ошибку?

Чтобы это было ясно: я новичок в C++, и я не ищу, чтобы кто-то сделал это для меня, но чтобы сообщить мне о моих erros в кодировании (перевод математики на код C++).

Answer 1

Вы были

double den[countN] {0.5}; 

это инициализирует первый элемент den к 0.5и все остальные элементы 0.0 (инициализация по умолчанию). Другими словами, вышесказанное эквивалентно

double den[countN] {0.5, 0.0, 0.0, 0.0}; 

с таким количеством нулей, сколько необходимо для заполнения массива. Вероятно, вы хотели инициализировать все элементы до 0.5. В вашем случае, самый простой способ сделать это, когда вы первый использовать этот элемент - или, поскольку вы получаете доступ только один элемент den[i] в течение всего срока службы den, сделать ее простой double, а не массив:

for (int i=0; i<countN; ++i) { 
    double micro {N % 2 ? 0.0 : 0.5}; // 0 for odd N, 1/2 for even N 
    double num{sin(PI/N * (i + 1 - micro))}; 
    double den{0.5}; 

    for (int m=1; m<4; ++m) { 
     num += pow(-1, m) * pow(q, m*(m + 1)) * sin((2 * m + 1) * PI/N * (i + 1 - micro)); 
     den += pow(-1, m) * pow(q, m*m) * cos(2 * m * PI/N * (i + 1 - micro)); 
    } 
    Omega[i] = fabs(pow(q, 0.25) * num/den); 
}