ФедотовА.А. Численные методы интегрирования .. 2015 (А.А. Федотов, П.В. Храпов - Численные методы интегрирования, решения дифференциальных уравнений и задач оптимизации), страница 5

Поэтому,когда подобная задача встречается на практике, получают приближенное решение с помощью численных методов, в частности конечно-разностных методов.3.1. Численные методы решения задачи КошиПри изучении численных методов для задачи Коши будем считать, что она имеет единственное решение в замкнутой прямоугольной области D   a  x  b, c  u  d  . Пусть требуется найти решение задачи (3.1) на отрезке [ a, b]. Введем на отрезке [a, b] сетку hiследующим образом:34hi  {a  x0  x1  ...  xn 1  xn  b},где точки xi , i  0,1, ..., n называются узловыми точками, или узлами сетки, hi  xi  xi 1 , i  1, ..., n — шагами сетки; n — натуральное число. Если hi  h  const, то такую сетку будем обозначать h и называть равномерной. Сетку h можно задать так:h  { xi  x0  ih, x0  a , h  (b  a ) n , i  0,1, ..., n}.(3.2)В этой главе в дальнейшем будем пользоваться равномернойсеткой h с шагом h.Пусть ui  u ( xi ) — значение точного решения (3.1) в точке xi ,а yi — соответствующее приближенное значение, полученное спомощью рассматриваемого численного метода.3.1.1.

Явный метод ЭйлераПредположим, что функция f ( x, u ) в рассматриваемой области D имеет непрерывные частные производные f ( x , u ) x иf ( x , u ) u . В таком случае, как отмечалось выше, решение задачиКоши (3.1) имеет непрерывную вторую производнуюf ( x, u ) f ( x, u )u ( x) f ( x, u ).(3.3)xuЯвный метод Эйлера определяется формулами(3.4)yi 1  yi  hf ( xi , ui ), i  0, 1, ...

, n  1,где y0  u0 .Соотношения (3.4) метода Эйлера получают следующим образом. Функцию u ( x ) разлагаем по формуле Тейлора в окрестноститочки xi :h2u (i ) 2h2u (i ), u ( xi )  hf ( xi , u ( xi )) (3.5)2где i  [ xi , xi 1 ], i  0, 1, ..., n  1. Затем отбрасываем остаточныйчлен и заменяем значения u ( xi ) на yi .u ( xi 1 )  u ( xi )  hu ( xi ) 35На рис. 3.1 дана геометрическая интерпретация метода Эйлера.Изображены первые два шага метода, т. е.

проиллюстрировано вычисление значений y1 и y 2 при x  x1 и x  x2 . Интегральные кривые a0 , a1 и a2 описывают точныерешения уравнения u   f ( x, u ) сначальными условиями u ( x0 )  u0  y0 , u ( x1 )  y1 и u ( x2 )  y2соответственно. При этом криваяa0 соответствует точному решению задачи Коши (3.1), так как онапроходит через начальную точкуA0 ( x0 , u0 ). Точки A1 и А2 получены в результате численного решения задачи Коши методом Эйлера.Рис. 3.1. Геометрическая интер- Их отклонение от кривой a0 характеризует погрешность метода.претация метода ЭйлераУже при выполнении первого шагамы фактически сразу попадаем на другую интегральную кривую. Отрезок A0 A1 — отрезок касательной к кривой a0 в точке A0 .

Тангенсугла наклона касательной A0 A1 равен значению производнойu0  f ( x0 , u0 ). Касательная A1 A2 проводится уже к другой интегральной кривой a1. Таким образом, погрешность метода Эйлераприводит к тому, что на очередном шаге решение переходит на другую интегральную кривую.Пример.

