Good time of day.
There is a function Phi(alpha) = sinc(pi*a/l*sin(alpha)) = sin(pi*a/l*sin(alpha))/(pi*a/l*sin(alpha)), domain is [-pi/2;pi/2].
I need to solve equation Phi(alpha) = c using the secant method. Default c is 0.7. As Phi is even function, we mostly work with non-negative space and search positive root.
Method requires the sin() realization and the choise of two start points. Also we need arcsin() and sqrt() and some other auxiliary functions. Phi(), sin(), arcsin() and sqrt() work properly and precisely (were tested separately). Secant method iterations are limited by Garwick's method – checking stop condition | alpha_n+1 - alpha_n | > | alpha_n - alpha_n-1 |. It could be modified, considering equality (stopping immediately if ">" and stopping after 10 consecutive "=="). But this modification doesn't give big advabtages (and doesn't affect the problem described below).
The problem is, solving function (Solve_Sinc) give good result (errors of order 10^(-7) — 10^(-9)) for one test examples, and very bad for others. Start points choise doesn't seem to affect it. The program output be like "error : solving result", bad examples are marked with comments in code.
So, the code is:
# include <iostream>
# include <conio.h>
using namespace std;
const float Pi = 3.1415927;
float abs(float);
float fact(int);
float pot_sgn(int);
float potent(float,int);
float rad_to_deg(float);
float sqrt(float);
float arcsin(float, int iter = 6);
float sin(float, int iter = 4);
float Phi(float, float);
float Solve_Sinc(float, float, float c = 0.7);
int main(void){
float a_size = 0.1;
cout << "__________" << endl;
cout << Solve_Sinc(a_size, 0.001, 0.35) << endl; //Incorrect: must be 0.00713124
cout << Solve_Sinc(a_size, 0.01, 0.35) << endl; //Incorrect: must be 0.0713724
cout << Solve_Sinc(a_size, 0.015, 0.35) << endl; //Incorrect: must be 0.107173
cout << Solve_Sinc(a_size, 0.04, 0.35) << endl; //Incorrect: must be 0.289264
cout << Solve_Sinc(a_size, 0.05, 0.35) << endl; //Incorrect: must be 0.364582
cout << Solve_Sinc(a_size, 0.06, 0.35) << endl;
cout << Solve_Sinc(a_size, 0.08, 0.35) << endl;
cout << Solve_Sinc(a_size, 0.1, 0.35) << endl; //Incorrect: must be ???
cout << Solve_Sinc(a_size, 0.11, 0.35) << endl; //Incorrect: must be ???
cout << Solve_Sinc(a_size, 0.3, 0.35) << endl;
cout << Solve_Sinc(a_size, 0.5, 0.35) << endl;
cout << "__________" << endl;
cout << Solve_Sinc(a_size, 0.001) << endl;
cout << Solve_Sinc(a_size, 0.01) << endl;
cout << Solve_Sinc(a_size, 0.015) << endl;
cout << Solve_Sinc(a_size, 0.04) << endl; //Incorrect: must be 0.180529
cout << Solve_Sinc(a_size, 0.05) << endl; //Incorrect: must be 0.226366
cout << Solve_Sinc(a_size, 0.06) << endl;
cout << Solve_Sinc(a_size, 0.08) << endl;
cout << Solve_Sinc(a_size, 0.1) << endl;
cout << Solve_Sinc(a_size, 0.11) << endl;
cout << Solve_Sinc(a_size, 0.3) << endl;
cout << Solve_Sinc(a_size, 0.5) << endl;
cout << "__________" << endl;
cout << Solve_Sinc(a_size, 0.001, 0.99) << endl; //Incorrect: must be 0.000780871
cout << Solve_Sinc(a_size, 0.01, 0.99) << endl; //Incorrect: must be 0.00780879
cout << Solve_Sinc(a_size, 0.015, 0.99) << endl; //Incorrect: must be 0.0117133
cout << Solve_Sinc(a_size, 0.04, 0.99) << endl; //Incorrect: must be 0.0312399
cout << Solve_Sinc(a_size, 0.05, 0.99) << endl; //Incorrect: must be 0.0390535
cout << Solve_Sinc(a_size, 0.06, 0.99) << endl; //Incorrect: must be 0.0468694
cout << Solve_Sinc(a_size, 0.08, 0.99) << endl; //Incorrect: must be 0.0625104
cout << Solve_Sinc(a_size, 0.1, 0.99) << endl; //Incorrect: must be 0.0781667
cout << Solve_Sinc(a_size, 0.11, 0.99) << endl; //Incorrect: must be 0.0860018
cout << Solve_Sinc(a_size, 0.3, 0.99) << endl; //Incorrect: must be 0.236459
cout << Solve_Sinc(a_size, 0.5, 0.99) << endl;
cout << "__________" << endl;
getch();
return 0;
}
float Solve_Sinc(float a, float l, float c){
if(a <= 0.0 || l <= 0.0 || c < 0.0){
cout << "Piston size and wavelength must be positive;\ncriterion must be non-negative!" << endl;
return 0.0;
}
float alph_old, alph_0, alph_1, k;
k = a/l*Pi;
alph_old = 0.0;
alph_0 = 0.0;
if(l/a > 1.0){
if(Phi(k,Pi/2.0) - c > 0.0){return Pi/2.0;}
alph_1 = Pi/2.0 ;
} else if(l/a > 0.5){
alph_1 = arcsin(l/a);
} else {
alph_1 = arcsin(2.0*l/a);
}
bool rep = true;
int eq = 0;
while(rep){ //Modified Garwick's method
alph_old = alph_0;
alph_0 = alph_1;
alph_1 = alph_0 - (Phi(k,alph_0) - c) * (alph_0 - alph_old) / (Phi(k,alph_0) - Phi(k,alph_old));
if(abs(alph_1 - alph_0) > abs(alph_0 - alph_old)){rep = false;}
else if(abs(alph_0 - alph_old) - abs(alph_1 - alph_0) <= 1e-6){eq += 1;}
else{eq = 0;}
if(eq == 10){rep = false;}
}
/*
do{
alph_old = alph_0;
alph_0 = alph_1;
alph_1 = alph_0 - (Phi(k,alph_0) - c) * (alph_0 - alph_old) / (Phi(k,alph_0) - Phi(k,alph_old));
}
while(abs(alph_1 - alph_0) < abs(alph_0 - alph_old)); //Common Garwick's method
*/
cout << Phi(k, alph_0) - c << " : ";
return alph_0;
}
float Phi(float k, float x){
if(abs(x) <= 1e-6){return 1.0;}
return sin(k*sin(x))/(k*sin(x));
}
float sin(float x, int iter){
if(x < 0.0){return -sin(-x);}
while (x >= 2.0*Pi){
x -= 2.0*Pi;
}
if(x > Pi){return -sin(2.0*Pi - x);}
if(x > Pi/2.0){return sin(Pi-x);}
float z = 0.0;
for(int i=iter;i>=0;i--){
z += pot_sgn(i)*(potent(x,2*i+1)/fact(2*i+1));
}
return z;
}
float arcsin(float x, int iter){
if(x < 0.0){return -arcsin(-x);}
if(x - (float)1.0 >= 1e-6){
cout << x << ":" << endl;
cout << "Absolute value of arcsin() argument can't be more than 1!" << endl;
return 0.0;
}
if(x >= 0.7071068){return Pi/2.0 - arcsin(sqrt(1.0 - x*x));}
if(x > 0.5 && x < 0.7071068){iter *= 2;}
float z = 0.0;
for(int i=iter;i>=0;i--){
z += (potent(x,2*i+1)*fact(2*i))/(fact(i)*fact(i)*potent(4.0,i)*(2*i+1));
}
return z;
}
float sqrt(float x){
float prec = 1e-6;
if(abs(x) < prec*10.0){return 0.0;}
if(x < 0.0){
cout << "Argument for sqrt() must be non-negative!" << endl;
return 0.0;
}
float y = x;
float delta;
float cnt = x;
if(x >= 1.0){
do{
cnt *= 0.1;
prec *= 10.0;
}
while(cnt >= 1.0);
}
do{
y = (y + x/y) / 2.0;
delta = x - y*y;}
while(delta < -prec || delta > prec);
return y;
}
float rad_to_deg(float x){
return (x/Pi)*180.0;
}
float potent(float x,int a){
float y = 1.0;
for(int i=1;i<=a;i++){
y *= x;
}
return y;
}
float pot_sgn(int n){
return (n%2 == 0) ? 1.0 : -1.0;
}
float fact(int n){
float f = 1.0;
while(n>1){
f *= n;
n -= 1;
}
return f;
}
float abs(float x){
return (x < 0.0) ? -x : x;
}
Why does it going on and how do I fix it (I mean making it give errors not more than 10^(-6) for all test examples in main())?
Thank you in advance!
