Метод Ромберга

В численном анализе , метод Ромберга ( Ромберга 1 955 ) используется для оценки определенного интеграла

${\displaystyle \int _{a}^{b}f(x)\,dx}$

многократно применяя экстраполяцию Ричардсона ( Richardson 1911 ) к правилу трапеции или правилу прямоугольника (правило средней точки). Оценки формируют треугольный массив . Метод Ромберга - это формула Ньютона – Котеса - он вычисляет подынтегральное выражение в равноотстоящих точках. Подынтегральная функция должна иметь непрерывные производные, хотя довольно хорошие результаты могут быть получены, если существует только несколько производных. Если можно оценить подинтегральный в разнесенных точках неодинакова, то другие методы , такие как гауссовая квадратура и Clenshaw-Curtis квадратура , как правило , более точные.

Метод назван в честь Вернера Ромберга (1909–2003), опубликовавшего метод в 1955 году.

Метод [ править ]

С использованием

${\displaystyle h_{n}={\tfrac {1}{2^{n}}}(b-a)}$

метод может быть определен индуктивно следующим образом:

${\displaystyle R(0,0)=h_{1}(f(a)+f(b))}$

${\displaystyle R(n,0)={\tfrac {1}{2}}R(n-1,0)+h_{n}\sum _{k=1}^{2^{n-1}}f(a+(2k-1)h_{n})}$

${\displaystyle R(n,m)=R(n,m-1)+{\tfrac {1}{4^{m}-1}}(R(n,m-1)-R(n-1,m-1))}$

или же

${\displaystyle R(n,m)={\tfrac {1}{4^{m}-1}}(4^{m}R(n,m-1)-R(n-1,m-1))}$

где и . В нотации большой буквы O ошибка для R ( n ,  m ) равна ( Мысовских, 2002 ) :${\displaystyle n\geq m\,}$${\displaystyle m\geq 1\,}$ harv error: no target: CITEREFMysovskikh2002 (help)

${\displaystyle O\left(h_{n}^{2m+2}\right).\,}$

Нулевая экстраполяция R ( n , 0) эквивалентна правилу трапеций с 2 ⁿ  + 1 точкой; первая экстраполяция R ( n , 1) эквивалентна правилу Симпсона с 2 ⁿ  + 1 баллами. Вторая экстраполяция, R ( n , 2), эквивалентна правилу Буля с 2 ⁿ +1 балл. Дальнейшие экстраполяции отличаются от формул Ньютона-Котеса. В частности, дальнейшие экстраполяции Ромберга очень незначительно расширяют правило Буля, изменяя веса в соотношения, подобные как в правиле Буля. Напротив, дальнейшие методы Ньютона-Котеса производят все более разные веса, что в конечном итоге приводит к большим положительным и отрицательным весам. Это свидетельствует о том, что методы Ньютона-Котеса с интерполяцией полиномов большой степени не могут сходиться для многих интегралов, в то время как интегрирование Ромберга более стабильно.

Когда вычисления функций дороги, может быть предпочтительнее заменить полиномиальную интерполяцию Ричардсона рациональной интерполяцией, предложенной Bulirsch & Stoer (1967) .

Геометрический пример [ править ]

Чтобы оценить площадь под кривой, правило трапеции применяется сначала к цельному элементу, затем к двум, затем к четырем и так далее.

Один кусочек. Обратите внимание, поскольку оно начинается и заканчивается на нуле, это приближение дает нулевую площадь.

Два куска

Из четырех частей

Восемь частей

После получения оценок по правилу трапеций применяется экстраполяция Ричардсона .

• Для первой итерации двухчастные и однократные оценки используются в формуле (4 × (более точная) - (менее точная)) / 3. Затем та же формула используется для сравнения четырехчастной и двухчастной оценки, а также аналогичным образом. для высших оценок
• Для второй итерации значения первой итерации используются в формуле (16 (более точная) - менее точная)) / 15
• Третья итерация использует следующую степень 4: (64 (более точная) - менее точная)) / 63 для значений, полученных во второй итерации.
• Схема продолжается до тех пор, пока не будет одна оценка.

• MA означает более точный, LA означает менее точный

Пример [ править ]

Например, функция Гаусса интегрируется от 0 до 1, т.е. функция ошибок erf (1) ≈ 0,842700792949715. Треугольный массив вычисляется построчно, и вычисление прекращается, если две последние записи в последней строке отличаются менее чем на 10 ⁻⁸ .

 0,77174333 0,82526296 0,84310283 0,83836778 0,84273605 0,84271160 0,84161922 0,84270304 0,84270083 0,84270066 0,84243051 0,84270093 0,84270079 0,84270079 0,84270079

