Лекции - Решение обыкновенных дифференциальных уравнений

Подождите немного. Документ загружается.

Формулы Адамса являются сильно устойчивыми (наиболее

устойчивыми из всех многошаговых методов).

Формулы Нистрема и Милна являются слабо устойчивыми.

Ихиспользование может привести к искажению результата на

больших интервалах интегрирования.

7.Реализация неявных разностных схем.

Отыскание решения разностного уравнения сводится к решению

алгебраического или трансцендентного уравнения. Часто это

уравнение решается методом итераций.

Итерационный процесс вычисления y

n+k

по формуле (18)

можетбыть представлен в виде:

   

 

  













ininiiniknknkkn

yxfhyyxfhy





, (95)

0



, 1, 2,…,

или

     

 

  













knknknknkknkn

yxfyxfhyy

, (96)

0



, 1, 2,…,

Начальное приближение y

n+k

(0)

может быть выбрано различными

способами. Например, можно положить y

n+k

(0)

= y

n+k-1

. Вкачестве

начального приближения можно взять значение, полученное по

явной разностной формуле. Втаком случае явная формула

называется предсказывающей (или предиктором), неявная

формула—исправляющей (или корректором), а весь

комбинированный процесс называется предсказывающе-

исправляющим методом, или методом прогноза и коррекции (или

методом предиктор-корректор).

Количество итераций



, которое необходимо выполнить, зависит

от величины шага интегрирования, степеней явной и неявной

формул, требуемой точности решения неявного уравнения. Если

степень предсказывающей формулы меньше степени исправляющей

формулы, то каждая вновь выполняемая итерация (95) или (96)

увеличивает порядок точности очередного приближения на единицу

до тех пор, пока не будет достигнут порядок исправляющей

формулы. Для того чтобы точность приближенного решения

 





определялась степенью неявной формулы, число итераций должно

быть не меньше разности степеней неявной и явной формул.

Если интерполяционный метод Адамса использовать в разностной

форме (70), то итерационный процесс строится с помощью

соотношения:

   









inn

fhyy





. (101)

Вэтом случае для начала итерационного процесса необходимо

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

пусть

 

ff 



и тогда по формулам составления конечных разностей:

fff

1 







ki 

1k

,…,1, (103)

можно вычислить все разности, входящие в (101).

Итерационный процесс повторяется до тех пор, пока в пределах

заданной точности не установятся приближения

 





. Накаждой

итерации переопределяются конечные разности по формуле:

fff











ki 

1k

,…,1. (104)

Формулу (101) можно преобразовать к другому виду.

Воспользовавшись формулами (73), (74) и (104) можно выразить:

  

 























nin

fff

0 0



 

 

 











ink

iink

fffff

0 0



Подставив полученное выражение в (70), можно получить

1 









inn

fhfhfhyy



. (105)

Обозначая P

n+1

приближение к y

n+1

, полученное по явной формуле

Адамса (53):









inn

fhyP



формулу (105) можно представить так:

11 









knn

fhfhPy



. (106)

Тогда итерационный процесс принимает вид

   

 















nnk

knn

yxfhfhPy

. (107)

Формула (107) более удобна для счета по сравнению со (101),т.к.

здесь не требуется пересчитывать конечные разности на каждой

итерации. Конечные разности вычисляются только один раз после

окончания итерационного процесса.

8.Построение начальных значений.

Разностная формула (2), предназначенная для нахождения

решения y

n+k

в точке x

n+k

, содержит значения решения y

, y

n+1

,..., y

n+k-

в предыдущих узлах x

, x

n+1

,..., x

n+k-1

и соответствующие значения

правой части f(y

, x

), f(y

n+1

, x

n+1

),..., f(y

n+k-1

, x

n+k-1

). Предполагается, что

указанные величины заранее известны. Совокупность этих величин

составляет так называемый фронт многошагового метода.

Поэтому применение многошаговых методов требует

предварительного построения недостающих начальных значений или

конечных разностей. Недостающие начальные значения могут быть

найдены с помощью одношаговых методов, например, методом

Рунге-Кутта. Построение фронта многошагового метода в начале

интегрирования иногда называют разгоном.

11.Практические способы оценки погрешности приближенного

решения.

Если разностная формула Адамса используется в форме конечных

разностей, то оценка погрешности равна

 

