Численные методы решения задачи Коши для системы обыкновенных дифференциальных уравнений (метод Зонневельда)

Размещено на http://www.allbest.ru/

РЕФЕРАТ

Пояснительная записка - 30 страниц, источников - 5, рис. 12

Выполнение курсовой работы из курса «Численные методы» по теме «Численные методы решения задачи Коши для системы обыкновенных дифференциальных уравнений (метод Зонневельда)» имеет такую цель и задание:

Улучшить теоретические знания из курса «Численные методы» по теме «Численные методы решения задачи Коши для системы обыкновенных дифференциальных уравнений (метод Зонневельда)»;

Развивать навыки самостоятельной творческой работы;

Приобретать навык самостоятельной работы с литературными источниками;

Уметь применить знание из численных методов во время развязывания практических заданий.

Ключевые слова: СХОДИМОСТЬ, ПОРЯДОК, ПОГРЕШНОСТЬ, НЕВЯЗКА, АЛГОРИТМ, АППРОКСИМАЦИЯ, ИНТЕГРИРОВАНИЕ, НОРМА, КОЭФФИЦИЕНТ, НЕПРЕРЫВНОСТЬ.

ВВЕДЕНИЕ

Выполнение курсовой работы из курса «Численные методы» по теме «Численные методы решения задачи Коши для системы обыкновенных дифференциальных уравнений (метод Зонневельда)» имеет такую цель:

Улучшить теоретические знания из курса «Численные методы» по теме «Численные методы решения задачи Коши для системы обыкновенных дифференциальных уравнений (метод Зонневельда)»;

Состоит из двух частей: теоретической и практической. В теоретической части приводятся основные определения и теоремы, необходимые для полноценной реализации на практике. Во второй части приведены решения модель Лефевра-Пригожина - брюсселятор с использованием метода Зонневельда.

1. ТЕОРЕТИЧЕСКАЯ ЧАСТЬ

1.1 Постановка задачи Коши для ОДУ и систем ОДУ. Классификация методов решения

Рассматриваем задачу Коши для системы обыкновенных дифференциальных уравнений

, (1)

или, подробнее,

(2)

. (3)

Хорошо известны условия, гарантирующие существование и единственность решения задачи Коши. Предположим, что функция непрерывны по всем аргументам в замкнутой области

Из непрерывности функции следует их ограниченность, т. е. существование константы такой, что всюду в выполняются неравенства

Предположим, кроме того, что в функции удовлетворяют условию Липшеца по аргументам , т. е.

,

для любых точек и области .

Если выполнено вышесказанное, то существует единственное решение Системы (2), определенное при и принимающее при заданные начальные значения (3).

Методы решения задачи (1), (2) можно разделить на три класса: точные аналитические, приближенные аналитические, методы численные. Точные аналитические методы изучаются в дисциплине «Дифференциальные уравнения». Они позволяют получить решение в виде формулы, однако круг их применения достаточно узок. В приближенных аналитических методах строятся последовательности функций, сходящиеся (обычно равномерно) к точному решению (1), (2). Основная сложность при их реализации - это вычисление большого числа производных и интегралов, среди которых могут появиться «неберущиеся». Численные методы носят наиболее универсальный характер. При их реализации решения получается в виде массива чисел, являющих приближенным значением решения на некоторой системе точек.

Итак, численное решение задачи (1), (2) ищется в узлах сетки

,

где . Обозначаем , , расстояние между соседними узлами сетки. Если , то сетка называется равномерной (регулярной); в противном случае - неравномерной (нерегулярной). В случае равномерной сетки , ; для неравномерной сетки , где - параметр нерегулярности.

На выходе численного метода мы получаем последовательность значений , являющихся приближениями к значениям точного решения в узлах сетки . Набор чисел называется каркасом приближенного решении задачи Коши. По этим точкам, используя аппарат теории аппроксимации, можно построить непрерывную функцию, являющуюся приближенным решением задачи (1), (2). Набор чисел

.

Называют проекцией точного решения задачи Коши (1), (2) на сетку . Значение , определяется в узле , может вычисляться явно

,(1.3)

или неявно

,(1.4)

т.е. для нахождения нужно решить нелинейное уравнение.

