I am trying to implement this function: , but it's not working. [EDIT] After user2079303's suggestion, PeterT's pointing out my error, and Toby Speight's suggestion, a mini-code that can be compiled to verify the validity of the function looks like this:
#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];
double den[countN];
for (int m=0; m<4; ++m)
{
double tmp {0.5 + 0.5 * (m > 0)};
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)) * tmp;
std::cout << tmp << "t" << num[i] << "t" << den[i] << "n";
}
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;
}
There is a minor simplification compared to the original: I factored 2 in both numerator and denominator and I used only the q^0.25
outside of the fraction. Also, countN
is the r
from the original document, micro
is only the 1/2
for even N
or 0
for odd N
, and i
is 0 for array's index but i+1
for calculations, but these shouldn't matter overall.
I tried this with 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;
And the results, in 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
and in 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
To my surprise, the first term is good, for bith N
, so I can't tell what in my code is not right. Could someone help me spot the error?
To be clear about it: I am a beginner in C++ and I am not looking for someone to do it for me, but to let me know of my erros in coding (translating math to C++ code).
Aucun commentaire:
Enregistrer un commentaire