111







knn

hfhyxy



. (144)

Отсюда получается асимптотическая оценка

 

knn

fhyxy 





111



, (145)

которая представляет собой первый отброшенный член разностной

формулы (53).

Т.к. в формулах Адамса (53) и (70) коэффициенты



' имеют

противоположные знаки



> 0,



'< 0, i > 0, то при использовании

метода прогноза и коррекции, основанного на явной и неявной

формулах Адамса, предсказанное y

n+1

(0)

и исправленное

 





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

точному решению:

 

 

   

1111

,max





nnnnnn

yyyxyyxy



. (146)

Кроме (146) могут быть использованы более точные оценки:

     

 

111







knnkn

hyy





, (147)

     

 

1111

0/'







knnkn

hyy





. (148)

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

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

интегрирования для многошаговых методов имеет существенной

значение, однако использование конечно-разностных схем вносит

свою особенность в этот процесс.

12.2.5Переход к произвольному шагу.5

Пусть h



h новая длина шага интегрирования. Для метода

Адамса в конечно-разностной форме необходимо пересчитать

конечные разности. Если ограничится членами до четвертого

порядка включительно,то

        

 

432

321

' h

fff

nnn



















    

 

4232

222

7/1117

' h















 

333

' h











, (164)

 

5444

0' hff





4.ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЖЕСТКИХ СИСТЕМ.

(Введение)

(Чего не могут явные методы)

[набрано из разных разделов]

1.10.3.5Распространение ошибок в начальных данных при решении

обыкновенных дифференциальных уравнений.

Вкачестве примера рассматривается численное решение задачи

Коши для обыкновенного дифференциального уравнения

yy 



 

0 yy 

решением которой является функция y(x) = y

5exp(



x), для анализа

распространения ошибок при использовании метода Эйлера.

Приближенное начальное значение y

отличается от точного

начального значения y(x

) либо из-за погрешностей в задании

начальных данных, либо из-за ошибок округления в представлении

значения y(x

) в ЭВМ, либо по обеим причинам сразу. Поэтому в

реальном начальном значении y

содержится ошибка e

5= y(x

) - y

Необходимо также отметить, что применение метода Эйлера в

данном случае эквивалентно использованию выражения 15+ h



для

аппроксимации функции exp(



x).

Выражение для глобальной ошибки E

n+1

на (n+1)-ом шаге

интегрирования запишетсятак:

 



 111 nnn

yxyE

 

 

 

   

111exp ehxyhhn







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

Первая компонента представляет собой ошибку, вытекающую из

аппроксимации значения exp(



x) выражением 15+ h



, вторая

компонента отражает распространение ошибки e

в начальном

условии. Если 1+ h



> 1, то вторая компонента растет при

увеличении n и в конечном счете станет доминирующей частью

глобальной ошибки E

n+1

. Если



< 0, то, для того чтобы ограничить

распространение ошибки e

, необходимо обеспечить выполнение

условия

11 



или



/2h

Данное условие на величину шага интегрирования h называется

условием устойчивости, поскольку в случае его нарушения

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

образом, условие устойчивости накладывает ограничение на

максимальную величину шага интегрирования h

max

, значение

которого зависит как от параметров задачи, так и от применяемого

метода. Условие устойчивости h < 2/



получено для метода Эйлера.

Для методов Адамса (3.41’) это условие записывается в виде











и приводит к таким ограничениям на шаг интегрирования:

для метода второго порядка аппроксимации h



<1,

для метода третьего порядка аппроксимации h



< 0.545,

для метода четвертого порядка аппроксимации h



<0.3,

для метода пятого порядка аппроксимации h



< 0.163.

Для классического метода Рунге-Кутта (2.22) четвертого порядка

аппроксимации h



< 2.785.

Для явного метода необходимое условие асимптотической

устойчивости принимает вид:

consth



max/

. (4.8)

Если матрица системы дифференциальных уравнений имеет

большие по модулю отрицательные собственные значения



, то

ограничение (4.8) на шаг h является на больших интервалах

интегрирования слишком обременительным,т.к. оно требует

использования очень малого шага на протяжении всего интервала

интегрирования. Если матрица системы имеет большой разброс

собственных значений, например, в случае двух уравнений







и



< 0,



< 0, то первая составляющая точного решения

     

222111