Численные методы делятся на одношаговые () и m-шаговые (многошаговые). В одношаговых методах для получения точки (,) требуется лишь о последней рассчитанной точки (,). В m-шаговых методах для получении точки (,) требуется информация о предыдущих m расчетных точках.

1.2 Явные методы типа Рунге-Кутты (общие понятия, сходимость)

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

Методы Рунге-Кутта базируются на использовании сокращенного ряда Тейлора. Разложим точное решение задачи

,(4)

в кругу узла сетки в ряд Тейлора и ограничим ряд членами p-го порядка:

.

Основная идея построения явных методов Рунге-Кутта p-го порядка состоит в замене части ряда Тейлора функцией , которая содержит значение правой части уравнения в промежуточных точках отрезка и приближает с точностью до p-го порядка где С - константа, которая не зависит от шага .

В общем случае методы Рунге-Кутта записываются в виде:

(5)

где

,(6)

- константы, имеют такой рекурсивный вид:

(7)

Методы обозначенные формулами (5)-(7), называются явными методами Рунге-Кутта, поскольку в них следующие приближенные значения находят прямым расчетом по формулам (6) и (7). Порядок точности методов Рунге-Кутта зависит от количества членов в отрезках (6) и значений констант .

Обычно коэффициенты или, короче, .

Эти условия были приняты Куттой без каких-либо комментариев. Смысл их в том, что все точки, в которых вычисляется , являются приближениями первого порядка к решению. Эти условия сильно упрощают вывод условий, определяющих порядок аппроксимации для методов высокого порядка. Однако для методов низких порядков эти предложения не являются необходимыми.

Метод Рунге-Кутты (13) имеет порядок , если для достаточно гладких задач (4) , т.е. если ряды Тейлора для точного решения и для совпадают до члена включительно.

После статьи Бутчера (1964г) вошло в обычай символически представлять метод (13) посредством следующей таблицы:

Можно отметить, что методы Рунге-Кутта при p>5 не используются.

Остановимся более подробно на отдельных методах. При p=1 получаем схему Эйлера:

локальная погрешность

,

при p=2 получаем семейство методов

(8)

Исследуем погрешность аппроксимации методов (8) в зависимости от выбора параметров. Исключая из последнего уравнения функции k₁ и k₂, получаем

(9)

По определению погрешностью аппроксимации или невязкой метода (8) называется выражение

(10)

Полученное заменой в (9) приближенного решения u_i точным решением y_i=y(x_i).

Найдем порядок погрешности аппроксимации в предложении достаточной гладкости решения y(x) и функции f(x,y). Для этого разложим все величины, входящие в выражение (10), по формуле Тейлора в точке x_i Имеем

,

,

где . Далее, согласно уравнению (1), получим

,

поэтому

.

Отсюда видно что методы (8) имеют первый порядок аппроксимации, если

Если же дополнительно потребовать то получим методы второго порядка аппроксимации. Таким образом, имеется однопараметрическое семейство двухэтапных методов Рунге-Кутта второго порядка аппроксимации. Это семейство методов можно записать в виде

,(11)

где

При получаем метод второго порядка:

Двухэтапных методов третьего порядка аппроксимации не существует. Чтобы убедиться в этом, достаточно рассмотреть уравнение . Для него двухэтапный метод Рунге-Кутта (11) принимает вид

,

и погрешность аппроксимации равна

.

Разлагая по формуле Тейлора и учитывая, что , получим

Отсюда видно, что наивысший достижимый порядок аппроксимации равен двум.

При использовании численных методов основным является вопрос о сходимости. Рассмотрим сходимость при .

Фиксируем точку x и строим последовательность сеток таких, что и (тогда необходимо ). Говорят что численный метод сходится в точке x, если

при .

Метод сходится на отрезке , если он сходится на каждой точке x этого отрезка.

Рассмотрим основные виды погрешностей, возникающих при численном решении задачи (1), (2).

Может оказаться, что начальное значение известно не точно, а определяется в результате эксперимента, например, с помощью измерений или в результате решения какой-либо другой задачи. В этом случае вместо задачи (1), (2) будет решаться задача