Явным методом Эйлера решить задачу:u   u  x, u (0)  1.Решение. Используя (3.4), получаемy0  1;y1  y0  h ( y0  x0 )  1  h;y2  y1  h( y1  x1 )  1  h  h(1  h  h)  1  2h  2h 2 ;y3  y2  h( y2  x2 )  1  3h  6h 2  2h336и т. д. С другой стороны, точным решением является функцияu ( x)  1  x  2e x . Сравним точное решение с полученным приближенным:h 2 1 u1  u ( h)  1  h  2e h  1  h  2  1  h e 2 1  h  h 2 e1 , 0  1  h;4 h 2 2 u2  u (2h)  1  2h  2e 2 h  1  2h  2  1  2h e 2 1  2h  4h 2 e2 , 0  2  2h.Аналогичноu3  u (3h)  1  3h  9h 2 e3 , 0  3  3h.Видно, что численное решение отличается от точного на величину, содержащую члены второго порядка h 2 и выше, с коэффициентами, растущими с номером i. Ниже покажем, что ошибка вметоде Эйлера всегда не превышает значения Ch, где C — постоянная, не зависящая от h.Вообще главный вопрос для любого численного метода состоит в оценке точности приближенных значений yi .Определение 1.

Локальной ошибкой вычислений при x  xiназывается величинаi  u ( xi )  yi .Эта ошибка зависит от приближенного метода, использованного при подсчете yi , от функции f ( x, u ) и точности вычислений.Поэтому говорят, что локальная ошибка зависит от методическойошибки i( m ) (т. е. от ошибки вычислений, связанной с методомнахождения приближенного решения) и ошибки округления (напрактике все значения, полученные в результате вычислений, берутся с конечным числом знаков). Ошибку округления можноуменьшить, повышая точность арифметических вычислений, а методическая ошибка не зависит от точности вычислений, и поэтомуза счет повышения точности вычислений ее устранить нельзя.37Иными словами, методическая ошибка совпадает с локальнойошибкой для абсолютно точно найденных значений yi .Определение 2.

Глобальной методической ошибкой на отрезке[a; b] называется величина ( m )  max i( m ) .1 i  nДля использования на практике пригодны лишь те методы численного решения, для которых( m )  0 при h  0.Определение 3. Численный метод решения задачи (3.1) называется методом k-го порядка точности, если( m )  O ( h k ) при h  0 ,или, иначе,( m )  Ch k при h  0,(3.6)где C  const  0 зависит от f ( x, u ), a, b и численного методарешения задачи (3.1), но не зависит от h.Покажем, что метод Эйлера является методом первого порядкаточности.

Будем полагать, что арифметические вычисления проводятся точно, и поэтому локальная и методическая ошибки совпадают. Докажем следующее утверждение.Лемма. Пусть при любом i  0,1, ..., n  1 справедлива оценкаi 1  (1  c1h )i  c2 h k 1 ,(3.7)где h  (b  a ) n ; c1 и c2 — постоянные, не зависящие отh , c1  0, c2  0. Тогда соответствующий метод численного интегрирования задачи (3.1) является методом k-го порядка точности.Доказательство.

В соответствии с определением надо доказать оценку (3.6). Заметим, что0  u ( x0 )  y0  u ( x0 )  u0  0.Из оценки (3.7) имеемi 1  (1  c1h )i  c2 h k 1  (1  c1h )[(1  c1h )i 1  c2 h k 1 ]  c2 h k 1 38 (1  c1h )2 i 1  c2 h k 1[1  (1  c1h )]  ...  (1  c1h )i 1 0  c2 h k 1[1  (1  c1h )  ...  (1  c1h )i ]  c2 h k 1(1  c1h )i 1  1.(1  c1h )  1Заметим теперь, что(1  c1h )i 1  exp{(i  1) ln (1  c1h )}  exp{(i  1)c1h} ba(i  1)  ec1 ( ba ) .nС помощью этого последнего неравенства из предыдущей цепочки неравенств следует оценка exp c1 i 1 [c2 ( e c1 ( ba )  1)] kh ,c1выполняющаяся для всех i  0, 1, ..., n  1.