Результат в правом нижнем углу треугольного массива соответствует показанным цифрам. Примечательно, что этот результат получен из менее точных приближений, полученных с помощью правила трапеции в первом столбце треугольного массива.

Реализация [ править ]

Вот пример компьютерной реализации метода Ромберга (на языке программирования Си ).

#include  <stdio.h>#include  <math.h>void dump_row ( size_t  i ,  double  * R )  {  printf ( "R [% 2zu] =" ,  i );  для  ( size_t  j  =  0 ;  j  <=  i ;  ++ j ) {  printf ( "% f" ,  R [ j ]);  }  printf ( " \ п " ); }double romberg ( double  ( * f / * функция для интегрирования * / ) ( double ),  double  / * нижний предел * /  a ,  double  / * верхний предел * /  b ,  size_t  max_steps ,  double  / * желаемая точность * /  acc )  {  двойной  R1 [ max_steps ],  R2 [ max_steps ];  // буферизует  двойной  * Rp  =  & R1 [ 0 ],  * Rc  = & R2 [ 0 ];  // Rp - предыдущая строка, Rc - текущая строка  double  h  =  ( b - a );  // размер шага  Rp [ 0 ]  =  ( f ( a )  +  f ( b )) * h * .5 ;  // первая трапециевидная ступенька dump_row ( 0 ,  Rp ); for  ( size_t  i  =  1 ;  i  <  max_steps ;  ++ i )  {  h  / =  2 .;  двойной  c  =  0 ;  size_t  ep  =  1  <<  ( i -1 );  // 2 ^ (n-1)  for  ( size_t  j  =  1 ;  j  <=  ep ;  ++ j )  {  c  + =  f ( a + ( 2* j -1 ) * h );  }  Rc [ 0 ]  =  h * c  +  .5 * Rp [ 0 ];  // R (i, 0) для  ( size_t  j  =  1 ;  j  <=  i ;  ++ j )  {  двойной  n_k  =  pow ( 4 ,  j );  Rc [ j ]  =  ( n_k * Rc [ j -1 ]  -  Rp [ j -1 ]) / ( n_k -1 );  // вычисляем R (i, j)  } // Дамп i-й строки R, R [i, i] - лучшая оценка на данный момент  dump_row ( i ,  Rc ); if  ( i  >  1  &&  fabs ( Rp [ i -1 ] - Rc [ i ])  <  acc )  {  return  Rc [ i -1 ];  } // меняем местами Rn и Rc, поскольку нам нужна только последняя строка  double  * rt  =  Rp ;  Rp  =  Rc ;  Rc  =  rt ;  }  return  Rp [ max_steps -1 ];  // возвращаем наше лучшее предположение }

________________________________________________________

Вот пример компьютерной реализации метода Ромберга (на языке программирования javascript ).

