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

[aitap]

Добрый день!

Необходимо по экспериментальным данным подобрать параметры $A_0$ (нормировочный множитель, порядка 1e6), $d$ (сдвиг линии, порядка 1e-1..1e1), $\omega$ (ширина линии, того же порядка), $\alpha$ (параметр ионной восприимчивости, обычно между 0 и 1), $b_0$ (базовая линия, порядка 1e-1) в функции

$$\frac{j(\lambda)}{A_0} = b_0 + \frac{2}{\pi^2} \int\limits_0^\infty d \beta \int\limits_0^\infty dx \frac {\beta x \sin(\beta x) e^{-x^{\frac{3}{2}}} } { 1+ \left(\frac{(\lambda - \lambda_0)-d}{\omega} - \alpha^{\frac{4}{3}}\beta^2\right)^2 } $$

$\lambda_0$ - константа, которую можно получить в базе данных. Базовую линию $b_0$, в принципе, я могу вычесть и сам, если это будет нужно.

Экспериментальные данные - спектры, снятые в 1390 точках, сам пик может занимать от нескольких десятков точек до всего снимка. Функция была набита в MATLAB:

Используется синтаксис Matlab M

% does *not* accept vectors
function [val] = j_Holtsmark(lambda,baseline,lambda_0,d,width,alpha,A_0)
        z = (lambda-lambda_0-d)./width;
        val = baseline+2/pi^2.*A_0.*integral2( ...
                @(beta,x) ...
                        beta .* x .* sin(beta.*x) .* exp(-x.^(3/2)) ...
                        ./ (1 + (z-alpha.^(4/3).*beta.^2).^2), ...
                0, Inf, ...
                0, Inf, ...
                'RelTol', 1e-2 ... % даже с такой точностью очень медленно
        );
end

% usage: feed to lsqcurvefit
function [val] = j_Holtsmark_vector(param,wavelengths)
        val = arrayfun(@(x) ...
                ( ...
                        j_Holtsmark(x,param(1),param(2),param(3),param(4),param(5),param(6)) ...
                ), ...
                wavelengths ...
        );
end
 

Проблема: если не указывать желаемую точность, а воспользоваться значениями по-умолчанию, часто получается NaN и предупреждение вида "The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.". Даже если снизить относительную точность до 1e-2, расчёт одной суммы квадратов ошибок для 1390 точек занимает несколько сотен секунд (>150 в MATLAB, >600 в Octave 3.8.2; процессор Intel Core 2 Duo E4600; в данный момент ПО 32-битное), а lsqcurvefit не заканчивается за ночь. Смена интегратора в Octave между quadgk и quadcc даёт только выбор между "ещё медленнее" и "больше предупреждений об ошибках интегрирования".

По-видимому, либо я неадекватно использую свои инструменты, либо они для такого не предназначены.

Спектров сейчас 14, но будет больше, а интеграл, возможно, понадобится сменить на более сложный (о нём напишу, когда и если точно разберусь, что он из себя представляет). Подскажите, пожалуйста, направление для поиска более быстрых способов аппроксимировать эту функцию (по-видимому, нужно убыстрить расчёт интеграла). Я знаю немного C, и у меня есть книга по Fortran, если это понадобится.

[aitap]

Во-первых, избавился от параметров $A_0$ и $b_0$: вычет базовой линии и нормировку лучше сделать заранее, чем тратить время на подбор ещё двух параметров, которые мало помогут мне в работе. Во-вторых, объединил параметры $\lambda$, $\lambda_0$, $d$, $\omega$ в один безразмерный $z$, почти как это было записано в книге (Griem, Spectral Lines Broadening by Plasmas: Appendix IV.b, самое начало и Section II.2a, формула (38)).

В-третьих, переписал интеграл на Фортране (gfortran-4.9) с использованием порта QUADPACK на Fortran-90 (с расчётом воспользоваться в дальнейшем библиотекой MINPACK для запуска нелинейного метода меньших квадратов):

