оператором DIMENSION А (4, 4), а фактический размер равен в лай ной задаче трем. Тогда девять элементов матрицы будут размещець в 16 зарезервированных под массив А последовательных ячейках памяти следующим образом;

Ои «21

1-Й столбец 2-й Столбец 3-й столбец 4-й столбец

Пример показывает, что если фактический размер матрицы ..VI меньше указанного макси.мального размера МО в операторе DIMENSION, то элементы в ячейках памяти располагаются не подряд. Порядковый номер ячейки К. соответствующей элементу А (I, .J), определяется формулой К = (-1 - I) * МО + I. где МО - максимально возможная длина столбца.

Однако, как отмечалось выше, в вьоываемой стандартной подпрограмме матрица Л представлена в виде массива переменной длины, а элементы матрицы должны быть расположены в массиве А подряд без пропусков ячеек. Поэтому при обращении к стандартной подпрограмме, проводимом при М = 3, будут использованы числа, содержащиеся в первых девяти (М * М = 3 3 = 9) ячейках, зарезервированных под массив А в вызывающей программе. Очевидно, что эти числа не соответствуют коэффициентам построенной в вызывающей программе матрицы А (3, 3). Это происходит потому, что при описании матрицы в подпрограмме в виде массива А (Д\, М) переменной длины выбор номера К ячейки памяти, соответствующей элементу А (I, J), производится иа основе формулы К ~- (J - I) * М 4- I, где М - фактическая длина столбца, указанная при обращении к подпрограмме. Таким образом, при М ?ь д\0 матрица А будет передаваться в подпрограмму неправильно.

Чтобы избежать рассмотренного несоответствия, следует перейти к векторной форме хранения матриц, при которой каждый столбец располагается в памяти непосредственно за предыдущим. В этом случае все М * М элементов расположатся в памяти подряд, а последующие (МО * МО - М * М) ячеек из общего числа зарезервированных останутся пустыми. Для нашего примера элементы матрицы А при М - 3 будут при использовании векторной формы хранения расположены так:

i-й столбец 2-й столбец 3-й столбец свободная область

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

Можно сформировать матрицу в вызывающей программе в оной форме, а перед обращением к стандартной подпрограмме vniecTBHTb преобразование двумерной формы в векторную, т. е.

I зовать «сжатие» области памяти, в которой находятся элементы Р Эта операция может быть выполнена с помощью специаль-ойстандартной подпрограммы ARRAY, которая осуществляет пре-обоазование одномерного массива матрицы в векторной форме в двумерный и наоборот. Обращение к этой подпрограмме имеет вид

CALL ARRAY (MODE, N, М, NO, МО, В, А),

где MODE - код преобразования (I - одномерный массив в двумерный, 2 - двумерный в одномерный); N. М - фактическое число строк и столбцов матрицы; N0, МО - число строк и столбцов, определенных для матрицы в операторе DIMENSION; В - одномерный массив длиной N * М (вектор); А - двумерный массив (матрица).

Вектор В и матрица А могут быть в памяти совмещены. Например, рассмотренное выше преобразование матрицы А (4, 4) в векторную форму при М = 3 реализуется оператором

CALL ARRAY (2, 3, 3, 4. 4, А, А).

2, Д10ЖИ0 формировать матрицу в подпрограмме, входными параметрами которой являются массив А, число строк N и число столбцов М. При этом в подпрограмме следует описать массив А как двумерный массив переменного размера с помощью оператора DIMENSION А (N, М). Тогда при формировании матрицы все ее элементы будут записаны подряд без «пробелов» в памяти, как это н требует векторная форма. В головной вызываюшей программе надо лишь определить максимально возможную длину массива А, причем его можно омнсать как одномерный, например DIMENSION А (1000). Описанный путь реализован в большинстве приводимых в дайной книге программ.

3. Можно сразу формировать матрицу в виде одномерного массива, в котором индекс К. соответствующий индексам I, J матрицы с N строками и М столбцами вычисляется по формуле К (J - I) * N +

