книги / Метод продолжения решения по параметру и наилучшая параметризация в прикладной математике и механике

В конце ее обхода ошибки, подсчитанные по формулам (4.66), были равны Д[ = 0,036, Дг = -0,042.

Начальное значение вектора Z принималось в виде У^=(1,1,1)/л/3. Задача также решалась при помощи программы DA1IN. Время счета

составило 255 с., а ошибки Д | = -0,012, Дг = -0,037. Система уравнений (4.54) имела вид

’ (2Р - Ь2)уУ' - \ЬЧ - P(2t ~ а)]Т' =

= - [IP1у + (2Р - b2)Y]Y + \Ь2Т - P'(2t - а) - 2РТ]Т, 4yY' + хХ' + tT' = -4У 2 - X 2 - Т2,

, YY' + Х Х ' + ТТ 1 = 0,

где Р>= - 2 хХ - 6yY - аТ.

Система дифференциальных уравнений (4.18), (4.55) интегрирова­

Q = у/4(а2 + Ь2) + Зо2Ь2 - Ь4,

в которых значения функций Y , X, Т найдены из решения системы уравнений (4.67) при условиях (4.61).

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

Следует обратить внимание на использование при интегрировании задач, нелинейно зависящих от производных у,-, программы DA1EXP, которая при несколько большем по сравнению с программой DA1ILN

времени счета приводит к результатам практически той же точности. Но в данном случае линеаризация уравнений не требуется. Поэтому, если перевод задачи в расширенное пространство решений, о котором говорилось выше, не приводит к недоразумениям, по-видимому, разумно использовать программу DA1EXP.

Обратим внимание еще на один момент. В [77] показано, что ошибка из-за плохой обусловленности концентрируется в алгебраической пере­ менной и ее нет в дифференцируемой переменной. Это предлагается использовать для контроля точности вычислений. Анализ численных ре­ зультатов, полученных в данном исследовании, показывает, что отмечен­ ные ошибки близки. Это может указывать на хорошую обусловленность вычислительного процесса.

В качестве более сложной рассмотрим модельную задачу о раскрытии купола осесимметричного парашюта под действием заданного перепада давлений р. Аэродинамическое влияние строп не учитывается. Материал стропы и ленты радиального каркаса подчиняется при растяжении закону Гука

{О,

где N — усилие; Е — модуль упругости, имеющий размерность силы; dS и da — длина элемента дуги соответственно в деформированном и недеформированном состояниях.

Если обозначить через т о и тп массы единицы длины стропы или ленты радиального каркаса с прилегающей к ней тканью вме­ сте с кольцевым каркасом соответственно в недеформированном и де­ формированном состоянии, то с учетом диссипативных сил уравнения движения осесимметричного парашюта примут вид [54, 7]

(4.68)

М2 d s \ l d s

Здесь х, г — оси неподвижной системы координат (рис. 4.3) с началом в куполе парашюта; t — время; е — коэффициент диссипации; М — число строп парашюта; параметр S = 1 для ленты радиального каркаса и S = 0 для стропы.

где те — масса единицы длины стропы; Le — длина стропы; До — раскройный радиус купола парашюта.

С учетом введенных обозначений система уравнений (4.68) запишется в виде

Вэтих уравнениях звездочка при безразмерных величинах опушена;

v— Ее/Е, если I > 1, v = 0, если / < 1, Ес — модуль упругости материала стропы.

В силу симметрии рассматриваем только половину парашюта, по­ этому граничные условия будут иметь вид

вг(т, Lc+ 1)

1(т, £„ + !) =

ds

Для решения задачи (4.69)-(4.71) применялся метод прямых по ко­ ординате в, согласно которому безразмерная начальная длина ленты парг плотной системы Lc+ 1 делилась на (п - 1) равную часть с шагом h — (Lc+ l)/(n - I). Производные по координате 8 аппроксимирова­ лись центральными разностями, имеющими второй порядок точности, которые, записанные для функции Z(s) в fc-ой точке (1 < к <п), будут иметь вид

Значения производных, входящих в граничные условия, опреде­ лялись из следующих соображений. Запишем в окрестности начала координат (Z = Zi) представление функции Z(s) по формуле Тейлора, ограничившись членами второго порядка относительно h (здесь Z(s) это функция Xi(s) или r*(s))

В центре купола с учетом того, что dx/ds = 0, получаем выражение для хп = (4®„_| - ®„-2) /3.

согласуется с экспериментальными данными [7]. Таким образом, началь­ ная форма парашюта может быть полностью описана, например, через радиус его входного отверстия R» (рис. 4.3), R, = 0,15Ло [7]. Скорость точек парашютной системы в начальный момент полагается равной ну­ лю. Отсчитывая аргумент Л от начальной точки задачи (4.69)—(4.71), начальные условия примем в форме

где величины х ^, од описывают контур, показанный на рис. 4.3 штриховой линией.

Задача (4.74), (4.7S) интегрировалась с помощью программы РС1

при

Задача интегрировалась с начальным шагом Ад = 0,05 и точностью КГ4. На рис. 4.3 приводится конфигурация парашюта для пяти моментов безразмерного времени то = 0, т\ = 2,5, тг = 5,0, тз = 7,5, Г4 = 10,0. На рис. 4.4 для тех же моментов времени показано изменение усилий вдоль стропы и ленты купола парашюта. На рисунках звездочкой обозначено начало купола парашюта. Увеличение числа разбиений вдоль координаты S к существенным изменениям ситуации не привело.

Заметим, что при решении задачи (4.74), (4.75) без применения Л-преобразования выяснилось, что вычислительный процесс зависит от выбора начального шага интегрирования Лщ. Так, при А*о ^ 1 про­ исходил останов из-за переполнения памяти ЭВМ. Это, по-видимому, объясняется тем, что при поиске подходящего шага интегрирования А* получаются большие значения правых частей дифференциальных урав­ нений. При решении же А-преобразованной системы таких ситуаций не возникало даже при начальном шаге интегрирования Лдо = 10.

Изучим еще одну модельную задачу. Кинетические уравнения ней­ тронов в реакторе в одномерном случае можно представить в виде [97,41]

Здесь первое уравнение является кинетическим уравнением запаз­ дывающих нейтронов, а второе — кинетическое уравнение мгновенных нейтронов; u(x,t), v(x, t) — концентрации мгновенных и запаздыва­ ющих нейтронов; a(x,t), b(x,t) — заданные функции, е — малый параметр, равный обратному значению скорости мгновенных нейтронов. Исследуем вырожденный случай, когда параметр е = 0, тогда систе­ му (4.76) можно записать в виде

dv

— = av + bu+g(u, v, t),

Очевидно, что в начальный момент функция u(x, t) должна удов­ летворять равенству

02и(®, 0)

+ /(«(®,0),0,0) = 0.

дх2

Пусть функция /(«(ж, 0), 0,0) = 0, тогда

Эта система является системой дифференциально-алгебраических уравнений. Принимая во внимание соотношения (4.79), (4.80), началь­ ные условия будут иметь вид

Сформулируем задачу (4.81)—(4.83) относительно наилучшего аргу­ мента Л. Дифференцируя алгебраические соотношения и учитывая смысл наилучшего аргумента, получаем

Для приведения этой неявной системы дифференциальных урав­ нений к нормальной форме, разрешаем ее относительно производных, представляя последнее уравнение в виде

где индексом (fc - 1 ) помечена функция, найденная на предыдущем шаге процесса интегрирования, тем самым используя аргумент, близкий к наилучшему.

Так как начальная точка не является предельной по t, то начальное значение вектора Z = (г/1)А, и2,Л>. . . , v„iX, «2,А, ■ ••, «п-1,А, <а)Г мож но принять в виде (4.22). В этом случае получается точное решение системы (4.84), если на каждом шаге проводить нормировку в соответствии с формулами (4.21).

Очевидно, что размерность системы (4.84) может быть уменьшена, если представить ее в виде

uf \ \ - ( 2- h2k « )u$ +uf+u +ft2 [fjA aiVj+bj uj+ gj)+ kt] * !?= °.

