Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Как известно, при любом выборе узлов (j = 1, 2, ..., s) квадратурная формула точна в классе полиномов степени s — 1. Имея в своем распоряжении s свободных параметров %К можно выбрать их так, что квадратурная формула станет точной в классе полиномов степени (s — l) + s = 2s— 1. Такие узлы называются гауссовыми (для них существуют таблицы). Если использовать схему интерполяционной) типа именно с этими узлами, 5-устойчивость такой схемы доказывается просто и не требуется проверять свойства сложно определяемой матрицы Q. Другими словами, 5-устойчивость оказывается закономерным свойством схемы интерполяционного типа с гауссовыми узлами.
Предварительно докажем одно полезное свойство схем интерполяционного типа, справедливое при любых, в сущности, узлах \К Напомним, что в такой схеме фигурируют значения Yj и
240
ЧИСЛЕННЫЕ МЕТОДЫ МАТЕМАТИЧЕСКОЙ ФИЗИКИ
[Ч. II
fj = f(yi) (у = I, 2, ..., s), определяемые решением системы уравнений (38), и интерполяционный полином ??(t), вычисляемый по узлам V и значениям /А Расширим эту информацию, добавив узел t0 = tn, и зададим значения Y0 = хп, Yj (/' = 1, 2, s). По сеточ-
ной функции {/у, Yj} (у = 0, 1, ..., s) построим интерполяционный полином p(t) степени s: p(tJ) = YJ (j = 0, I, ..., s). Имеет место следующее утверждение.
Утверждение. Полином p(t) в точках Iі (У = 1 ,2 ,..., s) удовлетворяет исходному дифференциальному уравнению:
Ktj)= f[p{t}) ], /=1,2 ,...,s. (59)
Доказательство. Имеем очевидное тождество
P(Ii) - Р(*°) = J P(t) dt. t0
Ho p(t‘) = Yi, P(t°) = xn, p(t) — полином степени s — I; интеграл точно вычисляется по квадратурной формуле с коэффициентами Xaij. Итак,
5
Yi = Xn +xj Uij p(tj), і= 1,2, ...,s.
J = і
Видно, что величины p(tJ) и /7 = f(Y*) — f[p(V)\ удовлетворяют одной и той же системе уравнений (38) (линейной относительно этих величин). Предполагая невырожденность матрицы А, получаем совпадение:
p(tJ) = fj = f[p(tJ)], У = I, 2, ..., s.
Кроме того, очевидно, что p(t) так как р и -2* — полиномы
степени s — 1, совпадающие в s точках Iі (У = 1, 2, ..., s).
Таким образом, схему интерполяционного типа можно получить так называемым методом коллокации. Решение на интервале [tn, tn + т] ищется (приближенно) в виде полинома степени s, для определения которого используются условия коллокации (т.е. выполнения уравнения в нескольких точках tj, число которых, естественно, связывается с числом свободных параметров полинома):
P(tj) = f[p(tj)\, / = 1,2,...,5.
К этому добавляется, конечно, условие p(t°) = хп.
ЖЕСТКИЕ СИСТЕМЫ ОДУ
241
Пусть теперь Xli й хп — два разных решения, полученных по
схеме (36), (38) с гауссовыми узлами Iі. Используем полиномы p(t) и і?(0і определенные выше для решения хп (аналогичные полиномы для хп обозначим p(t) и &(t)). Как оказалось, для этих
полиномов точки V являются точками коллокации, т.е. выполняются соотношения (/= 1, 2,____,
P(tn) = Xn, p(tn + Vx) = Yj, 'p(tn + Ут) = /(У').
Такие же соотношения, естественно, выполняются и для второго решения.
Рассмотрим величину
rO) = IIpCO -p(0II2 = р0) - КО)-
Нас интересует сравнение величин r(tn) и r(tn + г). 5-устойчивость будет установлена, если будет показано (ограничимся классом диссипативных систем), что r(tn + т) ^ r(tn). Очевидно,
/„+г
rOn + о -rOn) = Siidt
гп
и надо оценить интеграл от г:
'rO) = 2(р(О-рО), рО) — рО))-
Заметим, что р и р — полиномы степени s; следовательно, r(t) — полином степени 2s, a r(t) — полином степени 2s — 1. Для таких полиномов гауссова квадратура точна и интеграл ^ г dt можно вычислить по квадратурной формуле:
t+x s s
! Г dt = ^bJ r(tn + %h) =^bJ Cp-Р, P-р) \ і +^r
J=I J=1
В силу (59) имеем (у = 1, 2, ..., s)
rOn + V*) = 2(Р-Р, Р~Р) І ,+& = (fiYJ)-fiYJ), Yi-Yj) 0.
Здесь, конечно, использовалась диссипативность системы уравнений к = f(x). Так как коэффициенты гауссовых квадратур положительны, получаем требуемый результат — аттрактивность разностной схемы.
242
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. И
§ 18. Жесткие линейные краевые задачи
Рассмотрим численное решение специального класса линейных краевых задач, неявно содержащих большой параметр. ,Формально задача ставится стандартно: на интервале [О, Г] решается линейная задача
^ = Ax + a(t), xERn, (1)
с краевыми условиями (при t = 0 и t = T соответственно)
(Ii, х(0)) = bt, і= I, 2, ...,k< п,
(li,x(T)) = bi, і = к + 1, к+ 2, ..., п, (2)
IiERn, IIZiII = I, г= 1,2, ...,п.
Матрицу А в теоретических рассуждениях будем считать постоянной, хотя все последующее относится и к случаю, когда А зависит от t, но изменяется со временем в некотором смысле (ниже будет уточнено, в каком именно) «медленно».
Жесткость системы (1) означает, что собственные значения матрицы А можно разделить на три части:
1) левый жесткий спектр (собственные значения AJ) характеризуется соотношениями
ReA ~^—L, I Im Ar \ < L, і = 1, 2, ..., Г\
2) правый жесткий спектр (собственные значения At) характеризуется неравенствами