expexp exCexCxy 



где e

, e

— собственные вектора матрицы, очень быстро затухает на

отрезке, длина которого имеет порядок



= 1/



, а затем становится

ничтожно малой. Именно на этом отрезке она вносит свой вклад в

решение y(x). Вторая составляющая заметно изменяется на отрезке

длины 1/



, причем длина второго промежутка значительно

превосходит длину первого. Навтором отрезке именно вторая

составляющая определяет поведение решения.

Таким образом, на интервале интегрирования выделяются два

промежутка с разным характером поведения решения. Напервом

промежутке, называемом пограничным слоем, для воспроизведения

быстро изменяющегося решения с приемлемой точностью необходим

шаг интегрирования, удовлетворяющий условию h 1/



. Навтором

промежутке, где решение характеризуется малой скоростью

изменения, шаг нельзя увеличить,т.к. при нарушении условия (4.8)

соответствующая составляющая решения y

не будет затухать.

Таким образом, на всем интервале определения решения

необходимо применять малый шаг интегрирования, что приводит к

огромному числу шагов на больших промежутках интегрирования и

чрезмерному возрастанию машинного времени.

Описанная ситуация встречается при интегрировании жестких

систем уравнений. Жесткие системы характеризуются тем, что среди

собственных чисел матрицы Якоби



f/y имеются большие по

абсолютной величине, которые обязательно обладают большой по

модулю отрицательной действительной частью, а собственные числа

с положительной вещественной частью имеют малую величину. Для

жестких применяются специально сконструированные численные

методы.

Если же



> 0, то неравенство 1+ h



> 1выполняется для

любого h > 0,т.е. на каждом шаге процесса интегрирования

наблюдается увеличение локальных ошибок. Однако на этот раз

неустойчивость не столь очевидна, поскольку доминирующей

является первая компонента глобальной ошибки E

n+1

: каким бы

малым ни был выбран шаг h, разность (exp[(n + 1)h



] —(1+ h



)

n+1

) при

достаточно большом n всегда будет превосходить (1+ h



)

n+1

. Такого

рода задачи называют плохо обусловленными.5

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

ошибок на общий случай y’ = f(x, y), то тогда имеют место

следующие утверждения:

а)если f/y < 0, то влияние локальных ошибок уменьшается при

значениях h, удовлетворяющих условиям устойчивости;

б)если f/y > 0, то влияние локальных ошибок увеличивается

независимо от того, насколько малым выбран шагh.

Вомногих задачах знак f/y меняется на интервале

интегрирования,т.е. в процессе интегрирования локальные ошибки

то увеличиваются, то уменьшаются. Имеет смысл выполнять

интегрирование так, чтобы на каждом шаге знать знак f/y и тем

самым по крайней мере контролировать ситуацию. Впростейшем

случае y’ =



> 0, лучше интегрировать в направлении

уменьшения x,т.е. с отрицательным шагом.

21.ГантмахерФ.Р. Теория матриц. —М.: Гл. ред. физ.-

мат.лит.

изд-ва Наука, 1966. —576с.

Применение функций от матриц к интегрированию системы

линейных

дифференциальных уравнений с постоянными

коэффициентами

Рассмотрим систему однородных линейных дифференциальных

уравнений с постоянными коэффициентами первого порядка:













nnnnnn

xaxaxadtdx

.../

2211

32321313

22221212

12121111



, (38)

где t —независимое переменное;

, x

,..., x

—неизвестные функции переменнойt;

(i, k = 1, 2,..., n) —комплексные числа.

Система уравнений (38) может быть записана в виде одного

матричного дифференциального уравнения

Axdtdx /

. (39)

Будем искать решение системы дифференциальных уравнений,

удовлетворяющее заданным начальным условиям:



. (40)

Разложим искомый столбец x в ряд Маклорена по степенямt:

...

0000



xtxxx



, где









и т.д. (41)

Ноиз (39) почленным дифференцированием находим:



















, … (42)

Подставляя в (39) и (42) значение t = 0, получим:

xAx 



xAx 



,…

Теперь ряд (41) запишетсятак:

...

!3!2

xexA

xAtxx



. (43)

Таким образом, интегрирование данной системы

дифференциальных уравнений сведено к вычислению элементов

матрицы e

Если в качестве начального значения аргумента взять значение

t=t