Сравнивая последнее неравенство с неравенством (3.6) и полагая C   c2 ( e c1 ( b  a )  1)  c1 ,убеждаемся в справедливости леммы.Теперь, в соответствии с леммой, чтобы доказать, что методЭйлера является методом первого порядка точности, достаточнопроверить неравенствоi 1  (1  c1h)i  c2 h k 1при всех i  0,1, ..., n  1.Вычитая (3.4) из равенств (3.5), получаемui 1  yi 1  ui  yi  h[ f ( xi , u ( xi ))  f ( xi , yi )] h2u( i ). (3.8)2По теореме Лагранжа о конечном приращении функции имеемf ( xi , u( xi ))  f ( xi , yi ) f ( xi , i )(ui  yi ),uгде точка i лежит между точками ui и yi .С помощью последнего равенства из соотношения (3.8) находим оценкуi 1  i  c1hi  c2 h 2  i (1  c1h )  c2 h 2 ,39гдеc1  max fu( x, u) , c2  0,5 max  f x( x, u )  fu( x, u ) f ( x, u ) ( x, u )D( x, u )D(см.

формулу (3.3)).Таким образом, доказано, что явный метод Эйлера имеет первый порядок точности.3.1.2. Методы Рунге — КуттыРассмотрим теперь методы, погрешность которых при стремлении h к нулю убывает с более высокой скоростью.Метод Рунге — Кутты второго порядка точности. Его расчетные формулы:k1  hf ( xi , yi ),k2  hf ( xi  h, yi  k1 ),1yi 1  yi  ( k1  k2 ), i  0, 1, ..., n  1.2(3.9)Метод Рунге — Кутты четвертого порядка точности. Вычисления с помощью этого метода проводят по формулам:k1  hf ( xi , yi ),hkk2  hf  xi  , yi  1  ,22hkk3  hf  xi  , yi  222,(3.10)k4  hf ( xi  h, yi  k3 ),yi 1  yi 1(k1  2k2  2k3  k4 ), i  0,1, ..., n  1.6Мы рассмотрели методы, для которых при вычислении yi 1нужно знать лишь значение yi , а значения приближенного решения в предшествующих точках не входят в расчетные формулы.Иными словами, одношаговые методы — это методы с «короткойпамятью».

Если «память метода» получше, то его называют мно40гошаговым. Более точно метод численного интегрирования задачиназывать l-шаговым, если при вычислении значения yi 1 используют l величин yi l 1 , yi l  2 , ..., yi yi 1.3.1.3. Многошаговые методы АдамсаИз (3.1) следует, чтоu( xi 1 )  u( xi ) xi 1xiu( x )dx xi 1f ( x, u( x ))dx xixi 1p( x )dx,xiгде p ( x ) — полином, аппроксимирующий f ( x, u ( x)) .Пусть f i  f ( xi , yi ), где yi — приближенное решение задачи(3.1), и в качестве p ( x ) возьмем интерполяционный полином,проходящий через l ранее найденных точек ( x j , f j ),( j  (i  l  1), (i  l  2), (i  l  3), ..., i ), включая текущую точку( xi , f i ).

Если l  1, то имеем явный метод Эйлера (3.4). Если l  2,то p ( x ) — линейная функция (рис. 3.2, a), проходящая через дветочки ( xi 1 , fi 1 ) и ( xi , f i ) :p( x ) ( xi  x )( x  xi 1 )fi 1 f i , p( xi 1 )  f i 1 , p( xi )  fi .hhРис. 3.2. Метод Адамса — Башфорта:а — второго порядка точности; б — третьего порядка точности41Интегрируя полином от xi до xi 1 , получаем двухшаговый метод Адамса второго порядка точности (он также называется методом Адамса — Башфорта):hyi 1  yi  (3 fi  fi 1 ).2(3.11)Если l  3, то p( x ) — парабола, проходящая через точки( xi  2 , f i  2 ), ( xi 1 , f i 1 ) и ( xi , f i ) (рис.