вычислить параметрический несобственный интеграл

[illuminates]

Добрый день. Я в теме topic77175.html спрашивал как со стороны математики вычислить интеграл:

$$\int_{x_0}^{\infty} (e^{-4x}(32/x+32/x^2+16/x^3+4/x^4)/\sqrt{1-l^2/a^2x^2+2(1+1/x)e^{-2x}}) dx $$

Разобрался, преобразовал интеграл, разбил на два, но к сожалению, созданная программа не верно выдаёт ответ. А именно, то что выдаётся в ответе:


то что должно быть:


Если требуется, то могу подробно описать условие физической задачи и метод преобразования интеграла.

Вот код программы:

double f (double x, double l);//подкоренное выражение
double f2 (double t, double c, double l);//преобразованное подкоренное выражение, разложенное в ряд
double Fun1 (double c, double A, double l);//Метод Симпсона, когда нижний предел  особая точка
double Fun2 (double A, double B, double l);//Метод Симпсона, когда верхний предел бесконечность
double F1 (double x, double c, double l);//преобразованная функция
double F2 (double x, double l);//обычная функция
double lowLim (double a, double b, double l);//Нахождение нижнего предела интеграла

int main(void)
{   
   const short n = 10;
   double x[n], Ed[n], c, A, B, e0;
   double a, h = 0.5, x0 = 0.5; 
   double C = 3 * pow(10.0,8), H = 1.0546 * pow(10.0,-34), e = 1.6 * pow(10.0,-19), m = 9.11 * pow(10.0, -31);
   //a = pow(H, 2)/ (m * pow(e, 2));
   
   a=0.529*pow(10.0,-11);
   //e0 = pow(e,8)/(3 * a * pow(H,3) * pow(C,3));
   e0 = pow(e,5)/(3 * a * pow(H,3) * pow(C,3));
   short i;

   // для первой точки
   x[0]=x0;
   c = lowLim (0.133758, 0.133760, pow((x[0]), 2)); 

   A = c + 0.05;  
   B = A + 1; //первое приближение верхнего предела
   Ed[0] = Fun1 (c, A, pow(x[0], 2)) + Fun2 (A, B, pow(x[0], 2));

   //для последующих точек
   for (i = 1; i<n; i++)
   {
      x[i] = x[i-1] + h; 
      c = lowLim (0.64, 5.01, pow((x[i]), 2));

      A = c + 0.1;
      B = A + 1;
      Ed[i] = Fun1 (c, A, pow(x[0], 2)) + Fun2 (A, B, pow(x[i], 2));
   }
   
   ofstream out;
    out.open ("D:\\Задача 151.txt");
   for (i = 0; i<n; i++)
      //out << x[i] <<"\t" << Ed[i] <<"\n";
      out << x[i] <<"\t" << log(Ed[i]/e0) <<"\n";
   out.close();
   return 0;
}

double lowLim (double a, double b, double l)
{
   double с, E = 0.00001;
   do
   {
      с = (b + a)/ 2;
      if (f(a, l) * f(с, l) > 0) a = с;
      else b = с;
   }
   while (fabs(b-a) > E);
   return (a+b)/2;
}

double f(double x,double l)
{
   return 1-l/pow(x,2)+2*(1+1/x)*exp(-2*x);
}

double Fun1 (double c, double A, double l)
{
   short i, n = 20;
   double x, I, h;
   h = (A - c)/ n; //шаг

    I = h * (F1(c, c, l) + F1(c + n * h, c, l))/ 3; 
   for (i = 1; i<n-1; i++)//n-1!
   {
       x = c + i * h; 
      if (i%2 == 0) I += h * 2 * F1 (x, c, l)/ 3; 
      if (i%2 == 1) I += h * 4 * F1 (x, c, l)/ 3;
   }
   return I;

}

double F1 (double t, double c, double l)
{
   return (2*exp((-4)*(pow(t,2)+c))*(16+32*pow((pow(t,2)+c),-1)+32*pow((pow(t,2)+c),-2)+16*pow((pow(t,2)+c),-3)+4*pow((pow(t,2)+c),-4)))/pow(f2(t, c, l), 0.5);
      
   
}