Используется синтаксис Fortran

module constants
implicit none
  real, parameter :: epsabs = 0.0E+00
  real, parameter :: epsrel = 0.001E+00
  real, parameter :: PI = 3.141592653589793E+00
end module

module inner_int_param
implicit none
  real beta_param
end module

module outer_int_param
implicit none
  real z_param
  real alpha_param
end module

function j_holtsmark (lambda, lambda_0, d, omega, alpha)
  use outer_int_param
  use constants
  implicit none

  real :: j_holtsmark, lambda, lambda_0, d, omega, alpha
  real, external :: outer_integrand
  real :: abserr, result
  integer :: ier = 0
  integer neval

  ! precalculate params which will be used once per integration
  z_param = (lambda - lambda_0 - d)/omega 
  alpha_param = alpha**(4./3.)

  call qagi( &
   outer_integrand, &
   0, & ! bound
   1, & ! (bound,+inf)
   epsabs, &
   epsrel, &
   result, &
   abserr, &
   neval, &
   ier &
  )

  j_holtsmark = result * 2/PI**2
  return
end

function inner_integrand(x)
  use inner_int_param
  implicit none
  real :: inner_integrand, x

  inner_integrand = x * sin(beta_param*x) * exp(-(x**1.5))

! FIXME: the function is effectively 0 when x>10, which causes IEEE_UNDERFLOW for big x (~1e2)
! FIXME: ocsillations make the function very hard to integrate for big beta (~1e2)

  return
end

function outer_integrand(beta)
  use outer_int_param ! read them
  use inner_int_param ! fill them
  use constants
  implicit none

  real :: outer_integrand, beta

  real, external :: inner_integrand
  real :: result, abserr
  integer :: neval, ier

  beta_param = beta

  call qagi( &
   inner_integrand, &
   0, & ! bound 
   1, & ! (bound, +inf)
   epsabs, &
   epsrel, &
   result, & ! to be used
   abserr, & ! estimated error
   neval, & ! number of calls
   ier & ! error code, should be 0
  )
  ! FIXME: result can be obtained without integration? 

  outer_integrand = beta * result / ( &
    1 + (z_param - alpha_param * beta**2 )**2 &
  )
  return
end

program main
  implicit none
  real, external :: j_holtsmark
  real :: storage

   write ( *, '(g14.6)' ) j_holtsmark(435.33, 435.6, 0.09, 0.05, 0.95)

  stop
end
 

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

Во-первых, когда qagi вызывает inner_integrand для больших $x$ (например, 233), она либо получает SIGFPE, если исключения в операциях с плавающей точкой включены (сборка с -ffpe-trap=underflow), либо при выходе показывает "Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG" (если они выключены). Учитывая, что $\beta x \sin(\beta x) e^{-x^{\frac{3}{2}}}$ весьма близка к 0 (много меньше своего пика и много меньше интересных мне значений) уже при $x>10$, я бы предпочёл, чтобы я просто получал 0. Как правильно обрабатывать этот случай? Можно ли просто собирать без -ffpe-trap и игнорировать это сообщение?

Во-вторых, интегратор qagi при расчёте $\int\limits_0^\infty \beta x \sin(\beta x) e^{-x^{\frac{3}{2}}} dx$ при больших $\beta$ (например, 233) возвращает код ошибки 2 "the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved", а для ещё больших (например, 935) - 1 ("maximum number of subdivisions allowed has been achieved") (и $|\mathrm{abserr}| > |\mathrm{result}|$, на порядок). Все эти безобразия, видимо, из-за того, что функция очень осциллирует:

Как это можно обойти?

Согласно Wolfram Alpha, аналитическое решение существует:

