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

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

МолчановИ.Н. Машинные методы решения прикладных

задач.

Дифференциальные уравнения. —Киев: Наукова думка, 1988.

— 344с.

Монография посвящена численным методам решения задач

Коши для систем ОДУ, краевых задач для уравнений и систем

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

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

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

Устойчивость решений ОДУ.

Рассмотрим задачу с начальными условиями для системыОДУ

 

0 uu



. (1.29)

Будем полагать, что все собственные значения матрицы A различны.

Тогда существует такая матрица C, что C

-1

AC =



,где

















.00

....

0.0

, (1.30)

(



— i-е собственное значение матрицыA).

Если ввести замену переменных

   

xzCxu 

то из системы (1.29) получаем систему уравнений

 





, где

ni ,...2,1

(1.31)

Таким образом, изучить устойчивость решения системы (1.29) можно

с помощью системы (1.31). Более того, устойчивость решения

системы дифференциальных уравнений (1.31) можно рассматривать

на примере уравнения





, (1.32)

которое называют тестовым.

Предположим, что





+ j



—постоянная из поля комплексных

чисел ( j

= -1). Решение уравнения асимптотически устойчиво, если

Re  < 0, устойчиво, если Re  = 0, и неустойчиво, если Re  > 0. ( Re 

—реальная часть комплексного числа).

Рассмотрим теперь задачу Коши для системыОДУ

 

uxF

,

, (1.33)

где F(x, u) —вектор-функция, F(x, u) = (f

(x,u),.., f

(x,u))

. Качественное

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

решения u

(x) можно провести следующим образом. Разложим F(x,u) в

окрестности u

(x) в ряд Тейлора, в результате получим:

       

...,,,

111

 uuuxJuxFuxF

Тогда:

   

, uuuxJ



, (1.34)

где J(x, u

) —матрица Якоби. Если изменение J(x, u

) на некотором

интервале x достаточно мало, то J(x, u

) можно заменить локально

постоянной матрицейA. После приведения к главным осям матрица A

принимает вид (1.30). Таким образом, и в этом случае исследование

устойчивости решения системы (1.33) может быть сведено к

исследованию устойчивости решения уравнения (1.32).

Понятие о жестких дифференциальных уравнениях.

РассмотримОДУ

  

 

0 uu

xpuA





, (1.35)

где A = const. Точное решение этого уравнения:

      

xpepuxu

 0

. (1.36)

Если A —большая положительная величина, то решение задачи (1.35)

неустойчиво. Если A —малая величина, то решение задачи (1.35)

нейтрально устойчиво. Если A —большая по модулю и отрицательная

величина, то каким бы ни было выбрано u

через достаточно малый

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

слоем), кривая решения становится тождественной кривой частного

решения p(x) уравнения (1.35). Эта сверхустойчивость решения

дифференциальной задачи является идеальной в смысле

распространения наследственной ошибки в дифференциальном

уравнении, но она создает ряд трудностей численного

интегрирования на ЭВМ. Одна из них состоит в том, что хотя решение

(1.36) за пределами переходного участка ведет себя как и функция

p(x) и практически не зависит от A (при A < 0), тем не менее из

условия устойчивости шаг численного интегрирования приходится

выбирать зависящим отA. Чем большеA, тем меньше шаг

интегрирования (жестче ограничения на него). Такие задачи

называют жесткими.

Определение. Задача Коши (I.33) называется жесткой на

некотором интервале [x

, X], если для x из этого интервала

выполняются условия:

1)

0Re 



ni ,...,2,1

;

2)

 

Remin

Remax









где



—собственные значения матрицы Якоби, вычисленные на

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

Величину s(x) называют локальным “коэффициентом жесткости”

задачи. Если s(x) является величиной 0(10), то задачу можно считать

жесткой. Вряде прикладных задач (химическая кинетика, процессы

управления, теория электрических цепей) коэффициент жесткости