с измененным начальным условием, причем . Решение этой задачи зависит от и не совпадает с решением задачи (1), (2). Разность

,

называется не устойчивой погрешностью решения .

Разность

,

называется локальной погрешностью или погрешность метода. Здесь -значение точного решения при , а -значение приближенного решения, получаемого по формуле (1.3) или (1.4) при условии, что вместо приближенного значения используются значения, соответствующие точному решению, т.е. значения .

Глобальной ошибкой численного метода называется величина

где -значение, получаемое по формулам (1.3) или (1.4) при k=N-1.

Глобальная ошибка определяется:

а) ошибками округления и ошибками арифметических действий, обусловленными числом разрядов ЭВМ и характером выполняемых операций для расчета ;

б) методическими ошибками, определяемыми выбранным алгоритмом;

в) переходными ошибками, обусловленными тем, что при расчете вместо точных значений используются приближенные значения , полученные на предыдущих шагах;

г) ошибками в задании начальных данных.

Численный метод имеет p-й порядок точности, если существует число p>0 такое, что

Можно показать, что если локальная погрешность имеет порядок p+1, т.е. при , то глобальная погрешность имеет на единицу меньший порядок, т.е. при .

На практике в качестве характеристики метода часто используют величину

Полная (накопленная) погрешность приближенного решения в точке равна сумме неустранимой погрешности, погрешности метода и вычислительной погрешности. Можно показать, что для одношаговых методов

,

где -погрешность задания начального условия, - локальная погрешность метода, - погрешность округления при переходе от точки к точке заключено между и -интегральная кривая, проходящая через точку (),- фактическое получение значения по формуле одношагового метода .

Как видно из формулы для ,

1) если , т.е. , то полная погрешность равна сумме локальных погрешностей;

2) если (т.е. интегральные кривые расходятся), то влияние локальных погрешностей, полученных на предыдущих шагах, возрастает и полная погрешность больше суммы локальных погрешностей;

3) если (т.е. интегральные кривые сближаются), то влияние локальных погрешностей ослабевает и полная погрешность меньше суммы локальных погрешностей.

1.3 Практические способы оценки погрешности приближенного решения

1.3.1 Правило Рунге

Непосредственное вычисление интеграла

.(1)

Характеризующего главный член погрешности, довольно обременительно.

В практической работе обычно идут по несколько иному пути.

Зададимся некоторыми M₁ и M₂ , M₁ < M₂ и произведем дважды численное интегрирование с разбиением отрезка интегрирования на части и соответственно. В результате получим приближенные значения интеграла и .

,(2)

где . Предположим . Если бы в этих соотношениях отсутствовали величины и , то получилось бы два уравнения с двумя неизвестными и , которые можно было бы определить; в противном случае можно получить лишь некоторые приближения к ним. Исключая , получаем соотношение для определения :

.

Отсюда .

Таким образом, погрешность квадратуры может быть представлена в виде

,

первое слагаемое .

Можно вычислить. Если оно является главным членом погрешности, то можно считать, что , и использовать величину как оценку погрешности. Оценка погрешности через величину носит название правила Рунге. Чтобы было главным членом погрешности , имеющей порядок , необходимо выполнение условия

Если, например, ~, где , то

,

и требуемое условие выполняется.

Может случиться, что величина , определяемая равенством (1), равна нулю или очень мала. Тогда при реально используемых М слагаемое не будет существенно больше по модулю, чем , и проводимые нами построения будут выглядеть как довольно сомнительные. В этом случае при оценке погрешности через величину зачастую будет создаваться превратное представление о величине погрешности.

Если является главным членом величины и , мы имеем право написать равенство

,

т.е. величина дает лучшее по порядку приближение к исходному значению интервала.

Существуют различные процедуры вычисления интервалов, опирающиеся на использование правила Рунге. В одних случаях задаются некоторыми и последовательно вычисляют при после каждого вычисления величины , , определяют соответствующую величину

,

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

Против этой методики часто выдвигается такое возражение: если значение при различных вычисляется независимо, то происходит большой перерасход работы вследствие вычисления большого числа «ненужных» значений .