function  auto_integrator_trap_romb_hnm ( func , a , b , nmax , tol_ae , tol_rae ) // ВХОДЫ // func = интегрируемая функция // a = нижний предел интегрирования // b = верхний предел интегрирования // nmax = количество разделов, n = 2 ^ nmax // tol_ae = максимальная допустимая абсолютная приблизительная ошибка (должна быть> = 0) // tol_rae = максимальная допустимая абсолютная относительная приблизительная ошибка (должна быть> = 0) // ВЫХОДЫ // Integ_value = оценочное значение интеграла{ if  ( typeof  a  ! ==  'number' )  {  throw  new  TypeError ( '<a> должно быть числом' ); }  if  ( typeof  b  ! ==  'number' )  {  throw  new  TypeError ( '<b> должно быть числом' ); }  if  (( ! Number . isInteger ( nmax ))  ||  ( nmax < 1 )) {  throw  new  TypeError ('<nmax> должно быть целым числом, большим или равным единице.' ); } if  (( typeof  tol_ae  ! ==  'number' )  ||  ( tol_ae < 0 ))  {  throw  new  TypeError ( '<tole_ae> должно быть числом больше или равным нулю' ); } if  (( typeof  tol_rae  ! ==  'number' )  ||  ( tol_rae <= 0 ))  {  throw  new  TypeError ('<tole_ae> должно быть числом больше или равным нулю' ); } var  h = b - a // инициализируем матрицу, в которой хранятся значения интегралавар  Ромб  =  [];  // строки для  ( var  i  =  0 ;  i  <  nmax + 1 ;  i ++ )  { Romb . push ([]); для  ( var  j  =  0 ;  j  <  nmax + 1 ;  j ++ )  { Romb [ i ]. толчок ( математика . bignumber ( 0 ));  } }// вычисление значения по односегментной трапециевидной линейке Romb [ 0 ] [ 0 ] = 0,5 * h * ( func ( a ) + func ( b )) var  integ_val = Romb [ 0 ] [ 0 ]for  ( var  i = 1 ;  i <= nmax ;  i ++ ) // обновление значения с удвоением количества сегментов // с использованием только тех значений, где они должны быть вычислены // См. https://autarkaw.org / 2009/02/28 / эффективная-формула-для-автоматического интегратора на основе правила трапеции / { h = 0,5 * h var  inte = 0 для  ( var  j = 1 ;  j <= 2 ** i - 1 ;  j+ = 2 ) { var  inte = inte + func ( a + j * h ) }Romb [ i ] [ 0 ] = 0,5 * Romb [ i - 1 ] [ 0 ] + integer * h // Использование метода Ромберга для вычисления следующего экстраполируемого значения // См. Https://young.physics.ucsc.edu/115/ romberg.pdf для  ( k = 1 ;  k <= i ;  k ++ ) {  var  addterm = Romb [ i ] [ k - 1 ] - Romb [i - 1 ] [ k - 1 ] addterm = addterm / ( 4 ** k - 1.0 ) Ромб [ i ] [ k ] = Ромб [ i ] [ k - 1 ] + addterm// Расчет абсолютной приблизительной погрешности var  Ea = math . abs ( Ромб [ i ] [ k ] - Ромб [ i ] [ k - 1 ])// Расчет абсолютной относительной приблизительной погрешности var  epsa = math . абс ( Ea / Romb [ i ] [ k ]) * 100,0// Присваиваем самое последнее значение возвращаемой переменной inte_val = Romb [ i ] [ k ]// возвращаем значение, если соблюдается любой из допусков if  (( epsa < tol_rae )  ||  ( Ea < tol_ae )) { return ( integ_val ) } } } // возвращаем последнее вычисленное значение интеграла независимо от того, соблюден допуск или нет return ( inte_val ) }

Ссылки [ править ]

• Ричардсон, Л.Ф. (1911), "Приближенное арифметическое решение конечными разностями физических задач, связанных с дифференциальными уравнениями, с приложением к напряжениям в каменной дамбе", Философские труды Королевского общества A , 210 (459–470): 307 -357, DOI : 10.1098 / rsta.1911.0009 , JSTOR  90994
• Ромберг, В. (1955), "Vereinfachte numerische Integration", Det Kongelige Norske Videnskabers Selskab Forhandlinger , Тронхейм, 28 (7): 30–36.
• Тэчер младший, Генри С. (июль 1964), "Замечание об алгоритме 60: Ромберга интеграции" , коммуникаций ACM , 7 (7): 420-421, DOI : 10,1145 / 364520,364542
• Бауэр, Флорида; Rutishauser, H .; Штифель, Э. (1963), Метрополис, Северная Каролина; и другие. (ред.), «Новые аспекты в числовой квадратуре», Экспериментальная арифметика, высокоскоростные вычисления и математика, Труды симпозиумов по прикладной математике , AMS (15): 199–218
• Булирш, Роланд; Stoer, Йозеф (1967), "Справочник по серии Численное интегрирование Численное квадратурной экстраполяцией." , Numerische Mathematik , 9 : 271-278, DOI : 10.1007 / bf02162420
• Мысовских, И.П. (2002) [1994], "Метод Ромберга" , в Hazewinkel, Michiel (ed.), Encyclopedia of Mathematics , Springer-Verlag , ISBN 1-4020-0609-8
• Нажмите, WH; Теукольский, С.А. Феттерлинг, Вашингтон; Фланнери, Б.П. (2007), «Раздел 4.3. Интеграция Ромберга» , Численные рецепты: Искусство научных вычислений (3-е изд.), Нью-Йорк: Cambridge University Press, ISBN 978-0-521-88068-8

Внешние ссылки [ править ]

• ROMBINT - код для MATLAB (автор: Мартин Каценак)
• Бесплатный онлайн-инструмент интеграции с использованием методов Ромберга, Фокса – Ромберга, Гаусса – Лежандра и других численных методов.
• SciPy реализация метода Ромберга