$\int\limits_0^\infty \beta x \sin(\beta x) e^{-x^{\frac{3}{2}}}dx = $ $\beta \frac{2}{729} \Bigg($ $ 243 \beta {}_3\mathrm{F}_4 \left( \frac{3}{4}, 1, \frac{5}{4}; \frac{1}{3}, \frac{2}{3}, \frac{5}{6}, \frac{7}{6}; -\frac{4 \beta^6}{729} \right)$ $+ 22 \beta^5 \Gamma(\frac{2}{3}) {}_2\mathrm{F}_3 \left( \frac{17}{12}, \frac{23}{12}; \frac{4}{3}, \frac{3}{2}, \frac{11}{6}; -\frac{4 \beta^6}{729}} \right) + $ $28 \beta^3 \Gamma(-\frac{2}{3}) {}_2\mathrm{F}_3 \left( \frac{13}{23}, \frac{19}{12}; \frac{2}{3}, \frac{7}{6}, \frac{3}{2}; -\frac{4 \beta^6}{729} \right)\Bigg)$

Где $\Gamma$ - гамма-функция, а ${}_p\mathrm{F}_q$ - гипергеометрическая функция.

Гамма-функцию я точно вычислю (она должна быть в Fortran 2008), для гипергеометрической код в интернете тоже находится. Попробую покопать в этом направлении.

[aitap]

Попытки рассчитать $ {}_p \mathrm{F} _q $ для больших по модулю аргументов $z$ (куда входит $\beta^6$) при помощи библиотеки PFQ не увенчались успехом ("ERROR - VALUE OF EXPONENT REQUIRED FOR SUMMATION WAS LARGER THAN THE MAXIMUM MACHINE EXPONENT" и "ERROR IN IPREMAX--DID NOT FIND MAXIMUM EXPONENT" в разных ситуациях); я пошёл на сделку с совестью и проигнорировал возвращаемый из внутреннего интеграла $ \int \limits_0^\infty \beta x \sin(\beta x) e^{-\frac{3}{2}} dx$ код ошибки 2 (даже при $\beta \approx 0$, как ни странно), поскольку несмотря на это, график полученного интеграла от $\lambda$ выглядит точь-в-точь, как экспериментальные спектры.

При помощи Maxima взял частные производные от подинтегрального выражения по параметрам и написал для библиотеки MINPACK минимизируемую функцию ($ j_{\text{теор}} - j_{\text{эксп}} $ по всем точкам, в которых снят спектр) и якобиан по параметрам:

Используется синтаксис Fortran

module holtsmark_errf_param
 implicit none
 ! minpack fits real*8
 real (kind=8) :: lambda_0
 real (kind=8), dimension(:,:), pointer :: spectrum
end module

module j_deriv_d
 use holtsmark_errf_param
 use constants
 implicit none
 real (kind=8) :: z_param, alpha_param
!$omp threadprivate(z_param, alpha_param)
 contains

 function d_j_d_d(lambda, d, omega, alpha, A_0)
  ! this is nonsense because integrands are real*4
  ! but minpack fits real*8
  real (kind=8) :: d_j_d_d, omega, alpha, lambda, d, A_0
  ! quadpack integrates real*4
  real (kind=4) :: integral, abserr
  integer (kind=4) :: neval, ier


  z_param = (lambda - d - lambda_0)/omega
  alpha_param = alpha**(4./3.)

  call qagi ( &
   d_j_d_d_outer_int, &
   0.0_4, 1_4, &
   epsabs, epsrel, &
   integral, &
   abserr, neval, ier &
  )

  d_j_d_d = A_0 * 4 / (pi**2 * omega) * integral
  return
 end function

 function d_j_d_d_outer_int(beta)
  use holtsmark, only: beta_inner => beta_param, &
                       holtsmark_inner => inner_integrand
  ! quadpack integrates real*4
  real (kind=4) :: integral, abserr, d_j_d_d_outer_int, beta
  integer (kind=4) :: neval, ier


  beta_inner = beta

  call qagi( &
   holtsmark_inner, &
   0.0_4, 1_4, &
   epsabs, epsrel, &
   integral, &
   abserr, neval, ier &
  )

  d_j_d_d_outer_int = integral &
   * beta * (real(z_param,4)-beta**2.*real(alpha_param,4)) &
   / ( ((real(z_param,4)-beta**2.*real(alpha_param,4))**2.+1)**2. )
  return
 end function