На самом деле, этот перерасход не является столь катастрофическим, как это кажется на первый взгляд. Объем работы при вычислении каждого значения пропорционален ~ общий объем работы пропорционален

~.

Если относить к «полезной» работе только работу по вычислению ,то доля «полезной»работы будет равна отношению к , т.е. . Если относить к «полезной» работе также работу по вычислению значения необходимо для оценки погрешности, то доля полезной работы равна отношению к , т.е. .

Другой путь использования правила Рунге заключается в следующем. Производим интегрирование с относительно небольшим числом узлов М₁ и М₂ и затем определяем величину . В одномерном случае обычно берут . Затем из уравнения определяем требуемое число и осуществляем интегрирование с отрезками разбиения. затем контролируем точность результата, применяя правило Рунге (если оно является законным) по результатам расчетов с и узлами. Если при данных и применение этого правила является сомнительным, то с целью получения представления о погрешности производим расчет еще при каком-либо числе отрезков разбиения . В случае если точность оказывается неудовлетворительной, по результатам и пересчитываем новое значение и затем новое требуемое число отрезков разбиения . При выборе оптимальной стратегии этого алгоритма надо иметь в виду следующее: может быть, стоит завышать значение по сравнению с решением уравнения , чтобы избежать лишних пересчетов.

Сформулируем правило выбора критерия практической оценки погрешности. Пусть - решение задачи, - получаемые приближения, . За погрешность приближенного решения принимается некоторое выражение , удовлетворяющее условию при . Очевидно, ~.

В рассматриваемом случае задачи численного интегрирования за меру погрешности, удовлетворяющую этому условию, мы можем принять величину если только . Однако соотношение ~ является лишь асимптотическим, а применяется при конкретных М₁ и М₂; поэтому в особо ответственных случаях прибегают к способам проверки его реальной выполнимости. Например, в описанном выше случае, когда последовательно вычисляют величины при , имеем

.

Если число узлов настолько велико , что член существенно меньше, чем , то

.

По близости величины к можно судить о возможности ограничиться при оценке погрешности величиной .

1.3.2 Комбинация специально подобранных формул (метод вложенных форм)

Идея их состоит в том, чтобы вместо пользования экстраполяцией Ричардсона построить такие формулы РК, которые сами содержали бы кроме численного приближенного значения некоторое выражение более высокого порядка (или каким-либо иным образом более точное, чем ). Последнее могло бы тогда служить для управления погрешностью и длинной шага на каждом шаге. В частности, это удешевило бы выбраковку шагов.

Итак, нам надо найти такую таблицу коэффициентов

Чтобы величина

,(1)

имела порядок р, а

,(2)

- порядок q (обычно или ).

Согласно теореме 2.13 (Чтобы метод Рунге-Кутты (13) имел порядок р, необходимо и достаточно выполнения равенств для всех деревьев t порядка, меньшего или равного р), мы должны удовлетворить условиям

 - для всех деревьев порядка ,(3)

и

- для всех деревьев порядка .(4)

Первые методы такого типа предложил Мерсон (1957), Ческино (1962) и Зонневельд (1963) (см. табл. 4.1, 4.2 и 4.3, в которых заголовок вида «фамилия p (q)» указывает, что порядок равен р, а порядок «оценщика погрешности» равен q)

Однако эти методы еще не полностью удовлетворяют сформулированным нами выше требованиям. У Мерсона имеет порядок 5 только для линейных уравнений с постоянными коэффициентами. Поэтому этот метод переоценивает погрешность при малых h, но все же работает вполне хорошо и применяется весьма широко. Точно так же у Зонневельда оценивает не погрешность аппроксимации, а «последний учетный член». Метод Ческино не экономичен, потому что он оценивает погрешность со слишком высокой точностью.

1.4 Интегрирование с переменным шагом. Автоматический выбор шага интегрирования

До сих пор мы рассматривали такой процесс решения задач

,

когда формула Рунге-Кутта применялась с одной и той же величиной шага интегрирования во всей области вычисления решения вне всякой зависимости от характера поведения решения. Это был либо шаг h, либо h/2 в случае повторного счета для оценки погрешности по правилу Рунге

Либо