double f2 (double t, double c, double l)
{
   return (2*(l-c*(2*pow(c,2)+2*c+1)*exp(-2*c)))/(pow(c,3))+((2*c*(2*pow(c,3)+2*pow(c,2)+2*c+1)*exp(-2*c)-3*l)*pow(t,2))/(pow(c,4))+(2*(6*l-c*(4*pow(c,4)+4*pow(c,3)+6*pow(c,2)+6*c+3)*exp(-2*c)*pow(t,4)))/(3*pow(c,5))+((exp(-2*c)*(4*pow(c,6)+4*pow(c,5)+8*pow(c,4)+12*pow(c,3)+12*(c,2)-15*exp(2*c)*l+6*c)*pow(t,6))/(3*pow(c,6)));
   
}

double Fun2 (double A, double B, double l)
{
   
   double Ed1, Ed2, /*интегралы-приближения к интегралу правого предела*/ epsilon = 0.0001;
   short i, n = 20;
   double x, h;

   Ed1 = 0;
   do
   {
      Ed2 = Ed1;
      h = (B - A)/ n;

      Ed1 = h * (F2(A, l) + F2(A + n * h, l))/ 3;
      for (i = 1; i<n; i++)
      {
         x = A + i * h; 
         if (i%2 == 0) Ed1 += h * 2 * F2 (x, l)/ 3;
         if (i%2 == 1) Ed1 += h * 4 * F2 (x, l)/ 3;
      }

      n += 10;
      B += 10*h;

   } while (fabs(Ed2 - Ed1) > epsilon);

   return Ed1;
}

double F2(double x, double l)
{
   return (exp(-4*x)*(32/x+32/(pow(x,2))+16/(pow(x,3))+4/(pow(x,4))))/(pow((1-l*(1/(pow(x,2)))+2*(1+1/x)*exp(-2*x)),0.5));
}

Уже около двух недель решаю данное задание, и был-бы очень благодарен если-бы кто то нашел ошибку!

[Утундрий]

Вы рано бросились программировать. Интеграл надо бы слегка причесать. Поковырять поведение подкоренного в $x_0$, придумать в связи с оным подходящую замену...

[illuminates]

Утундрий
Интеграл был разбит на два. Левый из которых, считался в особой точке. Для него подкоренное выражение было разложено в ряд. Далее произведена замена $x-x_0=t^2$.. Далее после небольших преобразований, особая точка была устранена. Как еще его можно причесать, совсем не знаю.

[Утундрий]

А, это в параллельной вселенной. Ну, может ошибка где-то в преобразованиях. Кстати, а на кой там разложение в ряд?

[Pavia]

Утундрий

Утундрий в сообщении #784584 писал(а):

Кстати, а на кой там разложение в ряд?

Предел интегрирования бесконечность. Как будете брать такой интеграл? Я лично знаю только 2 метода.
1) Через разложения интеграла в ряд.
2) Или свести к известной функции $erfc(x)$. Для которой такое разложение в ряд проделано и хорошо известно.

А вот зачем раскладывать знаменатель в ряд я не знаю.

PS. illuminates, в виду того что не 1 не 2 метод не используется ваше решение, то его надо детально проверять.

-- Пн ноя 04, 2013 23:06:02 --

illuminates
У Вас на графике с верным ответом $ln(\epsilon_d/\epsilon_0)$
Что это такое? Из какой это книге?

[Утундрий]

Pavia в сообщении #784759 писал(а):

Предел интегрирования бесконечность. Как будете брать такой интеграл?

Дробно-линейно отображу, например.

[Alexu007]

А почему нельзя тупо вычислить то, что написано в формуле, для х = 0,1,2,3,4,5... и так далее? До бесконечности, или пока не надоест. Единственное что - там кроме х другие переменные, про которые ничего не сказано - l, e, d. "e" вроде как константа, что такое l и d?

[illuminates]

Утундрий в сообщении #784584 писал(а):

Кстати, а на кой там разложение в ряд?

Для устранения особой точки.

Pavia в сообщении #784759 писал(а):

У Вас на графике с верным ответом $ln(\epsilon_d/\epsilon_0)$ Что это такое? Из какой это книге?

Это зависимость энергии дипольного излучения от прицельного расстояния. Взято из - Алексеев "сборник задач по классической электродинамики" стр261

Alexu007 в сообщении #784889 писал(а):

А почему нельзя тупо вычислить то, что написано в формуле, для х = 0,1,2,3,4,5... и так далее? До бесконечности, или пока не надоест. Единственное что - там кроме х другие переменные, про которые ничего не сказано - l, e, d. "e" вроде как константа, что такое l и d?

Потому, что подынтегральная функция плохая. l пробегает от 0.5 до 5. Про d я вроде ничего не говорил.