end module

module j_deriv_omega
 use holtsmark_errf_param
 use constants
 implicit none
 real (kind=8) :: z_param, alpha_param
!$omp threadprivate(z_param, alpha_param)
 contains
 function d_j_d_omega(lambda, d, omega, alpha, A_0)
  ! this is nonsense because integrands are real*4
  ! but minpack fits real*8
  real (kind=8) :: d_j_d_omega, omega, alpha, lambda, d, A_0
  ! quadpack integrates real*4
  real (kind=4) :: integral, abserr
  integer (kind=4) :: neval, ier
  z_param = (lambda - d - lambda_0)/omega
  alpha_param = alpha**(4./3.)
  
  call qagi( &
   d_j_d_omega_outer_int, &
   0.0_4, 1_4, &
   epsabs, epsrel, &
   integral, abserr, &
   neval, ier &
  )

  d_j_d_omega = A_0 * 4 * z_param * integral / (pi**2 * omega)
  return
 end function

 function d_j_d_omega_outer_int(beta)
  use holtsmark, only: beta_inner => beta_param, &
                       holtsmark_inner => inner_integrand
  ! quadpack integrates real*4
  real (kind=4) :: integral, abserr, d_j_d_omega_outer_int, beta
  integer (kind=4) :: neval, ier
  beta_inner = beta

  call qagi( &
   holtsmark_inner, &
   0.0_4, 1_4, &
   epsabs, epsrel, &
   integral, &
   abserr, neval, ier &
  )

  d_j_d_omega_outer_int = integral * beta &
   * (real(z_param,4) - beta**2 * real(alpha_param,4)) &
   / ( 1 + (real(z_param,4) - beta**2 * real(alpha_param,4))**2 )**2
  return
 end function
end module

module j_deriv_alpha
 use holtsmark_errf_param
 use constants
 implicit none
 real (kind=8) :: z_param, alpha_param
!$omp threadprivate(z_param,alpha_param)
 contains
 function d_j_d_alpha(lambda, d, omega, alpha, A_0)
  ! this is nonsense because integrands are real*4
  ! but minpack fits real*8
  real (kind=8) :: d_j_d_alpha, omega, alpha, lambda, d, A_0
  ! quadpack integrates real*4
  real (kind=4) :: integral, abserr
  integer (kind=4) :: neval, ier
  z_param = (lambda - d - lambda_0)/omega
  alpha_param = alpha**(4./3.)
  
  call qagi( &
   d_j_d_alpha_outer_int, &
   0.0_4, 1_4, &
   epsabs, epsrel, &
   integral, abserr, &
   neval, ier &
  )

  d_j_d_alpha = A_0 * 16./(3. * PI**2) * alpha**(1./3.) * integral
  return
 end function

 function d_j_d_alpha_outer_int(beta)
  use holtsmark, only: beta_inner => beta_param, &
                       holtsmark_inner => inner_integrand
  ! quadpack integrates real*4
  real (kind=4) :: integral, abserr, d_j_d_alpha_outer_int, beta
  integer (kind=4) :: neval, ier
  beta_inner = beta

  call qagi( &
   holtsmark_inner, &
   0.0_4, 1_4, &
   epsabs, epsrel, &
   integral, &
   abserr, neval, ier &
  )

  d_j_d_alpha_outer_int = integral * beta**3 &
   * (real(z_param,4) - beta**2 * real(alpha_param,4)) &
   / ( 1 + (real(z_param,4) - beta**2 * real(alpha_param,4))**2 )**2

  return
 end function
end module

