Решение трансцендентного уравнения методом подвижных секущих

[Newbie_cpp]

Всем доброго времени суток.

Дана функция Phi(alpha) = sinc(pi*a/l*sin(alpha)) = sin(pi*a/l*sin(alpha))/(pi*a/l*sin(alpha)), область определения [-pi/2;pi/2].

Необходимо решить уравнение Phi(alpha) = c методом подвижных секущих (он же метод хорд). Параметр c по умолчанию равен 0,7. Т.к. функция чётная, работа в основном идёт с неотрицательной областью, и корень ищется положительный.

Метод предполагает вычисление синуса, а также выбор двух разгонных (т.е. начальных) точек. Также необходим арксинус, квадратный корень и некоторые другие вспомогательные функции. Реализации "фи", синуса, арксинуса и корня работают исправно и точно (проверялись по отдельности). Итерации метода хорд ограничены так называемым приёмом Гарвика – по сути, это проверка условия остановки | alpha_n+1 - alpha_n | > | alpha_n - alpha_n-1 |. Можно также его модифицировать, учитывая равенство (прекращая сразу при ">" и прекращая после 10 идущих подряд "="). Ярко выраженного эффекта такая модификация, правда, не даёт (и не влияет на нижеописанную проблему).

Проблема заключается в том, что в одних случаях программа решающая функция (Solve_Sinc) даёт хороший результат (погрешность порядка 10^(-7) — 10^(-9)), а в других – совсем неправильный. Выбор разгонных точек на это не влияет (ошибки наблюдаются при разных вариантах). Вывод программы осуществляется в формате "погрешность : результат решения", а в комментариях к коду помечены ошибочные примеры.

Итак, код:

# 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; //Неправильно: должно быть 0.00713124
cout « Solve_Sinc(a_size, 0.01, 0.35) « endl; //Неправильно: должно быть 0.0713724
cout « Solve_Sinc(a_size, 0.015, 0.35) « endl; //Неправильно: должно быть 0.107173
cout « Solve_Sinc(a_size, 0.04, 0.35) « endl; //Неправильно: должно быть 0.289264
cout « Solve_Sinc(a_size, 0.05, 0.35) « endl; //Неправильно: должно быть 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; //Неправильно: должно быть ???
cout « Solve_Sinc(a_size, 0.11, 0.35) « endl; //Неправильно: должно быть ???
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; //Неправильно: должно быть 0.180529
cout « Solve_Sinc(a_size, 0.05) « endl; //Неправильно: должно быть 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; //Неправильно: должно быть 0.000780871
cout « Solve_Sinc(a_size, 0.01, 0.99) « endl; //Неправильно: должно быть 0.00780879
cout « Solve_Sinc(a_size, 0.015, 0.99) « endl; //Неправильно: должно быть 0.0117133
cout « Solve_Sinc(a_size, 0.04, 0.99) « endl; //Неправильно: должно быть 0.0312399
cout « Solve_Sinc(a_size, 0.05, 0.99) « endl; //Неправильно: должно быть 0.0390535
cout « Solve_Sinc(a_size, 0.06, 0.99) « endl; //Неправильно: должно быть 0.0468694
cout « Solve_Sinc(a_size, 0.08, 0.99) « endl; //Неправильно: должно быть 0.0625104
cout « Solve_Sinc(a_size, 0.1, 0.99) « endl;  //Неправильно: должно быть 0.0781667
cout « Solve_Sinc(a_size, 0.11, 0.99) « endl; //Неправильно: должно быть 0.0860018
cout « Solve_Sinc(a_size, 0.3, 0.99) « endl; //Неправильно: должно быть 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){
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)); // Обычный приём Гарвика
*/

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)2*i)/(fact(i)i4.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;
}

Почему возникает эта проблема и как её устранить (т.е. чтобы для всех тестовых примеров в main() погрешность была не более 10^(-6))?
Заранее спасибо!

[Алексей1153]

для начала - исправить явные косяки:

0)
добавить форматирование кода

1)

Сообщение от Newbie_cpp

cout « Solve_Sinc

2)
float -> double

3)

Сообщение от Newbie_cpp

const float Pi = 3.1415927;
float abs(float);
float sqrt(float);
float arcsin(float, int iter = 6);
float sin(float, int iter = 4);

вместо костылей использовать стандартные функции из std:: (#include <cmath>)

[Newbie_cpp]

Алексей1153, Насчёт форматирования: если Вы имеете в виду отступы, то это сам сайт так сделал, я копировал из IDE, и там всё было нормально.
1. Я не понимаю, к чему это? "cout « Solve_Sinc" – команда вывода, что с ней не так?
2. К сожалению, я должен реализовать решение задачи минимальными средствами, и использование double не рекомендуется.
3. Это не костыли, это просто самостоятельно реализованный математический аппарат, который хорошо работает.

[Алексей1153]

Newbie_cpp,

Сообщение от Newbie_cpp

команда вывода, что с ней не так

сейчас там открывающая кавычка, а не оператор << . Если вдруг кто-то решит поразбираться с кодом, ему придётся исправлять все синтаксические ошибки

float - это потеря точности (и скорости исполнения, кстати, тоже)

стандартные функции лучше велосипедов (всё по той же причине)

[Newbie_cpp]

Алексей1153, открывающаяся кавычка, скорее всего, тоже из-за особенностей сайта. В таком случае, вот ссылка на исходник.
Насчёт float см. выше: реализация минимальными средствами; таковы условия задачи.
Также я считаю, что математические функции реализованы мною достаточно хорошо, чтобы в условиях данной задачи заменить библиотечные.