Здесь учтены условия (4.82), поэтому Ш(>Л= щп,а = 0.

После решения системы (4.85) относительно производных и нор­

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

принимая во внимание граничные условия (4.82), при помощи програм­ мы DA1EXP.

В качестве примера рассматривалась задача

dv

— = v + и + 4г sin хх, dt

При уменьшении начального шага интегрирования точность вычи­ слений повышается. Вычисленные функции v(x,t), u(x,t) симметричны

140 Глава 4. Дифференциально-алгебраические уравнения

по х относительно середины отрезка [0,1] изменения х. При увеличении числа точек разбиения п (рассматривалось n = II и п = 16) отличие между функциями v(x, t) и tt(x, t) уменьшалось.

Так, при начальном шаге интегрирования Адо = 0,05 и при п = 11 ошибки Д„ и Д0( равные разностям между численными значениями функции v и и и аналитическими, подсчитанными по формулам (4.91), в момент t = 1 были равны Д„ = 0,163 и Д„ = 0,202. При Адо = 0,025 Д„ = 0,091 и Дм = 0,128. При Ад0 = 0,0125 Д„ = 0,0386 и Д„ = 0,0755. При п = 16 и Адо = 0,025 Д„ = 0,075 и Д„ = 0,094.

Как и следовало ожидать, ббльшая ошибка накапливается в пере­ менной, которая не дифференцируется, т.е. в функции и.