, то формула (43) заменится формулой:

  

exp ttAxx 

. (46)

7.ХайрерЭ., НерсеттС., ВаннерГ. Решение обыкновенных

дифференциальных уравнений. Нежесткие задачи: Пер. с англ.

—М.: Мир, 1990. —512с.

III.6. Методы Нордсика

Нордсик предложил класс многошаговых методов, позволяющих

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

методы эквивалентны в некотором смысле неявным методам Адамса.

Идея Нордсика состоит в том, чтобы представить многочлен,

аппроксимирующий y

(x), через производные от нулевого доk-го

порядка включительно, т.е. с помощью вектора (“вектора Нордсика”):

 

nnnn

hyyz















,...,''

,',

. (6.1)

Величины y

(j)

имеют смысл приближенных значений для y

(j)

) где

y(x) — точное решение дифференциального уравнения

 

yxfy ,'

. (6.2)

Чтобы определить процедуру интегрирования, необходимо задать

правило нахождения Z

n+1

по известным Z

и дифференциальному

уравнению (6.2). При использовании разложения в ряд Тейлора

(например, при k = 3) такое правило имеет вид:

yhyy

nnnnn





'''

432

, (6.3)

yhyh

nnnn





4'''

3''

2''

432

nnn





6'''

3''

432





4'''

'''

где значение e выбирается таким образом, чтобы выполнялось

равенство

 

111





nnn

yxfy

. (6.4)

Подстановка (6.4) во второе соотношение из (6.3) дает

 

nnn

fyxfhe



 11

где

'''

3''

nnn

yhfh 

Сучетом этого выражения для e метод принимает вид:

 

nnnnnnnn

fyxf

yhyy 

 11

'''

 

nnnnnnn

fyxfhy

yhyh 

 11

,'''

3''

2''

, (6.6)

 

nnnnnn

fyxf







 11

'''

3''

 

nnnnn

fyxfhy



 11

,'''

