задача такого сорта - нужно численно решить уравнение вида
интеграл, понятное дело, берется в элементарных функциях, но это просто для проверки алгоритма, изначально там стоял неберущийся. Поэтому логарифмическая часть вычисялется явно, а первую я считаю методом трапеций, с чем-то более серьезным нет желания возиться при вычислении погрешностей, но вряд ли это как-то способно изменить суть дела.
А она вот какая: начиная с определенной точности, для примера берется
, а именно с 1e-6 ответ начинает отличаться от настоящего где-то на три сотых, что в моем случае существенно. Здесь есть два момента, которые мне абсолютно не ясны. Если считать с точностью меньше, а именно с 1е-5, то ошибки почти нет. И, кроме того, разница в значениях нормально посчитанного интеграла и посчитанного программой сначала медленно падает, оставаясь на крайне низком уровне, потом же происходит резкий скачок в несколько порядков, после чего это значение начинает неумолимо набирать обороты. Буду крайне благодарен, если объясните, в чем здесь дело.
Исходный код программы (вычисление погрешности осталось от другой функции, но разница в значеняих получается вообще несущественная):
Код:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
double check(double x)
{
return x*atan(x) - 0.5*log(x*x + 1);
}
struct pare
{
double x1, x2;
};
double func(double t)
{
return 1/(t*t + 1);
}
double d2(double p)
{
double res, t;
res = 6*p*(p*p - 1);
t = sqrt(p*p + 1);
t = t*t*t*t*t;
return res/t;
}
double d2_max(double x)
{
double r1;
r1 = 0.5*sqrt(7 - sqrt(41));
if (r1 - x > nil)
return fabs(d2(x));
return fabs(d2(r1));
}
double Integral(double a, double b, double h, double* err, int flag)
{
int i=1;
double res, s;
s = 0.5*(func(a) + func(b));
while (a + i*h < b)
{
s = s + func(a + i*h);
i++;
}
res = h*s;
if (flag == 1)
s = d2_max(b)*h*h*b*b/12;
else
s = d2_max(a)*h*h*a*a/12;
*err = s;
return res;
}
struct pare solution(double T, double h, double accuracy)
{
struct pare res;
int n1, n2 = 0;
double fI, I, dI, ln, x, err;
double dif, real;
while (((double)n2)*h < sqrt(T))
n2++;
I = Integral(0, n2*h, h/10, &err, 1);
x = n2*h;
fI = x*I - 0.5*log(x*x + 1);
real = check(x);
dif = fI - dif;
printf("x = %0.15lf :: I = %0.15lf :: dif = %lf\n", x, fI, dif);
while (1)
{
dI = Integral(x, x + h, h/10, &err, 1);
I = I + dI;
n2++;
fI = x*I - 0.5*log(x*x + 1);
real = check(x);
dif = fI - real;
printf("x = %.15lf :: fI = %0.15lf :: dif = %0.15lf\n", x, fI, dif);
x = x + h;
if ( fI - err - T > 0 )
{
n1 = n2-1;
x = x - h;
while ( fI + err - T > 0 )
{
dI = Integral(x - h, x, h/10, &err, 0);
I = I - dI;
n1--;
fI = x*I - 0.5*log(x*x + 1);
printf("x = %.15lf :: fI = %0.15lf :: dif = %0.15lf\n", x, fI, dif);
x = x - h;
}
res.x1 = n1*h;
res.x2 = n2*h;
printf("x1 = %.25lf\nx2 = %.25lf\n", res.x1, res.x2);
if (h*(n2-n1) < accuracy)
return res;
printf("Not enougth accuracy\n");
system("PAUSE");
return solution(T, h/10, accuracy);
}
}
}
int main(int argc, char *argv[])
{
double T, acc;
struct pare slt;
unsigned int t1, t2;
double p1, p2;
printf("T, accuracy:\n");
scanf("%lf", &T);
scanf("%lf", &acc);
t1 = clock();
p1 = ((double)t1)/((double)CLOCKS_PER_SEC);
slt = solution(T, acc/4, acc);
t2 = clock();
p2 = ((double)t2)/((double)CLOCKS_PER_SEC);
printf("The result is approached\n");
printf("Expired time: %lf seconds \n", p2-p1);
system("PAUSE");
return 0;
}
08.05.2009, 18:21 
и другая -- для последовательного уточнения численного нахождения интеграла, задающего функцию 
. И уж конечно, надо использовать формулу не трапеций, а Симпсона -- она ненамного сложнее, но на два порядка точнее. Кроме того, погрешность надо оценивать не по явдым формулам с производными, а по правилу Рунге.
достигается где-то на 7 итерациях внешней процедуры (максимальное к-во отрезков интегрирования 1280) и, соотв., на тридцати с чем-то шагах внутренней (поскольку там всего лишь метод половинного деления). Фактически оценка поглешности почти на порядок завышена (лень было возиться с точной реализацией правила Рунге).