Слаботочка Книги

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 [119] 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208

§ 1. Решение задачи Коши с помош;ью формулы Тейлора 361

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

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

§ 1. Решение задачи Коши с помош,ью формулы Тейлора

Один из простейших по своему описанию методов решения задачи Коши основан на использовании формулы Тейлора.

Пусть требуется найти на отрезке [хр, хр + X] решение дифференциального уравнения

у---fix, у) (1)

при начальном условии у(жо) = Уо; /(ж, у) -функция, аналитическая в точке {хо,уо). Дифференцируя (1) по ж, имеем соотношения

у = /.(ж, у) + fyix, у)у,

у = /хх(ж, у) + 2/ Дж, у)у + fyy{x, у)у + fyix, у)у ,....

Подставляя х = хр и у ~ Уо в (1) ив последние соотношения, последовательно получаем значения

yixo), /Ы, у {хо),...;

Таким образом, можно написать приближенное равенство

Если значение ж - жо больше радиуса сходимости ряда

-{х-хо) ,

ТО погрешность (2) не стремится к нулю при п - оо и предлагаемый метод неприменим.

Иногда целесообразно поступить следующим образом. Разобьем отрезок [xq, xq + X] на отрезки [xj-i, Xj], j = 1,... ,N. Будем последовательно Получать приближения yi к значениям решения y{xj), j = 1,... ,N, по следующему правилу. Пусть значение уг уже найдено, вычисляем значения в точке Xj производных у* решения исходного дифференциального

уравнения, проходящего через точку {xj, yj). На отрезке [xj, Xji] полагаем , ч

у{х) Zj{x.) = -i-(rc - XjY (3)

И соответственно берем

Уу+1= Zj(Xj+i). (4)

Рассмотрим случай, когда Xji - Xj = h. Если бы значение yj совпадало со значением точного значения y{xj), то погрешность от замены Уj-.l на Zi(Xjji) имела бы порядок 0(/i ). Поскольку мы вносим погрешность на 0(/i ) отрезках, то можно ожидать, что при уменьшении шага сетки будет вьшолняться соотношение

max\yj~y{x,)\ = 0{hn.

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

При использовании этого метода нужно вычислять значения функции / и всех ее производных fxjym-J при т < т.е. вычислять п(п ~Ь1)/2 значений различных функций. Это требует написания большого числа блоков вычисления производных, что противоречит основной тенденции упрощения отношений между пользователем и ЭВМ.

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

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

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

y{x + h)=y{x)+ / y{x+t)dt. (2)

При замене интеграла в правой части на величину hy{x) погрешность имеет порядок 0{h), т.е.

ij{x + h)= у{х) + hyix) + 0{h).

Поскольку у{х) = /(ж, у{х)), отсюда имеем

у{х + h)= у{х) -Ь /г/(ж, у{х)) Л- Oih ).

Отбрасывая член порядка 0{h) и обозначая х = Xj, х + h = Xj+i, получим расчетную формулу Эйлера (1). Для получения более точной расчетной формулы нужно точнее аппроксимировать интеграл в правой части (2). Воспользовавшись квадратурной формулой трапеции, получим

у{х + h)= у{х) + - [у{х) + у{х + /г)) + Oih),

иначе.

у{х + h)= у{х) + I (/(ж, у{х)) + fix + h, yix + h))) + Oih% (3) соответствующая расчетная формула

Уз+i = Уз + \{fij Уз) + fixj+ъ yj+i)) + 0{1?) (4)

называется неявной формулой Ада.иса второго порядка точности. В некоторых случаях, в частности, когда / линейна по у, это уравнение может быть разрешено относительно yj+i. Обычно же это уравнение неразре-пшмо явно относительно yj+i, поэтому произведем дальнейшее преобразование алгоритма.

Заменим у{х + h) в правой части (3) на некоторую величину

у* = у{х + h) + Oih). (5)

Тогда правая часть изменится на величину

I (/(ж + /I, у*) - fix + h, yix + h))) = fyix + h, y) [y* - yix + h)),

§ 2- Методы Рунге-Кутта

В частном случае п = 1 формула (1.3) имеет вид

yj+i = yj + hf{xj,yj), h = Xj+i~Xj. (1)

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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 [119] 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208