module holtsmark_errf
 !use, intrinsic :: ieee_arithmetic, only: ieee_value !missing until gfortran 5, sigh
 use holtsmark_errf_param
 implicit none
 contains
 subroutine error_function(m, n, x, fvec, fjac, ldfjac, iflag)
   use j_deriv_d
   use j_deriv_omega
   use j_deriv_alpha
   use holtsmark
   implicit none
   integer (kind = 4) ldfjac 
   integer (kind = 4) n, m ! parameters, points
   real (kind=8) fjac(ldfjac,n)
   real (kind=8) fvec(m)
   integer (kind = 4) iflag ! 1 => functions->fvec; 2 => jacobian->fjac
   real (kind=8) x(n)
   
   integer (kind=4) index
 
   print*, 'x:', x(1:4)
   if (iflag == 1) then
     if ((x(2)<0).or.(x(3)<0).or.(x(4)<0)) then
     ! omega, alpha, A0 should not be <0; let's not waste time on that
       fvec(1:m) = sqrt(-real(n))!ieee_value(fvec,ieee_quiet_nan)
       return
     end if
!$omp parallel do shared(fvec,spectrum) private(index) schedule(guided)
     do index=1,m
      fvec(index) = j_holtsmark(spectrum(index,1),lambda_0,x(1),x(2),x(3),x(4))-spectrum(index,2)
      !print*, 'err(',index,',',spectrum(index,1),'):', fvec(index)
     end do
!$omp end parallel do
     print*, 'sumsq(err):', sum(fvec**2) 
   else if (iflag == 2) then
     if ((x(2)<0).or.(x(3)<0).or.(x(4)<0)) then
     ! omega, alpha, A0 should not be <0; let's not waste time on that
       fjac(1:ldfjac,1:n) = sqrt(-real(n))!ieee_value(fvec,ieee_quiet_nan)
       return
     end if
!$omp parallel do shared(fjac,spectrum,x) private(index) schedule(dynamic)
     do index=1,m
!$omp parallel sections
!$omp section
       fjac(index,1) = d_j_d_d(spectrum(index,1),x(1),x(2),x(3),x(4)) ! d
!$omp section
       fjac(index,2) = d_j_d_omega(spectrum(index,1),x(1),x(2),x(3),x(4)) ! omega
!$omp section
       fjac(index,3) = d_j_d_alpha(spectrum(index,1),x(1),x(2),x(3),x(4)) ! alpha
!$omp section
       fjac(index,4) = j_holtsmark(spectrum(index,1),lambda_0,x(1),x(2),x(3),x(4))/x(4) ! A_0
!$omp end parallel sections
       !print*, 'J(',index,',',spectrum(index,1),'):', fjac(index,:)
     end do
!$omp end parallel do
     print*, 'sumsq(J):', sum(fjac**2)
   end if

   return
 end subroutine
end module
 

(процедуру error_function я передаю в LMDER1 библиотеки MINPACK после считывания спектра и начальных значений в память)

После добавления OpenMP и правильной расстановки атрибута threadprivate по переменным модулей фит даже стал иногда заканчиваться за ночь. Больше всего мне помог совет доброго человека равномерно отобрать из 1390 экспериметнальных точек $N_{\text{параметров}} \cdot 20 = 80$ штук и предварительно фитить только по ним; это позволило закончить за ночь все расчёты.

Зато оказалось, что базовая линия имеет вид наклонной, а не горизонтальной прямой, так что её придётся вернуть в аппроксимацию в виде слагаемого $b_0 + b_1 \lambda$, благо, частные производные по новым параметрам будут проще некуда, $1$ и $\lambda$ соответственно. Параметр $A_0$ пришлось тоже вернуть, потому что моя теоретическая функция нормирована на $\int\limits_{-\infty}^{\infty}j(z)dz=1$, а с экспериментальными данными такое с достаточной точностью сделать не удаётся.

На данный момент хороший фит (единственный, окончания которого по 1390 точкам я дождался) выглядит вот так:

а плохой (нуждающийся в нормальной поправке на базовую линию) - так:


Ещё один путь, которым я попробую ускорить расчёты - табулировать $\mathrm{H}(\beta) = \int \limits_0^\infty \beta x \sin(\beta x) e^{-\frac{3}{2}} dx$ через определённый интервал $\Delta \beta$ (пока не решил, какой); вне пределов табулирования возвращать 0, а внутри - линейно интерполировать.

[Евгений Машеров]

А знаменатель из-под внутреннего интеграла не выносится?

-- 18 май 2016, 06:39 --

А то, что остаётся, нельзя как-то через ТФКП, скажем, получить?

[Евгений Машеров]

Вольфрам выдал что-то страшненькое через гамма- и гипергеметрические функции...
http://www.wolframalpha.com/input/?i=in ... b+x+sin(bx)+exp(-x%5E(3%2F2))+from+0+to+infinity

-- 18 май 2016, 09:31 --

integrate b x sin(bx)exp(-x^(3/2)) from 0 to infinity

[Vince Diesel]

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

$$ f(x)=\int_0^{\infty } \frac{\beta \sin (\beta x)}{\left(c-a \beta ^2\right)^2+1} \, d\beta. $$

Для $a,c,x>0$ математика дает ответ

$$ \frac{i \pi \left(e^{-x\sqrt{\frac{-c+i}{a}}}-e^{i x \sqrt{\frac{c+i}{a}}}\right)}{4 a}= \frac{\pi}{2a} e^{-\frac{\sqrt[4]{c^2+1} x \sin \left(\frac{1}{2} \cot ^{-1}(c)\right)}{\sqrt{a}}} \sin \left(\frac{\sqrt[4]{c^2+1} x \cos \left(\frac{1}{2} \cot ^{-1}(c)\right)}{\sqrt{a}}\right).$$

Здесь $\cot ^{-1}$ - арккотангенс. Внешний интеграл уже несложно посчитать численно.

[aitap]

Евгений Машеров в сообщении #1124279 писал(а):

А знаменатель из-под внутреннего интеграла не выносится?

Выносится. В расчётных программах вынесен, чтобы не делать лишнего. Изначально был записан с объединением, потому что я пользовался функцией MATLAB'а integral2, которая принимает одно большое подынтегральное выражение и две пары пределов.

Евгений Машеров в сообщении #1124290 писал(а):

Вольфрам выдал что-то страшненькое через гамма- и гипергеметрические функции...

Ага, я приводил это решение. Мне не удалось воспользоваться им лучше, чем численным интегрированием.

Vince Diesel в сообщении #1124309 писал(а):

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

$$ f(x)=\int_0^{\infty } \frac{\beta \sin (\beta x)}{\left(c-a \beta ^2\right)^2+1} \, d\beta. $$

Ага! Производные под интеграл вношу, когда якобиан считаю, а порядок интегрирования попробовать поменять не догадался. Большое спасибо!

Vince Diesel в сообщении #1124309 писал(а):

Для $a,c,x>0$ математика дает ответ

$$ \frac{i \pi \left(e^{-x\sqrt{\frac{-c+i}{a}}}-e^{i x \sqrt{\frac{c+i}{a}}}\right)}{4 a}= \frac{\pi}{2a} e^{-\frac{\sqrt[4]{c^2+1} x \sin \left(\frac{1}{2} \cot ^{-1}(c)\right)}{\sqrt{a}}} \sin \left(\frac{\sqrt[4]{c^2+1} x \cos \left(\frac{1}{2} \cot ^{-1}(c)\right)}{\sqrt{a}}\right).$$

Здесь $\cot ^{-1}$ - арккотангенс. Внешний интеграл уже несложно посчитать численно.

Единственное - у меня ограничение $c > 0$ не выполняется в половине случаев. Постараюсь разобраться, как был посчитан такой интеграл и что менятся при $ -10 \le c \le 10 $.