достигает величины 0(10

Это определение “жесткости” системы нам кажется не совсем

точным. Действительно, если max Re(-



) = 10

-2

, а min Re(-



) = 10

-8

так что их отношение равно 10

, то по определению система является

“жесткой”. Вто же время такая система позволяет из условий

устойчивости выбрать шаг интегрирования порядка100.

Определение, возможно, надо дополнить словами: большие по

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

Многошаговые методы (Методы Адамса).

Имеются два важных отличия между явными и неявными

методами Адамса. Первое заключается в том, что поскольку модули

коэффициентов неявных методов меньше, чем явных, это приводит

не только к меньшим ошибкам округления, но и к меньшей

погрешности аппроксимации (погрешности метода) при одном и том

же порядке точности метода.

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

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

количестве точек. Порядок неявного метода на единице выше.

Самым простым способом решения неявного уравнения является

метод предиктор-корректор (комбинация явных и неявных методов).

АрушанянО.Б., ЗалеткинС.Ф. Численное решение

обыкновенных дифференциальных уравнений на Фортране.М.:

Изд-во МГУ, 1990. —336с.

Монография посвящена процессам организации вычисления

на ЭВМ решений ОДУ. Рассматриваются способы решения

задачи Коши для уравнений первого и второго порядков

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

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

интегрирования, эффективные матричные методы решения

линейных систем с постоянными коэффициентами, методы

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

Фортран, реализующих различные численные алгоритмы и

составляющих функциональное наполнение пакета программ

для решения ОДУ, разработанного в НИВЦМГУ.

I. Машинная арифметика и простейшие вычисления.

При выполнении расчетов на ЭВМ постоянно приходится иметь

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

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

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

процессов конечными (погрешность численного метода), например в

результате отбрасывания остаточного члена формулы; в-третьих,

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

машине (ошибки округления). Каждый из этих типов ошибок

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

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

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

окончательного результата…

II. Одношаговые методы решения ОДУ.

...будем применять численные методы, которые позволяют

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

1.Общее представление одношаговых методов.

Пусть требуется найти решение задачи

 

yxfy ,'

, (1)

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

 

yxy 

, (2)

на отрезке [ x

, x

5+ X]. Возьмем разбиение отрезка точками:

Xxxxxx



0210

...

Этот набор точек называется сеткой, а точки x

—узлами сетки.

Одношаговые численные методы последовательно дают

приближения y

к значениям точного решения y(x

) в каждом узле x

сетки на основе известного приближения y

n-1

к решению в

предыдущем узле x

n-1

. Вобщем виде их можно представитьтак:

 

,,,,





nnnnn

yyxxfFy

. (4)

Для явных одношаговых методов функция F не зависит от y

 

,,,





nnn

yxhfFy

, (5)

где h = x

- x

n-1

— шаг сетки.

2.Метод рядов Тейлора. Пусть правая часть f(x,y)

дифференциального уравнения (1) имеет непрерывные частные

производные до порядкаs. Тогда искомое решение y(x) имеет

непрерывные производные до (s+1)-го порядка включительно.

Точное значение решения в узле x

по формуле Тейлора:

 

)(

!1!

...''

)1(

)(

001











yhyxy

, (6)

где

)(

xyy



xxh 

xx 



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

требуется использовать все члены формулы (6). Производные,

входящие в правую часть формулы (6) могут быть фактически

найдены:

 

000

,' yxfy 

 

''''

fffy 

 

'''''''2'''''

yyxyyxyxx

fffffffffy 

...

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

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

является существенным недостатком данного метода.

3.Явные методы типа Рунге - Кутта. Рунге предложил

следующую идею, основанную на вычислении приближенного

решения y

в узле x

+h в виде линейной комбинации с постоянными

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

     

hkphkphkpyy

qqqqq

 ...

221101

, ( 7)

где

   

001

, yxfhhk 

    

hkyhxfhhk

1210202

, 



…

     

 

hkhkyhxfhhk

qqqqqq 11,1100

...,







Числа





и p

выбираются так, чтобы разложение выражения (7) по

степеням h совпадало с разложением (6) до максимально возможной

степени при произвольной правой части f(x, y) и произвольном

шагеh.

Это эквивалентно следующему. Если ввести вспомогательную

функцию:

     







iqiq

hkpyhxyh



, (8)

то ее разложение по степеням h должно начинаться с максимально

возможной степени:

 

)1(











, (9)

Если можно определить эти постоянные так, чтобы разложение



(h)

имело вид (9), то говорят, что формула (7) с выбранными

коэффициентами имеет порядок точности5s.

Величина

   

101

yhxyh





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

погрешностью метода, а первое слагаемое в(9)

 

)1(









, (10)

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

Доказано, что если q = 1, 2, 3, 4, то всегда можно выбрать

коэффициенты





, p

так, чтобы получить метод типа Рунге-Кутта

порядка точностиq. При q = 5невозможно построить метод типа

Рунге-Кутта (7) пятого порядка точности, необходимо брать в

комбинации (7) более пяти членов.

3.1.5Одночленная формула. Приближение в точке x

5= x

5+ h

ищется ввиде:

 

hkpyy

11101



где

   

001

, yxfhhk 

Коэффициент p

подбирается так, чтобы в разложении по

степеням h функции (8), в данном случае равной:

     

hkpyhxyh

111001





максимальное количество членов обратилось в нуль. Этому

требованию удовлетворяет единственное значение:

p

для которого разложение (9) имеетвид:

 

0''

hfff





Таким образом

 

0001

, yxfhyy 

Это известная формула Эйлера. Метод Рунге-Кутта при q = 1есть

метод Эйлера.

3.2.5Двухчленные формулы. Приближение в точке x

5= x

5+ h

ищется ввиде

   

hkphkpyy

22212101



где

   

001

, yxfhhk 

    

hkyhxfhhk

1210202

, 



Ввыписанной формуле четыре неизвестных параметра:





, p

Чтобы их найти, необходимо построить вспомогательную функцию

(8),равную:

       

hkphkpyhxyh

222121002





и подобрать коэффициенты так, чтобы в разложении этой функции

по степеням h максимальное количество членов обратилось в нуль.

Чтобы обратить в нуль первые две производные:

     

0022212

,10' yxfpp 



     

 

22212222

'21'210''

ffpfp 



необходимо, чтобы искомые коэффициенты удовлетворяли системе

уравнений

.021

,021

,01

2221

222

2221









(11)

Тогда разложение (9) имеетвид

   

 

00'''

h 



Из(11) следует, что













. Если принять



за

свободный параметр,то

221









p





p

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

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

1.Пусть



= 1. Тогда

 

1002

001

kyhxfhk

yxfhk









(14)

Погрешность формулы(14)

 

   

0'''2/''''2''

hffffffff

yyxyyxyxx

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



Данный метод называется методом Хойна. (или методом Эйлера с

коррекцией по средней производной).

(???)Если правая часть формулы (1) не зависит от y, то формула

(14)переходит в квадратурную формулу трапеций.(???)

2.Пусть 

= 1/2. Тогда

 



















002

001

201

xfhk

yxfhk

kyy

(16)

Погрешность формулы (16):

 

   

0'''4/''''2''

hffffffff

yyxyyxyxx





Метод называется методом Эйлера с коррекцией по средней

точке

Если правая часть формулы (1) не зависит от y, то формула (16)

переходит в квадратурную формулу средних прямоугольников.

3.Пусть 

= 2/3. Тогда

 





















1002

001

kyhxfhk

yxfhk

(18)

Погрешность формулы(18)

  

 

0'''

hffff

yyx





3.3.5Трехчленные формулы. Приближение в точке x

5= x

5+ h

ищется ввиде:

     

hkhkphkpyy

323213201



где

   

001

, yxfhk 

    

hkyhxfhk

1210202

, 



      

hkhkyhxfhk

2321310303

, 



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

система из шести уравнений:

212





32313





333231

 ppp

 

333322

 pp



 

332

 pp



33322

 p



Тогда разложение (9) имеетвид:

 

h 



Пример. Пусть



= 1/2,



= 1. Тогда:

 

.2,

21003

002

001

321

kkyhxfhk

xfhk

yxfhk

kkk



































. (21)

3.4.5Четырехчленные формулы. Существует семейство формул

типа Рунге-Кутта (7) четвертого порядка точности, для которых

q=s=4.

Пример: (классический метод Рунге-Кутта):

 

3004

003

002

001

4321

kyhxfhk

xfhk

yxfhk

kkkk

















































(22)

Тогда разложение (9) имеетвид:

 

120

h 



3.5.5Формулы порядка выше четвертого. Чтобы построить метод

типа Рунге-Кутта (7) порядка выше четвертого, для которого

погрешность на шаге имеет порядок выше пятого для произвольного

дифференциального уравнения (1) с достаточно гладкой правой

частью f(x, y), необходимо, чтобы в формуле (7) число слагаемых k

(h)

было больше пяти: q>5.

Здесь существенно то, что значения правой части k

(h)

дифференциального уравнения входят в формулу (7) в виде

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

метод пятого порядка точности y

5= y



y, использующий пять

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

(h)

будут входить в



y нелинейным образом.

Общим для всех рассмотренных методов численного

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

очередного шага точка, изображающая полученное приближенное

решение на плоскости (x,y), переходит с одной интегральной кривой

на другую.

4.Сходимость явных одношаговых методов.

Предположим, что требуется найти функцию y(x), которая

является решением дифференциального уравнения (1) и принимает в

точке x

0

некоторое определенное значение

 

yxy 

Может оказаться, что начальное условие y(x

) известно неточно, а

определяется в результате эксперимента, например, с помощью

измерений или в результате решения какой-либо задачи. Вэтом

случае вместо точного начального условия y(x

) приходится

использовать его приближенное значение y

, а вместо задачи Коши

(1) решать задачу

    

xyxfxy

,'

, (25)

 





yxy

, (26)

с измененным начальным условием

000





Ryy

, (27)

Решение задачи (25), (26) зависит от y

и не совпадает с искомым

решением задачи (1), (2). Разность

   

nnn

xyxy





(28)

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

(x).

Разность между значением решения y

) задачи (25), (26) и его

