Д.П. Костомаров, А.П. Фаворский - Вводные лекции по численным методам (pdf) (1113729), страница 20
Текст из файла (страница 20)
При этом результаты расчетов с шагами h2 и h3 в промежуточных точках xi ,которые не вошли в первый столбец, опущены. В последнем пятом столбце приведеныдля сравнения значения функции u ( x ) (53), дающей аналитическое решение задачи.Из таблицы видно, как по мере уменьшения шага повышается точность. В то же времяследует отметить, что даже при маленьком шаге h3 = 0.01 метод не может обеспечитьрешению хорошую точность: ошибка в последней точке x = 2 составляетz = -0.027059 .Результаты проведенных расчетов представлены также на рис. 1. На немприведены три кривые, соответствующие численному решению задачи по схемеЭйлера с шагами h1 , h2 , h3 .
При выбранном масштабе кривая III практическисовпадает с графиком аналитического решения задачи (53) (пунктирная линия).Рисунок наглядно показывает повышение точности приближенного решения по мереуменьшения шага h .Мы подробно разобрали метод Эйлера, поскольку на примере простойразностной схемы (30) он позволяет поставить и обсудить все основные вопросычисленного решения задачи Коши методом конечных разностей. Однако следуетотметить, что полученные в этом разделе результаты представляют прежде всеготеоретический интерес. Для решения реальных задач разностную схему Эйлераобычно не применяют из-за ее низкой точности: погрешность с уменьшением hубывает как O ( h ) .
В следующих разделах мы обсудим пути построения разностныхсхем более высокого порядка точности.2.2. Повышение точности разностного метода.Оценка погрешности решения через погрешность аппроксимации уравнения вметоде Эйлера (45) приводит к вполне естественному выводу: чтобы повыситьточность метода, нужно улучшить аппроксимацию дифференциального уравненияразностным. Рассмотрим возможные пути реализации этой идеи.- 97 -Предположим, что решение дифференциального уравнения u ( x ) имеетпроизводные достаточно высокого порядка и напишем для него разложение поформуле Тейлора11ui +1 = u ( xi + h ) = ui + u′ ( xi ) h + u′′ ( xi ) h 2 + u′′′ ( xi ) h 3 + L .(54)26Если его оборвать на члене порядка h и положить в соответствии сдифференциальным уравнением (27) u′ ( xi ) = f ( xi , ui ) , то мы придем к схеме Эйлера.Сделаем следующий шаг.
Оборвем разложение (54) на члене порядка h2 ивоспользуемся для вычисления производной u′′ ( xi ) формулой (46). В результатеполучим новое рекуррентное соотношение, более сложное чем (32),⎫1 ⎧ ∂f∂f(55)yi +1 = yi + f ( xi , yi ) h + ⎨ ( xi , yi ) + ( xi , yi ) f ( xi , yi ) ⎬ h 2 ,∂y2 ⎩ ∂x⎭которое можно также записать в виде разностного уравнения⎫yi +1 − yi1 ⎧ ∂f∂f(56)= f ( xi , yi ) + ⎨ ( xi , yi ) + ( xi , yi ) f ( xi , yi ) ⎬ h .∂yh2 ⎩ ∂x⎭Здесь, как и в предыдущем разделе, мы обозначили искомую функцию в разностномуравнении (56) буквой y , а не u , чтобы подчеркнуть, что (27) и (56) – это два разныхуравнения.Уравнение (55), дополненное начальным условием (31), дает явную разностнуюсхему численного решения рассматриваемой задачи Коши.
По рекуррентной формулеможно последовательно рассчитать все значения сеточной функции yi , 0 ≤ i ≤ n иполучить таким образом приближенное решение задачи (27), (28). Исследованиепоказывает, что такая усложненная схема имеет второй порядок точностиотносительно h как для аппроксимации уравнения, так и для погрешности решения.Существенно то, что основная идея данного подхода допускает дальнейшее развитие.Если оборвать разложение (54) на члене порядка h3 , h4 и т. д., то получатсяразностные схемы третьего, четвертого и более высоких порядков точности.Однако у данного подхода есть существенный недостаток. При расчетах посхеме Эйлера требуется вычислять только значения функции f ( xi , yi ) . В схеме же (55)на каждом шаге приходится вычислять не только функцию f , но и ее первые∂f∂fпроизводные( xi , yi ) и ( xi , yi ) .
Если мы, оставив в разложении (54) члены до h4∂x∂yвключительно, построим схему четвертого порядка точности, то на каждом шагепридется вычислять десять величин: функцию f ( xi , yi ) , две ее первых производных,три вторых производных и четыре третьих производных. Это существенно усложнитразработку программы и нарушит важный принцип вычислительной математики –использовать в расчетах только те величины, которые задаются условиями задачи.Формулировка задачи Коши предполагает, что известен алгоритм вычисленияфункции f ( x, u ) по значениям ее аргументов. Если этот алгоритм сводится к расчетупо простой формуле, то вычисление производных не составляет труда.
Однако- 98 -возможны и такие варианты представления алгоритма, при которых вычислениепроизводных функции f ( x, u ) либо очень сложно, либо практически невозможно.Поэтому при разработке разностных схем высокого порядка точности стремятсязаменить вычисление производных функции f ( x, u ) вычислением самой функции внескольких точках. В следующих разделах мы рассмотрим, как это удается сделать.2.3. Метод Рунге-Кутта.Рассмотрим правую часть разностного уравнения (56), содержащую первыепроизводные от функции f ( x, u ) . Главная идея метода Рунге-Кутта состоит в том,чтобы приближенно заменить ее на сумму значений функции f в двух разных точкахс точностью до членов порядка h2 .
С этой целью положим:⎫∂f1 ⎧ ∂ff ( xi , yi ) + ⎨ ( xi , yi ) + ( xi , yi ) f ( xi , yi ) ⎬ h =∂y2 ⎩ ∂x⎭= β f ( xi , yi ) + α f ( xi + γ h, yi + δ h ) + O ( h 2 ) ,(57)где α , β , γ , δ - четыре свободных параметра, которые нужно подобрать так, чтобыправая часть равнялась левой с нужной степенью точности.Разложим функцию f ( xi + γ h, yi + δ h ) по степеням h :⎧ ∂f⎫∂f(58)f ( xi + γ h, yi + δ h ) = f ( xi , yi ) + ⎨γ( xi , yi ) + δ ( xi , yi )⎬ h + O ( h 2 ) ,∂y⎩ ∂x⎭подставим разложение (58) в формулу (57) и приравняем слева и справа члены, несодержащие h и содержащие h в первой степени.
В результате получим для четырехпараметров три уравнения11α + β = 1, αγ = , αδ = f ( xi , yi ) .(59)22Они позволяют выразить параметры β , γ , δ через α :11β =1−α , γ =,δ=f ( xi , yi ) .(60)2α2αЗаменяя с помощью (57) левую часть уравнения (56) и отбрасывая членыпорядка O ( h 2 ) , получим однопараметрическое семейство разностных схем Рунге-Кутта:yi +1 − yihh⎛⎞= (1 − α ) f ( xi , yi ) + α f ⎜ xi +, yi +f ( xi , yi ) ⎟ .(61)h2α2α⎝⎠Уравнение (61), как и (30), можно записать в виде удобного для расчетоврекуррентного соотношенияhh⎡⎛⎞⎤yi +1 = yi + ⎢(1 − α ) f ( xi , yi ) + α f ⎜ xi +f ( xi , yi ) ⎟ ⎥ h ., yi +(62)2α2α⎝⎠⎣⎦- 99 -Наиболее удобные разностные схемы этого семейства соответствуют двум11значениям параметра α : α =и α = 1 .
При α =рекуррентная формула (62)22принимает видhyi +1 = yi +f ( xi , yi ) + f ( xi + h, yi + hf ( xi , yi ) ) .(63)2Она определяет следующую процедуру расчета yi +1 . Сначала делается шаг h по схемеЭйлера и вычисляется величина(64)y% i +1 = yi + f ( xi , yi ) h .Затем находится значение функции f в точке ( xi +1 , y% i +1 ) , составляется полусумма{}f ( xi , yi ) + f ( xi +1 , y%i +1 )2и проводится окончательный расчет величиныf ( xi , yi ) + f ( xi +1 , y%i +1 )(65)yi +1 = yi +h.2Такая схема вычислений называется «предиктор-корректор» или буквально«предсказание-исправление». Вычисление y% i +1 по схеме Эйлера (64) – это грубоепредсказание результата.
Вторичный расчет (65), сделанный на основании первого,является уточнением результата, его коррекцией.При α = 1 рекуррентная формула (62) имеет видhh⎛⎞yi +1 = yi + f ⎜ xi + , yi + f ( xi , yi ) ⎟ .(66)22⎝⎠hЗдесь схема расчета заключается в следующем. Сначала делается половинный шаг :2по схеме Эйлера вычисляется величинаhyi + 1 = yi + f ( xi , yi ) .(67)22()Затем находится значение функции f в точке xi + 1 , yi + 1 . Оно определяет по формуле22(66) очередное значение yi +1 .Следует заметить, что процедура расчета приближенного решения задачи Коши(27), (28) по схеме (61) по сравнению со схемой Эйлера усложняется: теперь накаждом шаге функцию f ( x, u ) приходится считать не один, а два раза.
Однако такоеусложнение оказывается оправданным благодаря более высокой точности метода. Кисследованию проблемы точности мы теперь и переходим.Введем, как и в предыдущем разделе, две сеточные функции: погрешностьрешения z (33) и погрешность аппроксимации уравнения ψ . В рассматриваемомслучае она определяется формулойu −u ⎡hh⎛⎞⎤f ( xi , ui ) ⎟ ⎥ .ψ i = i +1 i − ⎢(1 − α ) f ( xi , ui ) + α f ⎜ xi +, ui +(68)h2α2α⎝⎠⎦⎣- 100 -Выразим yi по формуле (35) через ui и zi и подставим в разностное уравнение(56). В результате получимzi +1 − zi ui +1 − uihh⎛⎞+= (1 − α ) f ( xi , ui + zi ) + α f ⎜ xi +, ui + zi +f ( xi , ui + zi ) ⎟ .
(69)hh2α2α⎝⎠Формулу (69) можно переписать в видеzi +1 − zi ⎧ ⎡hh⎛⎞⎤= ⎨ ⎢(1 − α ) f ( xi , ui + zi ) + α f ⎜ xi +, ui + zi +f ( xi , ui + zi ) ⎟ ⎥ −2α2αh⎝⎠⎦⎩⎣⎡− ⎢(1 − α ) f ( xi , ui ) + α f⎣hh⎛⎞⎤ ⎫, ui +f ( xi , ui ) ⎟ ⎥ ⎬ −⎜ xi +2α2α⎝⎠⎦ ⎭(70)⎧u − u ⎡hh⎛⎞⎤ ⎫− ⎨ i +1 i − ⎢(1 − α ) f ( xi , ui ) + α f ⎜ xi +f ( xi , ui ) ⎟ ⎥ ⎬ ., ui +2α2α⎝⎠⎦ ⎭⎣⎩ hЗдесь мы перенесли член ( ui +1 − ui ) / h слева направо и в каждое из двух выражений,собранных в фигурных скобках, добавили одно и то же слагаемое.
Поскольку междуфигурными скобками стоит знак минус, значение правой части формулы (70) в целомпри этом не меняется. Однако благодаря таким преобразованиям мы собрали вовторых фигурных скобках члены, которые дают погрешность аппроксимациидифференциального уравнения ψ (68).Перейдем к дальнейшему исследованию соотношения (70). Рассмотримфункциюhh⎛⎞F ( v ) = (1 − α ) f ( xi , v ) + α f ⎜ xi +,v +f ( xi , v ) ⎟ .(71)2α2α⎝⎠Выражение, стоящее в первых фигурных скобках формулы (70), можно записать какразность значений этой функции при v = ui + zi и v = ui и преобразовать эту разность спомощью формулы конечных приращений Лагранжа(72)F ( ui + zi ) − F ( ui ) = F ′ ( ui + θ i zi ) zi , 0 < θ i < 1 ,где∂f∂f ⎛hhh ∂f⎞⎛F ′ ( v ) = (1 − α ) ( xi , v ) + α ⎜ xi +,v +f ( xi , v ) ⎟ ⎜ 1 +( xi , v ) ⎞⎟ .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