Какая бы из этих величин не использовалась, она не изменялась от точки к точке. В этом случае говорят, что решение задачи получено с постоянным шагом интегрирования. Применение переменного шага интегрирования позволяет учитывать характер поведения решения и уменьшить общее число шагов, сохранив при этом требуемую точность приближенного решения. Тем самым могут быть снижены объем работ и машинное время и замедлен рост вычислительной погрешности.

Имея в распоряжении

Оценки локальной погрешности метода, величину шага интегрирования можно выбирать автоматически в процессе счета. При этом можно исходить из того, что на каждый шаг приходилась приблизительно одинаковая погрешность.

2. ПРАКТИЧЕСКАЯ ЧАСТЬ

2.1 Содержательная и формальная постановка задачи

Содержательная постановка задачи.

Рассмотрим подробно модель Лефевера и Николиса (1971), так называемый «брюсселятор»: предположим, что шесть веществ, A, B, D, E, X, Y, участвуют в следующих реакциях:

(бимолекулярная реакция)(1)

(автокаталитическая трехмолекулярная реакция)

.

Если обозначить через A(x), B(x),… концентрации вещества A, B,…как функций времени x, то реакции (1), согласно закону действующих масс, описываются следующими дифференциальными уравнениями:

Формальная постановка задачи.

Теперь упростим эту систему: исключим из рассмотрения уравнения для D и E, так как они не влияют на остальные; предположим, что A и B поддерживаются постоянными (положительными) и возьмем все скорости реакций равными единице. Далее, введем обозначения , и в результате получим систему

2.2 Качественный анализ брюсселятора

Полученная система имеет одну особую точку при .

В окрестности этой точки линеаризованное уравнение неустойчиво тогда и только тогда, когда . Далее изучая область, где , или положительны или отрицательны, мы придем к выводу, что все решения остаются ограниченными.

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

Результаты исследования поведения системы при A=1, B=0,1,2,3 (зависимости компонент системы от времени и фазовые траектории) указаны на рис. 2.1 - 2.4 а)-в).

а)б)в)

Рис.2.1 - Поведение системы при B=0 (а - график функции , б - график функции , в - фазовый портрет)

а)б)в)

Рис.2.2 - Поведение системы при B=1 (а - график функции , б - график функции , в - фазовый портрет)

а)б)в)

Рис.2.3 - Поведение системы при B=2 (а - график функции , б - график функции , в - фазовый портрет)

а)б)в)

Рис.2.4 - Поведение системы при B=3 (а - график функции , б - график функции , в - фазовый портрет)

2.3 Применение метода Зонневельда для расчета брюсселятора

Расчетные формулы метода Зонневельда для данной системы уравнений имеют вид:



где .

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

2.4 Результаты вычислительного эксперимента

Приближенное решение



Приближенное решение

Фазовый портрет

График изменения шага интегрирования

График оценки локальной погрешности



График оценки локальной погрешности

ЗАКЛЮЧЕНИЕ

Во время выполнения курсовой работы из курса «Численные методы » по теме «Численные методы решения задачи Коши для системы обыкновенных дифференциальных уравнений (метод Зонневельда)» я:

- Улучшила теоретические знания из курса «Численные методы» по теме «Численные методы решения задачи Коши для системы обыкновенных дифференциальных уравнений (метод Зонневельда)»;

- Развила навыки самостоятельной творческой работы;

- Улучшила навыки применения знаний из численных методов во время решения практических заданий.

ПЕРЕЧЕНЬ ССЫЛОК

1. Бабушка И., Витасек Э., Прагер М. Численные процессы решения дифференциальных уравнений. - М.: Мир, 1969. - 368 с.

2. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. - М.: Лаборатория Базовых Знаний, 2001. - 632 с.

3. Хайрер Э., Нерсет С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. - М.: Мир. 1990. - 512 с.

4. Фельдман Л.П., Петренко А.I., Дмитрієва О. А. Чисельні методи в інформатиці. - К.: Видавнича група BHV, 2006ю - 480 с.

5. Самарский А. А., Гулин А. В. Численные методы: М.: Наука 1989 - 432 с.

ПРИЛОЖЕНИЕ

задача погрешность интегрирование брюсселятор зонневельд



Размещено на Allbest.ru