приближенным значением y

, полученным по формуле(5),

 

nnn

yxy 



(29)

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

метода. Погрешность метода на одном шаге интегрирования

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



Вдействительности же вследствие ошибок округления и

приближенного вычисления правой части f(x, y) дифференциального

уравнения вычисления значений y

n+1

по формуле (5) выполняются,

как правило, неточно. Фактически найденные значения y

удовлетворяют не соотношению (5), а условию

 

nnnn

yyxhfF





 11

,,,

(30)

Невязка



называется погрешностью округления на n-м шаге.5

Разность между точным решением задачи y(x

) задачи (1), (2) и

приближенным фактически найденным значением y

n

 





nnn

yxyR

(31)

называется полной погрешностью приближенного решения.

Величина





nnn



(32)

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

Изсоотношений (28), (29), (31) и (32) следует,что:

nnnn





, (33)

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

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

погрешности.

Присутствие в (33) слагаемого



означает, что погрешность

начального значения распространяется на все узлы сетки. Часть

полной погрешности приближенного решения, зависящая от ошибки

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

начальные условия заданы точно, то этого членанет.

Второй член получается за счет того, что мы находим не точное

решение задачи, а приближение к нему по формуле Рунге-Кутта. Это

погрешность метода, и она имеет порядок h

Третий член получается за счет ошибок округления. Впрактике

