Приближение функций и численное дифференцирование. Численное интегрирование.

Лекция 5. Стр. 1 Лекция 5 Приближение функций и численное дифференцирование Численное интегрирование Численное дифференцирование применяют в тех случаях, когда невозможно или очень сложно продифференцировать функцию аналитически (будем говорить точно). Например, если функция задана таблично или вычисление производной требует большого числа арифметических операций. В этом случае чрезвычайно эффективно использование сплайн-интерполяции. Поскольку производные сплайна на каждом интервале можно найти аналитически, то приближенное значение производной функции полагают равными значению производной интерполяционного сплайна. Интерполяцию кубическим сплайном используют также для вычисления вторых производных. Кроме того, формулы численного дифференцирования широко используются при разработке численных методов решения многих задач(решение дифференциальных уравнений, поиск решений нелинейных уравнений, поиск точек экстремума функций и др.) Простейшие формулы численного дифференцирования. Вычисление первой производной. Предположим, что в окрестности точки x функция f дифференцируема достаточное число раз. Исходя из определения производной f ( x)  lim x 0 f ( x  x)  f ( x) , естественно попытаться использовать для ее вычисления две x простейшие формулы: f ( x  h)  f ( x ) h (1), соответствующих выбору фиксированных значений f ( x )  f ( x  h) f ( x)  h x  h, x   h . Разностные отношения в правых частях формул называют правой и левой f ( x)  разностными производными. Для оценки погрешностей f ( x  h)  f ( x ) h (погрешностей аппроксимации) воспользуемся формулами f ( x )  f ( x  h) r ( x, h)  f ( x)  h r ( x, h)  f ( x)  Тейлора: f ( x  h)  f ( x)  f ( x)h  f (  ) 2 h , здесь и ниже   и   - некоторые точки, расположенные 2 на интервалах ( x, x  h), ( x  h, x ) соответственно. f (  ) 2 h  f ( x) 1 2 r ( x, h)  f ( x)    f (  )h h 2 f (  ) 2 f ( x)  f ( x)  f ( x)h  h 1 2 r ( x, h)  f ( x)   f (  )h h 2 f ( x)  f ( x)h  Следовательно 1 М 2 h, M 2  max f ( ) [ x, xh] 2 1 r ( x, h)  М 2 h, M 2  max f ( ) [ x h, x ] 2 r ( x, h)  Лекция 5. Стр. 2 Таким образом формулы (1) имеют первый порядок точности по h .Приведенные формулы численного дифференцирования имеют простую геометрическую интерпретацию (см.рис.) Естественно предположить, что лучшим по сравнению с tg  , tg  приближением к f ( x)  tg является тангенс угла наклона  0 секущей к графику функции, проведенной через точки N  , N  . Соответствующая приближенная формула имеет вид: f ( x)  f ( x  h )  f ( x  h) , величину правой части формулы называют центральной 2h разностной производной. Подставляя в выражение для погрешности r ( x, h)  f ( x)  f ( x  h)  f ( x  h) 2h соответствующие разложения по формуле Тейлора f ( x) 2 f (  ) 3 h  h получим, 2 6 f (  ) 2 f (  ) 3 f (  ) 2 f (  ) 3 f ( x)  f ( x)h  h  h  f ( x)  f ( x)h  h  h 2 6 2 6 r0 ( x, h)  f ( x)  2h f (  )  f (  ) 2 r0 ( x, h)   h . Следовательно, справедлива оценка погрешности 12 1 r0 ( x, h)  М 3 h 2 , M 3  max f (3) ( ) [ x h, x  h ] . 6 f ( x  h)  f ( x)  f ( x)h  Таким образом, центральная разностная производная аппроксимирует производную f (x) со вторым порядком точности относительно h . Вычисление второй производной. Наиболее простой и широко применяемой для приближенного вычисления второй производной является формула Лекция 5. Стр. 3 f ( x)  f ( x  h)  2 f ( x )  f ( x  h) (2), величину в правой части равенства называют h2 второй разностной производной. Подставляя в выражение для погрешности r ( x, h)  f ( x)  f ( x  h)  2 f ( x )  f ( x  h) h2 соответствующие разложения по формуле Тейлора f ( x) 2 f ( x) 3 f ( 4) (  ) 4 h  h  h получим, 2 6 24 f ( 4) (  )  f ( 4) (  ) 2 r0 ( x, h)   h . Следовательно, справедлива оценка погрешности 24 1 r ( x, h)  М 4 h 2 , M 4  max f ( 4) ( ) [ x h, x  h ] 12 f ( x  h)  f ( x)  f ( x)h  Таким образом, формула (2) имеет второй порядок точности. Обусловленность формул численного дифференцирования. Несмотря на внешнюю простоту формул численного дифференцирования, их применение требует особой осторожности. Используемые при численном дифференцировании значения f * ( x ) функции f (x) непременно содержат ошибки. Поэтому к погрешности аппроксимации формул численного дифференцирования добавляется неустранимая погрешность, вызванная погрешностями вычисления функции f . Для того, чтобы погрешность аппроксимации была достаточно малой, требуется уменьшить шаг h . Однако при малых шагах формулы численного дифференцирования становятся плохо обусловленными, и результат их применения может быть полностью искажен неустранимой погрешностью. Важно понимать, что действительная причина этого явления не в несовершенстве предложенных методов вычисления производных, а в некорректности самой операции дифференцирования приближенно заданной функции. Полная погрешность реально вычисляемого значения правой разностной производной представляет собой сумму погрешности аппроксимации и неустранимой погрешности: f * ( x  h)  f * ( x )  r ( x, h)  rн ( x, h), где h f ( x  h)  f ( x ) 1 r ( x, h)  f ( x)  , а _ rн ( x, h)  ( f ( x  h)  f * ( x  h)  ( f ( x)  f * ( x))) h h Пусть  - верхняя граница абсолютной погрешности  ( f * ( x))  f ( x )  f * ( x) r * ( x, h)  f ( x)  используемых значений функции. Тогда 2 2  , число обусловленности    . Т.к.    при h  0 , то формула правой h h разностной производной при малых h становится очень плохо обусловленной. Верхняя граница полной погрешности будет неограниченно возрастать при h  0 , хотя погрешность аппроксимации rH  стремится к нулю. 1 2 M 2 h   . Ниже приведен график функции r (h ) , полученной для функции 2 h x f ( x)  e в точке x  1 с разными значениями шага h . r ( h)  Лекция 5. Стр. 4 Выберем оптимальное значение шага hопт , при котором величина r (h ) достигает минимального значения. Приравняем производную r (h ) к нулю. Получаем:  , __ rmin  r (hопт )  2 M 2 . M2 hопт  2 Формулы для вычисления производных порядка k  1 обладают еще большей чувствительностью к ошибкам задания функций. Поэтому значения производных высокого порядка, найденные помощью таких формул, могут быть очень неточными. Численное интегрирование Введение в методы численного интегрирования: постановка задачи, обусловленность, простейшие квадратурные формулы, квадратурные формулы Гаусса b Постановка задачи. Пусть требуется вычислить определенный интеграл I   f ( x)dx от a непрерывной на отрезке a, b функции f (x) . Чаще всего при приближенном вычислении определенного интеграла заменяют подынтегральную функцию некоторым обобщенным интерполяционным многочленом: n f ( x)   f  xi  i ( x)  r ( x) , где r (x ) — погрешность (остаточный член) интерполяции. После i 0 несложных вычислений получается квадратурная формула интерполяционного типа I  n  c f x  : i 0 i n n  n  I   f ( x)dx     f  xi  i  x   r ( x)  dx   f  xi   i  x dx   r  x dx   ci f  xi   R . i 0 i 0  a a  i 0 a a Здесь xi — узлы квадратурной формулы, ci — коэффициенты квадратурной формулы, R — b b b b погрешность или остаточный член квадратурной формулы. Важно понимать, что узлы и коэффициенты квадратурной формулы не зависят от подынтегральной функции, т.е. вычисленные однажды они позволяют вычислять приближенное значение интеграла для широкого класса подынтегральных функций. Простейшие квадратурные формулы Простейшие квадратурные формулы получаются из геометрических соображений. ab  — формула прямоугольников.  2  Если применить линейную интерполяцию по точкам a, f (a), b, f (b) , то получим f ( a )  f (b) I  I тр  b  a  — формулу трапеций. 2 Если f ( x)  const на a, b , то I  I пр  b  a  f  i Лекция 5. Стр. 5 ab  ab , f( ) и 2   2 Если же аппроксимировать f (x) параболой, проходящей через точки, a, f ( a ) ,  b  a   f (a)  4 f  a  b      f (b)  — формулу Симпсона.  6   2   ( x  x0 ) f ( x)  f ( x0 )  f ( x0 )( x  x0 )  f ( x0 )  ... - формула Тейлора для функции f (x) в точке 2! ab , заменим производные функции f (x) в точке x 0 .формулами разностных x 0 . Пусть x0  2 f (b)  2 f (a  b 2)  f (a ) ab f (b)  f (a) ab ) , f ( ) производных: f ( . Интерполяционный 2 2 ba 2 ba 2 b, f (b) , то получим I  I Симпс    многочлен будет иметь вид: a  b  f (b)  2 f (a  b 2)  f (a )  ab  a  b  f (b)  f (a)  P2 ( x)  f   x   x   . 2 ba 2  2   2    2ba 2 2   Проинтегрировав его получаем формулу Симпсона: f (b)  2 f (a  b 2)  f (a) 1 f (b)  f (a) ab 0  (b  a) 3 (b  a)  2 2  ba 12 2ba 2 b  P ( x)dx  f  2 a   Эти три формулы часто называют простейшими квадратурными формулами. Погрешность простейших квадратурных формул. Теорема 1.Пусть функция f (x) дважды непрерывно дифференцируема (для формулы Симпсона — четырежды) на отрезке a, b ..Тогда для квадратурных формул справедливы оценки погрешностей: M Для формулы прямоугольников Rпр  I  I пр  2 b  a 3 . 24 ab I   f ( x)dx  I пр  Rпр  (b  a) f    Rпр , используем разложение Тейлора в окрестности  2  a b подходящей точки b Rпр  I  I пр   a  ab  a  b  f ( x)dx  b  a  f      f ( x)  f    dx   2  a  2  b 2  ab ab 1 ab  a  b    a  b     f  f ' x  f ' ' () x   f       dx   2  2  2! 2   2    2   a  2 b max f ' ' ()   a  b   M a  b 1 a  b    3 3  dx  ab     f ' x  f ' ' () x  b  a   2 b  a        2   2 2 ! 2 24 24     a M4 M Совершенно аналогично получаются оценки Rтр  I  I тр  2 b  a 3 и RСимпс  I  I Симпс  b  a 5 . 12 2880 b Лекция 5. Стр. 6 b Rтр  I  I тр   a b  f(a)  f(b) b  a     f ( x)  f(a)  f(b) dx  f ( x)dx  2 2  a b 1  2  f a   f ' a x  a   2! f ' ' ( )x  a  1 1 2  f b   f ' b  x  b   a b  1  2  f ' a x  a   2! f ' ' ( )x  a  1 2 1  f ' b  x  b   a 1  2 f ' ' ( 2 ) x  2   f a   f (b)  dx  2!  M 1 2 3 f ' ' ( 2 ) x  2    2 b  a  2!  12 Составные формулы прямоугольников, трапеций и Симпсона. Для повышения точности численного интегрирования применяют составные квадратурные формулы. Разбивают отрезок на четное число n  2m одинаковых отрезков длины h  ba точками xi  a  ih , i  0, 1, ..., n и на n каждом отрезке xi , xi 1  применяют квадратурную формулу. Полученные результаты суммируют. n 1   h  h2 M 2 b  a  .  2  24 i 0 h n 1 h2    f xi   f xi 1   M 2 b  a  . 2 i 0 12 Составная формула прямоугольников: I  I пр  Rпр  h f  xi  Составная формула трапеций: I  I тр  RТр Составная формула h m1 h4         f x  4 f x  f x  M 4 b  a  .  2i 2i 1 2i  2 3 i 0 2880 Оценка погрешности квадратурных формул по правилу Рунге. Любая попытка сравнить погрешности простейших квадратурных формул упирается в необходимость оценить значения соответствующей производной подынтегральной функции. Поэтому обычно, вместо приведенных формул погрешностей квадратурных формул, применяют правило Рунге. Покажем, как получается оценка по Рунге на примере формулы прямоугольников. Вычислим по формуле прямоугольников приближенное значение интеграла I h с шагом h и h приближенное значение интеграла I h / 2 с шагом . 2 2 2 h h Имеем: I  I h  M 2 b  a  и I  I h / 2  M 2 b  a  . 24 24  4 Вычтем одно равенство из другого: h2 h2 1 Ih/2  Ih  3 M 2 b  a  , и, следовательно, M 2 b  a   I h / 2  I h  . 24  4 24  4 3 h Таким образом, погрешность квадратурной формулы прямоугольников с шагом , равна 2 I  Ih Ih/2  Ih и уточненная по Рунге формула прямоугольников имеет вид I пр  I h / 2  h / 2 . 3 3 Аналогично получаются уточненные по Рунге формулы трапеций и Симпсона: I  Ih I  Ih , I Сим пс  I h / 2  h / 2 . При этом за погрешность вычисления интеграла для I тр  I h / 2  h / 2 3 15 I  Ih I  Ih формулы трапеций полагаем равной h / 2 ,а для формулы Симпсона — равной h / 2 . 3 15 Симпсона: I  I Cиимп  RCимп  Лекция 5. Стр. 7 Оценки по Рунге позволяют строить так называемые "адаптивные" алгоритмы. Адаптивные алгоритмы состоят в следующем. Исходя из некоторого начального значения шага h, вычисляем погрешность по правилу Рунге. Если величина погрешности получается больше требуемой, делим шаг пополам и повторяем вычисления. Делим шаг пополам до тех пор, пока величина погрешности не станет меньше заданной при постановке задачи погрешности. Более тонкие адаптивные алгоритмы не только уменьшают, но и увеличивают шаг так, чтобы проводить интегрирование с максимально возможным шагом, но так, чтобы погрешность интегрирования не превысила заданную погрешность. Квадратурные формулы Гаусса Легко видеть, что квадратурная формула прямоугольников точна для многочлена нулевой степени, формула трапеций — для многочлена первой степени, а формула Симпсона — второй. Гауссом были построены квадратурные формулы наивысшего алгебраического порядка точности. 1 Т.е. в квадратурной формуле  1 n f ( x)dx   ci f ( xi ) узлы x1 , x2 , ,..., xn и коэффициенты c1 , c2 , ,..., cn i 1 подбирались так, чтобы формула была точна для всех многочленов как можно более высокой степени, (степени, превышающей n ) Доказано, что эта наивысшая степень для n узлов — 2n  1 .Как правило, сначала строят формулы Гаусса для стандартного отрезка [-1,1].Затем с помощью замены переменной t  ab ba  x 2 2 осуществляют переход к формулам интегрирования на произвольном ba n ab ba  ci f   xi  Формула отрезке:  f (t ) dt   2 i 1 2  2  a b 1  1 n f ( x)dx   ci f ( xi ) точна для i 1 многочленов степени 2n  1 тогда и только тогда, когда она точна для функций f ( x)  1, x, x 2 ,..., x 2 n1 . Это эквивалентно тому, что узлы x1 , x2 , ,..., xn и коэффициенты c1 , c2 , ,..., cn 1   1 , k  0, 1, ..., 2n  1 . k  1 i 1 1 Пример. Получим квадратурную формулу Гаусса для двух узлов, т.е. n  2, 2n  1  3 . n должны удовлетворять системе уравнений 1  ci xik   x k dx  k 1 Соответствующая квадратурная формула Гаусса имеет вид 1  1 2 f ( x)dx   ci f ( xi )  c1 f ( x1 )  c 2 f ( x 2 ) и она точна для всех многочленов до третьей степени i 1 включительно. 1 Тогда:  1  dx  c1  c2  2 , 1 1  xdx  c1 x1  c2 x2  0 , 1 1 x 1 2 dx  c1 x12  c 2 x 22  2 , 3 1  x dx  c x 3 3 1 1  c 2 x 23  0 . Получим для коэффициентов и узлов квадратурной формулы систему 1 уравнений c1  c 2  2, c x  c x  0, 2 2  1 1 2  2 2 c1 x1  c 2 x 2  3 ,  c1 x13  c 2 x 2 3  0, Лекция 5. Стр. 8 1 . Таким образом получаем квадратурную формулу 3 решение которой — c1  c2  1 и x 2   x1  1 Гаусса  f ( x)dx  f ( 1 3 1 1 ) f( 3 ) , точную для многочленов третьей степени. Замечательное свойство квадратурных формул Гаусса — возможность вычислять несобственные интегралы от неограниченных функций, поскольку узлы квадратурных формул Гаусса лежат строго 1 внутри отрезка интегрирования. Например, dx  1 x 1 2 интеграла равно  ).  2 1 1 3  6  2.5 (точное значение Обусловленность квадратурных формул интерполяционного типа На первой лекции доказано, что задача численного интегрирования хорошо обусловлена. Исследуем обусловленность задачи приближенного вычисления определенного интеграла по квадратурным формулам интерполяционного типа. Пусть при вычислении интеграла вместо точных значений f (x) функции, в квадратурных формулах фигурируют приближенные значения функции f * ( x ) . Обозначим  f * — абсолютную погрешность подынтегральной функции, f ( x)  f * ( x)    f * , и I * — абсолютную погрешность квадратурной формулы I  n  c f x  , i 1 i i n I  I *  I * . Здесь I *   ci f * xi  . i 1   f * xi     ci  f *  I * .Т.е. квадратурная формула i 1 i 1  i 1   n  устойчива к погрешностям округления, I *    ci  f *     f * , и ее число  i 1  обусловленности     ci . Поскольку все квадратурные формулы точны для f ( x)  const , то Тогда I  I *  n n  c f x    c i i n i i b  1 dx   c a i i  b  a , то     ci   ci  b  a . Получили достаточно тривиальное i i утверждение — влияние погрешности функции на приближенное значение интеграла растет с увеличением промежутка интегрирования.