Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена (Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена.djvu), страница 8

В настоящее время широкое распространение получили четырех- этапные схемы Рунге — Кутта, имеющие четвертый порядок аппроксимации и называемые в связи с этим схемами Рунге — Кутта четвертого порядка. Наиболее употребительная из них имеет вид: и~+' — ит + (Л+2~,+21,+я 6 1 =7(, '), (.=ф~+ —, г+ — ~,), (1.48) Лт . Ьт 2 2 Лт . Лт И=[*[ з+ —, кг+ — )з[, ~, =1(т~+„и~+М,).

2 2 Здесь оценки производной вычисляются 4 раза: в точке тл дважды в точке т; + Лт!2 и в точке тььь а для вычисления интеграла используется квадратурная формула Симпсона. Схемы Рунге — Кутта (1.47) и (1.48) явные и обладают условной Устойчивостью. Например, у схемы четвертого порядка условие устойчивости, проверяемое на модельной задаче (1.38), выполняется только при Лт( 2,8/т0. Рассмотренные два вида схем обладают общей чертой — при определении значения и~+' в узле ть, используется только значение и' в предыдущем узле тл Такие схемы называют одношаговыми.

Ясно, что в общем случае возможно для определения и!+' использовать не только и!, но и ит-', ит-' и т. д. В этом случае мы приходим к многошаговым методам, среди которых наибольшее распространение получили линейные. Линейные многошаговые методы. При построении многошаговых схем, как и всхемах Рунге — Кутта, будем исходить изравенства (1 46). Однако для приближенного вычисления интеграла применим другой прием, а именно на основе значений функции 7 (т, и) в й Зак. 2217 33 точках т! тт-!, ..., т; х+! построим интерполяционный полипом, который принимает в этих точках уже известные значения ? (тп иг), ... !'(т~-х+!, и' 4+!). Далее, используя продолжение построенного полинома за точку т (экстраполяцию на отрезок [т;, т,+,)), можно вычислить интеграл в (1.45), и получить выражение для ш+!, в котором используются найденные ранее значения и!', ..., и!'-"+'.

и!+' — иг = ==- !!т ч'. [[„,~(т! эи и!' — +'), (1.49) т=1 где р — постоянные коэффициенты, зависящие от числа точек интерполяции. Например, если на основе значеРис. 1и ний иг и и!'-' построить линейную функцию (рис. !.5), продолжить ее на интервал [т,, т,э,! и провести интегрирование выражения (!.45), то получим следующую двухшаговую схему: и!+ ! — иг.= ~ — ) (тп и!) — — [(тз „и!' !)~ Лт. (1.50) Г 3 . ! ! 2 2 Схемы вида (1.49) явные. Однако несложно получить и неявные схемы.

Для этого следует использовать полипом, проходящий не только через известные точки иг, ..., и!-ь!!, но и через неизвестную точку иг!!. Тогда, интегрируя, получим уравнение неявной схемы относительно и!+! и!+! — и1 =Лт 'с', р ~(тт +ь ит +'). (1.51) Достаточно часто применяют следующую явную схему при й = =4: иГ+ ' — ит = — [5ЦI — 59~~ — !+ 3?г! — ' — 91! — з[ (1.52) 24 и неявную при й = 3 и!+ ! — и! = — [9~l+ ! + 19!! — 5!!' — ! +~l — з[; (1.53) 24 зДесь ~/ =- 1 (тд и!). При получении схем (1.49), (1.51) мы исходили из соотношения (!.45).

Если же брать за основу равенство ит+! — иг —" = ~ 1(т, и) бт, 0 < л < й — 1, (1.54) $ т-и 34 то можно аналогично получить схему з ш+~ иг — Ьт ~~Р р г(т~ + иг +~) Таким образом, в общем случае линейную й-шаговую схему можно записать в виде я=О В уравнении (1.55) предполагается, что значение и/+'вычисляется на основе л известных значений и~, ..., иг "+'. Если в (!.55) Ц = О, то уравнение разрешается относительно из+' — схема явная. Если (), чь О, то иг+' входит и в левую и в правую части, получается в общем случае нелинейное уравнение относительно и~+', схема— неявная.

Остановимся на ряде особенностей многошаговых (при й ) 2) методов по сравнению с одношаговыми (Й = 1). За счет введения в разностную схему значений функции ) (т, и) в я тачках, предшествующих искомой (1 + 1)-й точке, удается повысить порядок аппроксимации. Похожий прием использовался для повышения порядка аппроксимации в методе Рунге †Кут, но там вычисление значений )(т, и) проводилось в точках интервала (ть т~+,). С точки зрения сокращения затрат машинного времени для одного шага новый способ предпочтительнее, поскольку в схемах Рунге— Кутта вычисленные на промежутке (ть тз+,) значения )(т, и) не будут использованы на следующем шаге от т;+, до т з+„а основные затраты времени связаны с вычислением этих значений. В многошаговом же методе при вычислении иг+' мы не сможем использовать только значение ) (т; „+„иг-"+') в наиболее удаленной точке, участвовавшей в определении значения и~+'.

Остальные значения функции 1 в точках т,, ..., т; „+, можно использовать и при вычислении ш+'. Однако в целом сопоставление затрат машинного времени нужно проводить, учитывая общее число шагов 1, необходимое для достижения заданной погрешности. В отличие от одношаговых методов многошаговые не являются самостартующими. действительно, в одношаговых методах вычисления можно начать с и~ = Т„получить на его основе и', затем на основе и' найти и~ и т. д. При использовании же многошаговых методов при ~ 2 вычисления могут начаться только лишь с из при условии, что каким-то образом определены и', и', ..., и"-'.

Эти стартовые значения находят с помощью какого-либо одношагового метода, в качестве которого обычно используют одну из схем Рунге — Кутта. С точки зрения ограничения на шаг Лт, связанного с требованием наличия у схем устойчивости, явные многошаговые методы не имеют "Реимуществ по сравнению с явными методами Рунге — Кутта. 2' ЗЬ Рассмотрим особенности реализации неявных многошаговых методов. Перепишем схему (1.55) в виде и'+ ' = — г" ( ~+ и'+')+Ф~ (1.56) аь где величина й~ зависит от уже найденных значений и/, ..., иг-ь+' и известна. Для решения уравнения (1.56) применяют два метода: метод простой итерации и метод Ньютона. Рассмотрим первый метод.

В этом случае ие+' находим с помощью итерационного процесса (их+')' = — '~(ту+„(иг+')' — ')+йн з=1, 2, ..., (1.57) Яь где (иг+')ь — начальное приближение к искомому значению функции. Лля быстрой сходимости итерационного процесса нужно удачно выбрать начальное приближение (ие+')ь. Обычно его получают с помощью явной многошаговой формулы. Тогда в целом алгоритм расчета на каждом шаге выглядит так: сначала предсказывается значение (иг+')', а затем проводится одно или несколько уточнений начального значения по формуле (1.57), т. е.

получается предсказывающе-исправляющий метод, который часто называют методом предиктар — корректор. Используют два способа реализации этого метода. В первом случае итерации проводят до тех пор, пока два последовательно полученных значения (и~ь')'"', (ит+')' не станут достаточно близки. Тогда число итераций может меняться от шага к шагу. Во втором случае на каждом шаге исправляющая формула применяется фиксированное число раз. Обычно на практике метод предиктор — корректор применяется с числом итераций не более трех. Наиболее часто из линейных й-шаговых схем используют схемы (1.49), (1.51), называемые схемами Адамса.

Можно доказать, что явная схема Адамса имеет порядок аппроксимации равный й, а неявная — (й + !). При использовании метода предиктор — корректор обычно применяют предсказывающую и исправляющую схемы одного порядка точности. В частности, широко применяется метод предиктор — корректор со схемами Адамса четвертого порядка, в котором предсказание делается по формуле (1.52), а уточнение — по (1.53). Выбор шага интегрирования н оценка погрешности численного решения. Обычно при реализации на ЭВМ большинства численных схем решения обыкновенных дифференциальных уравнений предусматривается автоматический выбор величины шага Лт для обеспечения определенной погрешности расчета.

Этот выбор основан на оценке локальной погрешности численного решения на шаге, т. е. погрешности численного решения в точке те+„оцениваемой в предположении, что в начале шага в момент времени т; значение искомой функции было известно точно. 36 Для одношаговых методов существуют два основных способа оценки локальной погрешности на шаге. При первом — интегрирование от точки т; до точки тг+, осуществляется по схеме одного и того же порядка точности р дважды: с шагом Лт и с шагом Лт/2. Можно доказать, что в этом случае для локальной погрешности е!+„' Решения, найденного с шагом Лт/2, справедлива следующая приближенная оценка: з~'+ ' = (и/+ ' (Лт/2) — и/+ ' (Лт)1/(2г — 1), (1.58) где и/+' (Лт/2), ш+'(Лт) — решения в точке туэ м найденные на основе одного н того же решения иг в точке т; с шагами Лт/2 и Лт соответственно.

При втором способе интегрирование на отрезке (ть тэ+,1 также выполняется дважды: с помощью схемы порядка точности р и схемы порядка точности (р + ! ) при одном и том же шаге Лт, и находятся соответственно два решения в точке т;+,. иг+' (р) и и/+' (р + 1). В этом случае локальная погрешность оценивается по формуле е/+' ж и/+'(Р) — иг+' (р+1). (1.59) Обычно при расчетах указывается допустимое значение абсолютной или относительной локальной погрешности, в соответствии с которым автоматически подбирается величина шага. Например, при первом способе оценки погрешности это делается так.

В ходе расчетов сначала из некоторой точки т; делается шаг Лт, затем два шага длиной Лт/2, и на основе двух решений иг+' (Лт) и и/~' (Лт/2) оценивается погрешность е(+,' в точке т;+, по формуле (1.58). Если погрешность не превосходит допустимого значения, то и~+' (Л г/2 ) принимается в качестве найденного значения сеточной функции, и из точки тг+, продолжается расчет с шагами Лт и Лт/2. Если оценка локальной погрешности превосходит допустимое значение, то из точки тг делают два шага длиной Лт/4 и оценивают погрешность в точке тг+ Лт/2. При достижении требуемого уровня погрешности расчет из точки т; + Лт/2 продолжается с шагами Лт/2 и Лт/4. В противном случае проводится дальнейшее измельчение шага.

В том случае, когда при первой проверке погрешность оказывается существенно ниже допустимого значения, обычно шаг увеличивают путем его удвоения. Описанные приемы оценки локальной погрешности внешне выглядят очень просто. Сложности начинаются при их программной реализации с учетом минимизации затрат машинного времени. Алгоритм стремятся построить таким образом, чтобы при проведении вспомогательных шагов использовать по возможности большее число уже вычисленных значений функции / (т, и) и находить сравнительно небольшое число ее новых значений.

При использовании многошаговых схем по методу предиктор— корректор оценку локальной погрешности на шаге обычно удается 37 выразить через вычисляемые при реализации схемы предсказываемое и исправленное значения искомой функции. Таким образом, в этом случае для оценки погрешности не требуется дополнительных вычислений функции ~ (т, и), что является преимуществом по сравнению с одношаговыми методами. Однако у многошаговых методов сложнее реализуется процедура изменения шага. При уменьшении шага необходимо либо путем интерполяции определить значения искомой функции в новых предшествующих точках, либо получать по одношаговой схеме новые «разгонные» точки.

При любой стратегии, основанной на оценке локальной погрешности, не учитывается накопление погрешности в ходе всего расчета и, следовательно, фактическая погрешность численного решения остается неизвестной. Более надежный путь оценки погрешности — сравнение двух решений (иг)' и (и~)", полученных во всей расчетной области [О, т „„[ при двух сетках разной густоты — с шагами Лт' и Лт". В этом случае полная погрешность ег = Тг — (иг)" оценивается по следующей формуле: ег ж [(и~)" — (и1)'УЯМ'[Лт')» 1) (1.60) где р — порядок точности схемы.