'''

Первое уравнение— неявное относительно y

n+1

, а остальные—

явные. Следует отметить, что если величины y

(j)

достаточно точно

аппроксимируют y

(j)

), то значение e аппроксимирует y

(4)

). Сточки

зрения точности это свойство выглядит привлекательно. Но, к

сожалению, метод (6.6) является неустойчивым.

Чтобы преодолеть указанную трудность, Нордсик предложил

заменить постоянные 1/4, 1, 3/2, 1, стоящие в (6.6) перед скобками,

на произвольные значение (l

, l

) и использовать эту

дополнительную степень свободы для достижения устойчивости.

Втаблицах 6.1и 6.2даны векторы l, соответствующие неявным

методам Адамса и ФДН-методам.

Таблица6.1

Коэффициенты l

k-шаговых неявных методов Адамса

6

k = 1 ½ 1

k = 2 5/12 1 1/2

k = 3 3/8 1 ¾ 1/6

k = 4 251/720 1 11/12 1/3 1/24

k = 5 95/288 1 25/24 35/72 5/48 1/120

k = 6 19087/60480 1 137/120 5/8 17/96 1/40 1/720

Таблица6.2

Коэффициенты l

k-шаговых ФДН-методов

 l

6

k = 1 1  1

k = 2 2/3 1 1/3

k = 3 6/11 1 6/11 1/11

k = 4 12/25 1 7/10 1/5 1/50

k = 5 60/137 1 225/274 85/274 15/274 1/274

k = 6 20/49 1 58/63 5/12 25/252 1/84 1/1764

13.СамарскийА.А. Введение в численные методы. —М.:

Наука.

Гл. ред. физ.-мат. лит., 1982. —272с.

V. ЗАДАЧА КОШИ ДЛЯ ОДУ.

§1. Методы Рунге-Кутта.

4.Схемы Рунге-Кутта.

Весьма распространены на практике схемы Рунге-Кутта второго и

четвертого порядков точности.

Вычисления по схеме Рунге-Кутта второго порядка проводятся в

два этапа. Напервом этапе находится промежуточное значение y

по схеме Эйлера с шагом



 

nnn

yxfhyy ,



на втором этапе находится значение y

n+1

по формуле:

   

 

nnnnnn

yhxfhбyxfбhyy ,,1







где б > 0,



> 0—параметры. Исключая y

получим для y

n+1

схему

 

      

nnnnnn

yxfhyhxfбyxfб

,,,1









. (13)

Порядок точности схемы зависит от параметров б и



Выражение для невязки схемы (13) имеетвид

   

 

0''2/1 hxyбh





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

при выполнении условия

2/1



. (14)

21.Математическое моделирование высокотемпературных

процессов в энергосиловых

установках /В.Е.Алемасов,А.Ф.Дрегалин,В.Г.Крюков,

В.И.Наумов. —М.: Наука, 1989. —256с.

6.КИНЕТИЧЕСКАЯ ОДНОРЕАКТОРНАЯ МОДЕЛЬ.

6.2.5Алгоритм расчета параметров К-модели.5

Для интегрирования смешанной дифференциально-

алгебраической системы уравнений особенно удобна схема Ньютона.

Предварительно дифференциальные уравнения необходимо

представить в виде однородных конечно-разностных уравненийвида:

     

 

0,1,

111



 nnnnnn

yxfsyxfshyyF

, (6.31)

где h —шаг интегрирования, s —параметр аппроксимации (как

правило, s = 0.4).

Кпреобразованной в соответствии с (6.31) системе уравнений

применяется итерационная схема Ньютона

     



1 





nnn

yyy

, (6.32)

(где



—номер итерации,



n+1

(



)

—поправки на



-й итерации),

включающая решение системы линейных уравненийвида

 

 

   



nnn

FyyF /

, (6.33)

где [F/y

] —матрица частных производных (якобиан). Расчет

проводится до получения таких значений у, которые отличаются в

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

величину. Затем производится следующий шаг интегрирования. При

этом может быть использовано несколько модификаций схемы

Ньютона:

а)классическая— с пересчетом якобиана на каждом

приближении каждого шага интегрирования;

б)модифицированная— якобиан перевычисляется только на

первом приближении каждого шага интегрирования;

в)модифицированная с “замороженным” на нескольких шагах

якобианом.

Наиболее сложной частью схемы Ньютона является вычисление

частных производных якобиана, и от того, каким способом они

определяются, в решающей мере зависит эффективность этой схемы.

Применение аналитического вычисления частных производных

существенно эффективней, чем их численный расчет. Вэтом случае

при использовании модифицированных схем Ньютона

“замороженный”якобиан можно сохранять на большем числе шагов

интегрирования.

3.БабушкаИ., ВитасекЭ., ПрагерМ. Численные процессы

решения дифференциальных уравнений. Пер. с англ. —М.:

Мир, 1969. —368с.

3.ЗАДАЧИ КОШИ ДЛЯ ОДУ.

Проблемы построения численного алгоритма можно разделить на

три группы:

а) Локальные свойства алгоритма.5

Задача заключается в выборе такого алгоритма для определения

, чтобы разность y(x

) —y

была минимальной. Здесь y

обозначает

приближенное значение точного решения y(x); предполагается, что

y(x

n-j

), j = 1, 2,..., m, известны. Асимптотическая оценка точности

алгоритма для x

n-j

- x

n-j-1

= h 0определяется двумя числами



C,если

 



hhCyxy

0

Решение предполагается достаточно гладким, а

параметр



характеризует алгебраическую степень точности.

б) Глобальные свойства алгоритма.5

Численный метод включает в себя циклическое повторение

алгоритма. Его характерная особенность состоит в том, что

при x

-5x

n-1

0число повторений алгоритма, требуемое для

определения значения решения в точке x, возрастает.

Вследствие этого может произойти накопление ошибок,

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

расчет может стать совершенно бессмысленным. Основной

задачей становится выбор такого алгоритма, который бы не

приводил к накоплению ошибок,т.е. чтобы y

y(x

) при

x

-5x

n-1

0и xx

. в) Устойчивость алгоритма как

численный процесс.5

Некотрые сходящиеся алгоритмы могут сопровождаться

численной неустойчивостью.

3.1.2.5Оценки ошибок.

Говоря об оценке погрешностей, следует различать две

проблемы: априорные оценки, получаемые до расчета, например при

выборе метода, и апостериорные оценки, связанные с конкретным

счетом на самой машине.

1.ПолакЛ.С., ГольденбергМ.Я., ЛевицкийА.А. Вычислительные

методы в химической кинетике. —М.: Наука, 1984. —280с.