Если матрица имеет специальный вид. например является сим--четричной или ленточной, то для работы со стандартными подпрограм-ами коэффициенты матрицы всегда должны быть записаны подряд в одномерный массив в последовательности, зависящей от вида матри-ДЬ( и используемой подпрограммы. Примеры таких способов записи матриц будут рассмотрены ниже.

Перейдем к краткому описанию ряда наиболее важных подпро-1"Рамм.

Характеристика стандартных программ. Подпрограмма обращения пгрицы MINV реализует вычисление обратной матрицы А- ме-

годом Гаусса-Жордана [2. Вызов этой подпрограммы производится так:

CALL MINV(A, N, D, L, M),

где A - входная матрица размером N X N, которая иа выходе заменяется обратной матрицей; D - определитель матрицы, вычисляемый в подпрограмме; L и М - целые рабочие массивы длиной N, которые следует описать в вызывающей программе.

Подпрограмма решения системы линейных уравнений обиего вида GELG реализует решение методо.ч последовательного исключения Гаусса с выбором главного элемента. Эта и последующие подпрограммы предусматривают возможность решения N систем с одной и той же матрицей А, нос различными столбцами правых частей В. Для этого правые части задаются как матрица размером М / N. а N векторов решений также расположены в одном массиве последовательно по М элементов. Такая возможность реализована с целью экономии машинного времени, поскольку в случае N отдельных обращений к подпрограмме с разными правыми частями В над матрицей А будут (гроизводиться одни и те же операции исключения неизвестных. Обра-пеиие к подпрограмме имеет вид

CALL GELG (В. А, М, N, EPS, lER),

где В - матрица правых частей размером М х N, которая на выходе заменяется .матрицей решений; А - матрица коэффициентов системы уравнений размером М X М (в ходе вычислений эта матрица «портится»); М - число уравнений; N - число правых частей; EPS - входная константа, используемая при контроле потери точности; IER-выходной параметр ошибки.

В программе в процессе реализации процедуры исключения предусмотрен контроль потери точности. Для этого на каждом шаге k исключения главный элемент аи, на который производится деление в формулах (1,11), (L14), сравнивается с величиной EPS * аях. тде "max - макси.чальный но модулю элемент матрицы. Если абсолютная величина главного элемента оказывается меньше, чем это произведение, то параметр ошибки принимает значение {к + !) Обычно (адают EPS = (10~ -f- Ю-**). Выходной параметр ошибки IER = " i в случае, когда на каком-то шаге исключения главный эле-MCfiT равен нулю и, следовательно, решение получено быть не может. Значение IER = О свидетельствует об отсутствии ошибок. Параметр IER следует анализировать после обращения к подпрограмме GELG, и в случае IER О предусматривать соответствующее сообщение.

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

рации над нулевыми коэффициентами, что сокращает время решения системы.

Программа GELS реализует решение системы линейных урав-пений с симметричной матрицей коэффициентов. Обращение к ней

имеет вид

CALL GELS (В, А, М, N, EPS, 1ER, AUX).

Здесь А - массив, в который последовательно по столонам загш-саиа верхняя треугольная часть симметричной матрицы, л. е. этот массив содержит М * (М 1)/2 элементов в такой последователь-

ности: Йп. 121 22,

fljAf, .... ймм\ AUX - рабочий

О о Ог О

Рис !.3

массив длиной (М - 1), который должен быть описан в вызывающей программе; остальные параметры по смыслу полностью совпадают с одноименными параметрами подпрограммы GELG.

При формировании массива А коэффициент матрицы (при i < /) записывается в элемент массива А (К), но.чер которого определяется по формуле К = J * - - 1)/2+ L

Решение методом Гаусса систем линейных уравнений с ленточной матрицей реализует [юдпрограм.ча GELB. Существенной особенностью этой подпрограммы является способ записи коэффициентов матрицы. Для его пояснения рассмотрим матрицу с числом верхних диагоналей 2 и числом нижних диагоналей 1, показанную на рис. [.3. В массив А следует записать/70 ст/?ол:ал( все коэффициенты, лежащие в пределах ленты матрицы, т. е. для данного при.мера массив А должен содержать такую последовательность коэффициентов: йц, ai, й,з, Й22. 22., 0,3, Й24, Й54, Й55. Подпрограмма вызывается с помощью следующего обращения:

CALL GELB (В, А, М, N, MUD, MLD, EPS, IER),

гдеМ ~~ число уравнений; N - число правых частей; MUD верхних диагоналей, т. е. лежащих выше главной; MLD нижних диагоналей, т. е. лежащих ниже главной; А - матрица, записанная по строкам в одномерный массив, под который в вызывающей Программе должно быть отведено (М * МС ~ ML * (ML + 1)/2) 1; здесь ML = (МС - MLD - 1), МС = min (М, 1 + MUD J" MLD) - ширина ленты, остальные параметры совпадают с одно-

енными параметрами подпрограммы GELG. gycVn 3 н 4 будут рассмотрены еще две подпрограммы: МГНр ~ решения систем с трехдиагональиой матрицей и ° ~-для решения систем уравнений с ленточной симметричной Р"ией методом квадратного корня.

число число

§ 1.4. ПРОГРАММНАЯ РЕАПИЗАЦИЯ РАСЧЕТА СТАЦИОНАРНЫХ СРЕДНИХ ТЕМПЕРАТУР

При составлении программы для расчета стационарных средних температур в системе тел н потоков теплоносителей используем систему уравнений (1.2)-(1.5). Приравнивая производные по времени нулю и подстаатяя выражения (1.5) для средней температуры теплоносителя Vi в уравнения баланса (1-2), (КЗ), получим систему алгебраических уравнений относительно температур Ti (i 1..... N.,)

UT- Т" {I = 1, л/ж)- ДЛя составления программы целесообразно сгруппировать коэффициенты при неизвестных и записать эту систему в следующей форме:

уравнения для твердых тел ii ~ 1, Nj)

,, /= 1 (=1 j /= 1 i-1

уравнения для объемов с теплоносителем (/ = 1, .... N

1=1 \ - I

ЦВЫК

(Ь27)

уравнения для определения температур теплоносителя и а вхо-де Ul

с \ Ж

2 c,,G„„+2c,G,, 2 c.n0iVT-- 2 0,G,,e,;(].28)

здесь в уравиеиии (1.27) С, 2 G/ + 2G,t,.

Представим систему (1.26)-{1.28) в матричной форме

где А - матрица размером N х N, yV = yV -- 2 * Л,; X - столбец неизвестных температур, содержащий Ы членов; В - столбец из jV свободных членов.

Примем следующий порядок нумерации jV неизвестных температур: первые членов столбца X соответствуют температурам тел Л.....Тл.. далее парами следуют 2 * N. температур теплоносителей

на выходе и на входе Of, Будем использо-

системы: первые yV, далее следуют парами

также следующую нумерацию уравиеиии tuPHHfi совпадают с уравнениями (1.26). да RTypaBuenm 0.27), 0-28).

Приводимая инже программа расчета средних температур

, 1 4) оформлена в виде подпрограммы, параметрами которой вляются все необходимые исходные данные н рабочие массивы, писанные в комментариях к тексту. На основе принятой иу.мераи.ии лясистемы(1.26)-(1.28) в подпрограмме SYSTT произюдится формирование массивов А и Т, соответствующих матрице и столбцу свободных членов, а затем с помощью стаидартион подпрограммы GELG {см. § 1-3) решается система уравнений и находится столбец неизвестных. Отметим, что после решения системы найденные температуры размещаются в том же массиве Т, в котором находились свободные члены. В уравнениях баланса в общем случае могут быть связаны температуры Г,-, Tj, Ui с любыми ио.мерами /, /, /, поэтому матрица А является матрицей общего вида и для решения системы используется подпрограмма GELG.

Рассмотрим процедуру формирования матрицы А н столбца Т. Сначала двумерный массив А и одномерный массив Т обнуляются, а затем производится расчет их ненулевых элементов путем последовательного суммирования отдельных членов, входящих в формулы (1.26)-(1.28). Организация этой процедуры суммирования зависит от используемого способа описания теплового взаимодействия между элементами системы.

Для примера поясним способ описания теплового взаимодействия между телами. В случае, когда каждое тело взаимодействует со всеми остальными, все а]} отличны от нуля, и для описания тепловых связей можно использовать двумерный массив, элементами которого являются о]}. Однако обычно на практике число взаимодействующих пар тел значительно меньше максимально возможного значения, и поэтому большинство тепловых проводимостей равно нулю. С целью экономии машинной памяти и сокращения объема исходных данных целесообразно приводить информацию только об отличных от нуля Oi/- В приведенной программе эта информация задается с помощью двух одномерных массивов IJ и S1J. В массиве U длиной 2 * N1J, где NIJ - число связей между телами, парами записаны номера взаимодействующих тел: элементы [J {2 * М - I), [J {2 * М) указывают номера тел, участвующих в М-м взаи.модействни. Массив SIJ

длиной NIJ содержит значения тепловых проводимостей а", соответствующих этим тепловым связям.

Совершенно аналогичным образом представлена информация о тепловых связях ajf между телами и теплоносителями с помощью массивов IL и SIL, а также о связях gJI между телами и средами с ло.чощью массивов (К и SIK-

Движение потоков теплоносителей описывается с помощью массива ML длиной 2 * NG и массива GC длиной NG, где NG - число по-