$$ \int\limits_0^\infty \frac{ \beta \sin(\beta x) d\beta }{1 + (c - a \beta^2)^2} = \frac{1}{x^2}\int\limits_0^\infty \frac{z \sin z dz}{1+c^2 - 2\frac{c a z^2}{x^2} + \frac{a^2 z^4}{x^4}}$$

Табличный интеграл $ \int\limits_0^\infty \frac{x^m \sin(x)dx}{\mathbf{P}_n(x)}$ я не нашёл. WolframAlpha возвращает Ваш ответ, но даёт условия $ \Im \left( \frac{ \sqrt{c+i} x}{\sqrt{a}} \right) > 0 \wedge \Im \left( \frac{\sqrt{c-i}x}{\sqrt{a}} \right) > 0 $. Насколько я понимаю комплексную арифметику, при положительных действительных $x$,$a$ и действительном $c$ второе условие не выполняется? По крайней мере, с моими значениями $ \Im \left( \frac{\sqrt{2-i} * 1}{\sqrt{0,25}} \right) < 0 $.

Кроме того, теперь стоит попытаться найти частные производные интеграла по параметрам,

$$ \frac{2}{\pi^2} \frac{\partial}{\partial \alpha}\int\limits_0^\infty dx \int\limits_0^\infty d \beta \frac{ \beta x \sin(\beta x) e^{-x^{\frac{3}{2}}} }{ 1+ \left(\frac{(\lambda - \lambda_0)-d}{\omega} - \alpha^{\frac{4}{3}}\beta^2\right)^2 } = \frac{16 \alpha^\frac{1}{3}}{3\pi^2} \int\limits_0^\infty x e^{-x^\frac{3}{2}} dx \int\limits_0^\infty d \beta \frac{ \beta^3 \sin(\beta x) \left( \frac{\lambda-d-\lambda_0}{\omega} - \beta^2 \alpha^\frac{4}{3} \right) }{ \left( 1+ \left( \frac{\lambda-d-\lambda_0}{\omega} - \beta^2\alpha^\frac{4}{3} \right)^2 \right)^2 } $$

$$ \frac{2}{\pi^2} \frac{\partial}{\partial d}\int\limits_0^\infty dx \int\limits_0^\infty d \beta \frac{ \beta x \sin(\beta x) e^{-x^{\frac{3}{2}}} }{ 1+ \left(\frac{(\lambda - \lambda_0)-d}{\omega} - \alpha^{\frac{4}{3}}\beta^2\right)^2 } = \frac{4}{\pi^2 \omega}\int\limits_0^\infty x e^{-x^\frac{3}{2}} dx \int\limits_0^\infty d \beta \frac{ \beta \sin(\beta x) \left( \frac{\lambda-d-\lambda_0}{\omega} - \beta^2 \alpha^\frac{4}{3} \right) }{ \left( 1+ \left( \frac{\lambda-d-\lambda_0}{\omega} - \beta^2\alpha^\frac{4}{3} \right)^2 \right)^2 } $$

$$ \frac{2}{\pi^2} \frac{\partial}{\partial \omega}\int\limits_0^\infty dx \int\limits_0^\infty d \beta \frac{ \beta x \sin(\beta x) e^{-x^{\frac{3}{2}}} }{ 1+ \left(\frac{(\lambda - \lambda_0)-d}{\omega} - \alpha^{\frac{4}{3}}\beta^2\right)^2 } = \frac{4(\lambda-d-\lambda_0)}{\pi^2 \omega^2}\int\limits_0^\infty x e^{-x^\frac{3}{2}} dx \int\limits_0^\infty d \beta \frac{ \beta \sin(\beta x) \left( \frac{\lambda-d-\lambda_0}{\omega} - \beta^2 \alpha^\frac{4}{3} \right) }{ \left( 1+ \left( \frac{\lambda-d-\lambda_0}{\omega} - \beta^2\alpha^\frac{4}{3} \right)^2 \right)^2 } $$