вычислений на ЭВМ величина



ограничена снизу, если при этом

длина шага h слишком мала, а интервал интегрирования велик, то

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

Вдействительности же ошибки округления



могут иметь

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

знаки могут иметь и погрешности метода



. Отдельные

составляющие, входящие в полную погрешность R

(33), могут давать

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

Изприведенных оценок можно сделать следующий вывод,

подтверждаемый практикой численного решения на ЭВМ

дифференциальных уравнений. Если заранее обеспечена

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

погрешности преобладает либо погрешность метода, либо

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

сделана сколь угодно малой за счет уменьшения шагаh.

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

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

значащих цифр. Однако практические возможности в этом далеко не

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

Извышесказанного следует, что вычислительный процесс

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

баланс между всеми видами погрешности, составляющими полную

погрешность приближенного решения. Этот баланс может быть

достигнут, если надлежащим образом будут согласованы между

собой требуемая точность решения задачи, точность задания

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

интегрирования и используемая длина разрядной сеткиЭВМ.

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

решения.

5.1.5Апостериорная оценка глобальной погрешности метода.5

Правило Рунге апостериорной оценки погрешности состоит в том,

что решение задачи в некоторой точке x

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

вычисляется дважды по одной и той же формуле (7) с разными

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

получения апостериорной оценки погрешности. Обычно в качестве

шагов выбирают h иh/2.

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

можно получить в следующем виде:

 

nnn

yxyR





, (53)

 





yxyR

, (54)

где

— значения решения, вычисленные с шагами h/2иh

соответственно;

— значения погрешностей для решений с шагами h/2иh.

5.2.5Апостериорные оценки локальной погрешности метода.5

5.2.1.Оценка погрешности по правилу Рунге.

(аналогично оценке глобальной погрешности)

5.2.2.Оценка погрешности на основе комбинации формул

различных порядков точности.

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

значений решения в одной точке. Однако эти приближения в отличие

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

разных порядков точности p и s с одним и тем же шагом.

5.2.2.1.Комбинация двух независимых формул. Данный способ

основан на комбинации двух формул вида (7) разных порядков

точности p иs. Пусть p > s, тогда локальная погрешность будет

равна:

 





psps

hyy



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

sps





. (71)

5.2.2.2.Комбинация специально подобранных формул. Если

коэффициенты в формулах таковы,что:





ijij





,...,2,1

, (72)

то



,...,2,1

и для локальной погрешности получается выражение следующего

вида:







sps

kqyy



, (73)

где

iii

ppq



ri ,...,2,1

Такой подход к оценке локальной погрешности позволяет

уменьшить по сравнению с правилом Рунге и оценкой (71) количество

вычислений правой части уравнения(1).

Оценка (73), как и более общая оценка (71), является

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

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