Введение в теорию колебаний и волн. - Рабинович М.И.
Скачать (прямая ссылка):
Ех = max I x\+w - wx\+x - (1 - w) xi I. (6.4.14)
i 1 1
Аналогично определяется Ey. Такую же оценку можно получить и для метода
"предиктор-корректор".
c) В случае использования какой-то схемы линеаризации
(6.4.8) мы можем положить
Ех -- maxjf (дах{+1 + (1 - w)xi, wyi+l + (1 - w) г/() -
- ai - b\xi+l - с[у[+х |. (6.4.15)
Аналогичным образом записывается и оценка Еу. Отметим, что в данном
случае мы находим отклонение в значениях функций / и g, но не в
переменных состояния х и у.
О Обычно норму || || определяют так: ||u|| = max|u;|.-Прим.
ред.
6.4. Методы динамического моделирования
333
Каждую из величин Е = тах(?х, Еу)> приведенных выше, можно использовать
для регулировки шага по времени т. С этой щелью обычно выбирается
некоторое характерное значение е, причем шаг т уменьшается, если Е > е.
Если же Е существенно меньше, чем е, то шаг т увеличивается.
Мы ограничились здесь лишь идейной стороной алгоритма изменения шага т.
Очень вероятно, что действительная величина погрешности решения будет
заметно отличаться от е. Поэтому если мы хотим удостовериться в том, что
заданная точность •обеспечивается, следует провести вычисления при двух
различных значениях е (отличающихся, например, на порядок) и сравнить
полученные результаты. Таким образом, однако, можно •оценить лишь влияние
погрешностей аппроксимации "в направлении оси Ь>. Влияние погрешности
аппроксимации "в направлении координаты 2" и адаптационные алгоритмы для
выбора шага h будут обсуждаться в следующем пункте.
"6,4.3. Автоматическое изменение шага h
В этом пункте мы исследуем проблему автоматического изменения шага h на
равномерной сетке. Для решения вопроса о том, достаточно ли велико п (т.
е. достаточно ли мал шаг К), при заданной "характерной величине
погрешности" е, будем искать оценку погрешности аппроксимации в
направлении "оси г. В принципе мы могли бы провести вычисления для
момента tj дважды - с шагом h и 2h - и результаты сравнить. Такой
численный подход не был бы, однако, достаточно простым и эффективным.
Другой, более подходящий путь заключается в том, чтобы использовать
выражения для остаточного члена в соответствующих разностных формулах.
При замене второй производной в (6.4.3) имеем
Ц(г". ',) = TsrW-.-2*f + *!+,]-';)• <6А16>
где ?e(z*_i, zi+\). Для аппроксимации четвертой производной можно
воспользоваться разностной формулой вида
1F С1' */) ~ 1F M-а - М-i + Щ ~ ixUi + (6-4Л7>
при I = 2,3, ..., п - 2. Тогда оценку погрешности "в направлении оси г"
можно записать в виде
Ех==г=2ШаХп-2 wi ~ 4*'-> + 6xli ~~ 4*' + 1 + *J+2l- (6-4Л8)
334
Глава 6
Аналогичное выражение получается и для Еу (при i = 1 и i = = л- 1 для
четвертой производной можно использовать асимметричные формулы).
Обозначим Е = тах(?х, Еу) ¦
Тогда для изменения (регулировки) п можно использовать следующую
стратегию:
1. Е > е; шаг h слишком велик, и его следует уменьшить вдвое. Значения
решения в узловых точках, которые появятся между исходными узловыми
точками, после уменьшения шага, можно определить с помощью интерполяции,
например, положить
После проведения интерполяции перенумеруем заново узловые точки, удвоим
число п и вдвое уменьшим шаг h; затем продолжим вычисления.
2. Е < е/ш, где значение ю, как правило, выбирается между 4 и 8. Шаг h
излишне мал и может быть удвоен. При этом узловые точки с нечетными
индексами 1, 3, 5 п- 1 выбрасываются, оставшиеся точки
перенумеровываются, шаг h удваивается, а п уменьшается вдвое.
3. В остальных случаях шаг h остается неизменным и проводятся дальнейшие
вычисления.
При формулировке алгоритма удобно задавать максимальный и минимальный
допустимый шаг h.
6.4.4. Адаптивная неравномерная сетка
Пусть у нас имеется некоторая, вообще говоря, неравномерная сетка узловых
точек Zo = 0, Z\, ..., z" - l. Тогда вторую производную в направлении z в
точке (¦?,-, tj) можно заменить трехточечной разностной формулой
При этом разностные уравнения (6.4.3) изменяются только за счет
подстановки формулы (6.4.20) вместо простейшей фор-
*1 + 1/2
*{/2=т13*о + 6*{-4]>
*п-1/2 = If [3*" (r)*п-1 *л-г]-
(6.4.19)
6.4. Методы динамического моделирования
335
д^х
мулы1) для -fcp- Все остальное остается неизменным, включая
вычисления нового слоя ;+1. Погрешность аппроксимации в направлении оси z
оценивается на каждом подынтервале Zi) отдельно. Айгенбергер и Батт
[6.29] предложили следующий подход. Приближенное значение решения при t =
tj и z = (z,--i + zt)/2 отыскивается двумя способами (ниже индекс j
опущен). В первом случае мы получаем его интерполированием по узловым
точкам z,-_2, z,_i, z,-,
1-1/2 HZi - - Zi-2) l~2 + 4Zi-X ~ *i-2) г-' +
Zi+Zi-X ~2zi-2 4(Zi ~ Zi-2)
xit (6.4.21)
во втором случае - интерполированием по узловым точкам
-1" li
"(2) _ 2zi+t-zi~zi-1 2z.+l-zi-z._l
'-1/2 4zi+x-zi-1) г-1+ 4zi+i-zi) 1
(6Л22>
Для оценки погрешности аппроксимации на интервале (z,-i, Zi) используется
величина
Е1Х = | хг-1/2 - xf-i/21; (6-4.23)
?'у определяется аналогично. Тогда Е? = max (Е1Х, Е1у). Далее сетка