Что упирается в поиск

$$\int\limits_0^\infty d x \frac{ x^m \sin(x) \left( a - b x^2 \right) }{ \left( 1+ \left( a - b x^2 \right)^2 \right)^2 },\, m = 1,3 $$

Буду продолжать перечитывать Демидовича и напишу, если куда-то продвинусь.

[Vince Diesel]

aitap в сообщении #1124485 писал(а):

Насколько я понимаю комплексную арифметику, при положительных действительных $x$,$a$ и действительном $c$ второе условие не выполняется?

Да. Видимо, для такого количества параметров какие-то внутренние проверки не проходят. Однако для конкретных отрицательных значений $c$ математика считать не отказывается. И для них получается та же формула. В действительном виде упрощается однако до

$$ \frac{\pi}{2a} e^{-\frac{\sqrt[4]{c^2+1} x \sin \left(\frac{1}{2} \cot ^{-1}(c)\right)}{\sqrt{a}}} \sin \left(\frac{\sqrt[4]{c^2+1} x \sin \left(\frac{1}{2} \cot ^{-1}(-c)\right)}{\sqrt{a}}\right).$$

aitap в сообщении #1124485 писал(а):

Что упирается в поиск

$$\int\limits_0^\infty d x \frac{ x^m \sin(x) \left( a - b x^2 \right) }{ \left( 1+ \left( a - b x^2 \right)^2 \right)^2 },\, m = 1,3 $$

Буду продолжать перечитывать Демидовича и напишу, если куда-то продвинусь.

Тут нужен не Демидович, а ТФКП и вычисление интегралов с помощью вычетов.

Для $a>0,b>0$ и $m=1$ математика дает

$$ \frac{\pi \left(-\frac{e^{i \sqrt{\frac{a+i}{b}}}}{\sqrt{a+i}}-\frac{i e^{-\sqrt{\frac{-a+i}{b}}}}{\sqrt{-a+i}}\right)}{16 b^{3/2}}= -\frac{\pi \cos \left(\sqrt\frac{\sqrt{a^2+1}+a}{2b}-\frac{1}{2} \cot ^{-1}(a)\right) e^{-\frac{\sqrt[4]{a^2+1} \sin \left(\frac{1}{2} \cot ^{-1}(a)\right)}{\sqrt{b}}}}{8\sqrt[4]{a^2+1} {b^{3/2}}}, $$

а для $m=3$

$$\frac{i \pi e^{-\sqrt{\frac{-a+i}{b}}} \left(2 \sqrt{b} \left(-1+e^{\frac{\sqrt{-a+i}+i \sqrt{a+i}}{\sqrt{b}}}\right)+i \sqrt{a+i} e^{\frac{\sqrt{-a+i}+i \sqrt{a+i}}{\sqrt{b}}}+\sqrt{-a+i}\right)}{16 b^{5/2}}= $$

$$ =-\frac{\pi e^{-\frac{\sqrt[4]{a^2+1} \sin \left(\frac{1}{2} \cot ^{-1}(a)\right)}{\sqrt{b}}} \left(2 b \sin \left(\sqrt\frac{\sqrt{a^2+1}+a}{2b}\right)+\sqrt[4]{a^2+1} \sqrt{b} \cos \left(\sqrt\frac{\sqrt{a^2+1}+a}{2b}+\frac{1}{2} \cot ^{-1}(a)\right)\right)}{8b^3} $$

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

[aitap]

Большое спасибо!

Результаты численного интегрирования аналитического внутреннего интеграла сходятся с результатами численного двойного интегрирования, и фит одного спектра занимает около 30 секунд. Больше ни с какими проблемами не столкнулся, кроме опечаток при набивании формул. О том, как доставать после Левенберга-Маркуардта погрешность определения полученных параметров, нашёл в книге William H. Press et al. - Numerical Recipes.

Думаю, проблема решена.