Непонятного рода ошибка в вычислительной задаче

[smile]

задача такого сорта - нужно численно решить уравнение вида
$ \int_{0}^{x} \frac{x - t}{t^{2} + 1}dt = \alpha$
интеграл, понятное дело, берется в элементарных функциях, но это просто для проверки алгоритма, изначально там стоял неберущийся. Поэтому логарифмическая часть вычисялется явно, а первую я считаю методом трапеций, с чем-то более серьезным нет желания возиться при вычислении погрешностей, но вряд ли это как-то способно изменить суть дела.
А она вот какая: начиная с определенной точности, для примера берется $\alpha = 1$, а именно с 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;
}

[Xaositect]

smile в сообщении #212106 писал(а):

А она вот какая: начиная с определенной точности, для примера берется $\alpha = 1$, а именно с 1e-6 ответ начинает отличаться от настоящего где-то на три сотых, что в моем случае существенно. Здесь есть два момента, которые мне абсолютно не ясны. Если считать с точностью меньше, а именно с 1е-5, то ошибки почти нет. И, кроме того, разница в значениях нормально посчитанного интеграла и посчитанного программой сначала медленно падает, оставаясь на крайне низком уровне, потом же происходит резкий скачок в несколько порядков, после чего это значение начинает неумолимо набирать обороты. Буду крайне благодарен, если объясните, в чем здесь дело.

Такое поведение ошибки очень характерно для программ вычислительной математики. Есть два источника ошибок - ошибки алгоритма(теоретические) и ошибки реализации. При стремлении шага к нулю ошибка алгоритма обычно стремится к нулю, а ошибка реализации растет.
В данном конкретном случае дело в том, что при уменьшении шага в методе трапеций увеличивается число слагаемых. Так как действительные числа в компьютере хранятся с некоторой точностью, то у каждого слагаемого есть некоторая погрешность связанная с округлением. При сложении эта погрешность может как компенсироваться, так и складываться. Если слагаемых очень много - в результате может набежать довольно большая ошибка.

В качестве решений: попробуйте метод Симпсона, ему для достижения той же теоретической точности нужно меньше слагаемых. Для достижения еще большей точности используйте длинную арифметику(библиотеки типа MP) или вычисления в рациональных числах(последнее может очень сильно замедлить работу программы)

[Watchman]

Обратите внимание: если бы вы обильно писали комментарии, то люди смогли бы разобраться с вашей программой намного легче. Соответственно, было бы больше желающих помочь.

А так 90% времени уйдет на угадывание, что какая переменная у вас означает.

[ewert]

Вот и я не смог вникнуть в логику программы -- и не только потому, что не знаю Си. Скорее потому, что логика эта совершенно перемешана.

В этой задаче должны присутствовать две вложенные итерационные процедуры: одна -- для численного решения уравнения вида $F(x)=a$ и другая -- для последовательного уточнения численного нахождения интеграла, задающего функцию $F(x)$. И уж конечно, надо использовать формулу не трапеций, а Симпсона -- она ненамного сложнее, но на два порядка точнее. Кроме того, погрешность надо оценивать не по явдым формулам с производными, а по правилу Рунге.

В программе (пардон, что на Паскале) это и реализовано, причём без каких бы то ни было оптимизаций (в частности, уравнение решается просто методом пооловинного деления). И работает вполне надёжно. Точность в $10^{-10}$ достигается где-то на 7 итерациях внешней процедуры (максимальное к-во отрезков интегрирования 1280) и, соотв., на тридцати с чем-то шагах внутренней (поскольку там всего лишь метод половинного деления). Фактически оценка поглешности почти на порядок завышена (лень было возиться с точной реализацией правила Рунге).

type  float = extended;

{----------------------------------------------------------------}
function func(x,t: float): float;
begin    func:=(x-t)/(1+t*t);     end;
{----------------------------------------------------------------}
function integr(x: float; n: longint): float;
{Интегрирование от 0 до x формулой Симсона по n отрезкам.}
var  s: float;    i: longint;
begin
  s:=func(x,0)+func(x,x);
  for i:=1 to n div 2 do   s:=s+4*func(x, x*(2*i-1)/n);
  for i:=1 to n div 2 - 1 do   s:=s+2*func(x, x*(2*i)/n);
  integr:=s*x/(3*n);    end;
{----------------------------------------------------------------}
function iter(x1,x2, a, eps: float; n: longint): float;
{Решает уравнение integr(x)=a с точностью eps.              }
{[x1,x2] -- участок локализации (предполагается корректным).}
var  x0, y0,y1,y2: float;
begin
  while x2-x1>2*eps do begin
    x0:=(x1+x2)/2;
    y0:=integr(x0, n)-a;
    y1:=integr(x1, n)-a;
    y2:=integr(x2, n)-a;
    if y0*y1<0   then x2:=x0   else x1:=x0;
    end;
  iter:=(x1+x2)/2;
  end;
{----------------------------------------------------------------}
function iter1(a: float): float;
{Решает уравнение  для "точной" функции с максимальной машинной точностью.}
var  x0,x1,x2, d,dold, y0,y1,y2: float;
begin
  x1:=0;    x2:=1;
  while x2*arctan(x2)-ln(1+x2*x2)/2<=a do  x2:=x2*2;
  d:=abs(x2-x1);    dold:=d*2+1;
  while d<dold do begin
    dold:=d;
    x0:=(x1+x2)/2;
    y0:=x0*arctan(x0)-ln(1+x0*x0)/2-a;
    y1:=x1*arctan(x1)-ln(1+x1*x1)/2-a;
    y2:=x2*arctan(x2)-ln(1+x2*x2)/2-a;
    if y0*y1<0   then x2:=x0   else x1:=x0;
    d:=abs(x2-x1);
    end;
  iter1:=(x1+x2)/2;
  end;
{----------------------------------------------------------------}
{================================================================}
var  x1,x2, a, x,xold, eps: float;
     n, nrep: longint;
{================================================================}
begin
  writeln('=========================================');

  a:=5;
  eps:=1e-10;

  n:=20;
  nrep:=0;
  xold:=-1;
  x:=0;
  while abs(x-xold)>eps do begin
    xold:=x;
    inc(nrep);
    x1:=0;    x2:=1;
    while integr(x2, n)<=a do x2:=x2*2;

    x:=iter(x1,x2, a, eps/10, n);
    n:=n*2;
    end;

  writeln('eps      = ', eps);
  writeln('a        = ', a);
  writeln('x_exact  = ', iter1(a));
  writeln('x_num    = ', x);
  writeln('err      = ', iter1(a)-x);
  writeln('last_dif = ', x-xold);
  writeln(nrep);
  writeln(n);
  readln;

  end.