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

— разработаны специальные кубатурные формулы вида /=-~,с;1!х;), (6. 14) где с; — некоторые коэффициенты; /(х,) — значения функций в определенных точках, лежащих в области интегрирования. Коэффициенты с; и расположение точек выбираются так, чтобы на каком- либо классе функций погрешность формулы (6.14) была бы минимальной !2!. Применение подобных формул ограничено требованиями к форме областей.

Последовательное интегрирование. Рассмотрим этот метод на примере двумерного интеграла по прямоугольнику. Этот интеграл можно представить в виде аь 1 = ) ~ / (х, у) д хг(у = ~ у (У) «)У, (6.15) с а с При использовании (6.15), (6.! 6) сначала с помощью какой-либо квадратурной формулы проводят численное интегрирование по «горизонтальным прямым» и определяют значения функции Р (у) в точках разбиения по переменной у. Затем на основе этих значений также по квадратурной формуле вычисляется окончательное значение двумерного интеграла.

В принципе для разных направлений квадратурные формулы могут различаться. Однако обычно используют квадратурные формулы одного порядка точности. Метод последовательного интегрирования можно применять и для областей сложной формы, но в этом случае при его программной реализации возникают ббльшие по сравнению с методом ячеек сложности. Метод Монте-Карло. Хотя этот метод целесообразно применять именно для многомерных интегралов и, более того, его эффективность возрастает с увеличением размерности, для простоты сначала рассмотрим основную идею метода на примере одномерного интеграла ь 1 =- ~1 (х) дх.

а Представим этот интеграл в виде ! = 1 — р (х) цх, г 1(х) ,) р (х) а (6.17) где р (х) — некоторая функция, обладающая свойством ~р(х) 4х= 1, а (6.! 8) Е(Х) =-~хр(х) 4х, а (6. 19) а математическое ожидание Л вЂ” Е (Л) — по формуле ь Е(Л) =~ф(х) р(х) бх. (6.

20) 188 на которую пока никаких дополнительных ограничений не накладываем. Напомним ряд формул из теории вероятностей. Рассмотрим непрерывную случайную величину Х, принимающую значения только из промежутка [а, Ь) и имеющую функцию плотности распределения вероятности р (х), а также вторую случайную величину Л, связанную с Х функциональной зависимостью Л = ф (Х).

Математическое ожидание величины Х вЂ” Е (Х) рассчитывается по формуле ь Сравнивая формулы (6.17) н (6.20), видим, что интеграл (6.!7) можно трактовать как математическое ожидание непрерывной случайной величины Л, связанной с Х зависимостью (6.21) Такая трактовка позволяет указать оригинальный способ вычисления интеграла (6.17). Вспомним, что в математической статистике математическое ожидание случайной величины оценивается по среднеарифметическому значению из совокупности результатов ее наблюдений, которые берутся из эксперимента. В методе МонтеКарло применяется такая же оценка, но результаты наблюдений берут не из эксперимента, а получают путем статистического моделирования на ЭВМ.

Для этого реализуется специальная процедура генерирования последовательности значений независимых реализаций х„..., хя случайной величины Х с функцией плотности распределения р (к), Имея набор х,, ..., хл, рассчитывают значения 3,„..., Лл реализаций случайной величины Л: )!! =-1(х!)/р (х;) и далее находят оценку математического ожидания Л по формуле Е(Л) = — ~ )!!=1= — '~' )(х)(р(х!) =-7 = I.

(622) й!, 4 й! != ! 1=! Эту оценку считают приближенным значением интеграла (6.17) Процедура получения значений независимых реализаций х; слу чайной величины Х строится на основе стандартных программ для генерации псевдослучайных последовательностей. Определим погрешность формулы (6.22) для вычисления интеграла. Напомним, что погрешность оценки математического ожидания пропорциональна ее среднему квадратическому отклонению, которое убывает пропорционально 1!)~ !У.

Например, среднее квадратическое отклонение а-„выборочного среднего, определенного по выборке Л„ ..., Лл из нормального распределения, равно о-„ = = о!)/ !1!', где о — среднее квадратическое отклонение одного результата наблюдения. Для среднего квадратического отклонения оценки (6.22) ау справедливо аналогичное равенство (6. 23) о! = о (Л)Д' Ж. Таким образом, погрешность формулы (6.22) можно уменьшать не только за счет увеличения У, но и за счет понижения а (Л).

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

В этом случае любую координату можно моделировать независимо от остальных аналогично одномерному случаю. Остановимся для примера на расчете рассмотренного выше углового коэффициента рдв со(см. рис. 6.6) методом Монте-Карло. В качестве случайного вектора Х здесь выступает совокупность двух значений координат (х, у).

Для получения простейшей функции плотности распределения р (х, у) можно принять, что компоненты х и у статистически независимы и равномерно распределены на соответствующих интервалах своего изменения [а, Ь1 и [с, 4: р(х, у) =-1/[(Ь вЂ” а) (Й вЂ” с)1. Тогда при вычислении ч/ла со последовательно генерируются пары значений координат (хо у;), для каждой пары рассчитывается значение подынтегральной функции /(хо у;) и, наконец, после получения й/ реализаций определяется оценка срлв ср по формуле ю=-! Значения хь у; генерируются независимо друг от друга с помощью двух обращений к стандартной подпрограмме, рассматриваемой далее в 8 6.3.

При расчете / (х;, у;) остаются все упомянутые выше трудности, обусловленные наличием затененности. Мы рассмотрели простейший вариант реализации метода МонтеКарло для расчета интегралов. Существует ряд приемов повышения его точности, с которыми можно познакомиться в !2, 281. Проведем сопоставление скоростей сходимости методов ячеек и Монте-Карло. Из (6.23) вытекает, что погрешность определения многомерного интеграла с помощью метода Монте-Карло убывает пропорционально 1/)/ й/, где Л/ — число многомерных точек. Причем скорость сходимости не зависит от размерности интеграла. В методе ячеек„применяемом для кусочно-аналитических подынтегральных функций, которые, как было указано, часто встречаются при расчете угловых коэффициентов, скорость сходимости пропорциональна 1/и, где и — число отрезков разбиения по каждой координате. Поскольку в методе ячеек для расчета т-мерного интеграла необходимо рассчитывать й/ = и" многомерных точек, погрешность численного интегрирования в этом случае будет иметь порядок 1/$/ й/.

Отсюда вытекает, что для двумерных интегралов (т=2) 188 оба метода дают одинаковую скорость сходнмости, а начиная с трехмерных интегралов скорость сходимости метода Монте-Карло выше. Это позволяет рекомендовать его для численного расчета угловых коэффициентов. $ Ь.З. РАСЧЕТ УГЛОВЫХ КОЭФФИЦИЕНТОВ МЕТОДОМ СТАТИСТИЧЕСКОЙ ИМИТАЦИИ Основная идея метода. Имитация является одной из разновидностей метода Монте-Карло. Общую идею и схему применения этого метода несколько упрощенно можно сформулировать следующим образом. Для решаемой задачи, кото "".

состоит в определении некоторого параметра, конструируется случайная величина, распределение которой зависит от этого искомого параметра. С помощью ЭВМ проводится моделирование построенной случайной величины, в результате которого находится набор ее реализаций. Далее по этому набору вычисляется статистическая оценка искомого параметра, которая и принимается за решение исходной задачи. В $ 6.2 мы разобрали один из вариантов применения метода Монте-Карло, при котором задаче определения интеграла (6.1!) ставилась в соответствие случайная величина Л, математическое ожидание которой равнялось искомому интегралу.

Этот вариант можно назвать У статистическим интегрированием. Однако О,ь возможны и другие приемы. Ниже мы рассмотрим один из наиболее распространенных подходов, основанный на статистической имитации процесса переноса излучения. Проиллюстрируем его сначала на примере расчета углового коэффициента ~ры.

При подборе подлежащего имитации Рис. 6.7 соответствующего случайного эксперимента надо учитывать в данном случае два обстоятельства. Во-первых, все элементарные площадки 65о принадлежащие поверхности 5о дают одинаковый вклад в ее общее излучение. Во-вторых, для диффузно излучающей поверхности число фотонов, вылетающих с любой элементарной площадки в направлениях, задаваемых интервалами изменения азимутального угла (ф, ф + бф! н угла отклонения от нормали [О, 0+ ИО! (рис.

6.7), пропорционально произведению соз 0 з)п 0 бйбф. Это вытекает из закона Ламберта (0, 30), согласно которому интенсивность излучения 1 рассчитывается по формуле 6Р е (6.24) 6065 сои 6 и 189 где е — полусферическая плотность потока излучения; 0 — угол отклонения от нормали к поверхности; бР— поток, излучаемый в элементарном телесном угле дь) = яп Обйдф Указанные обстоятельства позволяют рассмотреть следующий случайный эксперимент.

На поверхности 5; случайным образом по равномерному распределению выбирается точка. Для этой точки также случайным образом по равномерному распределению для азимутального угла ф и распределению с плотностью вероятности, пропорциональной з(п 9 соз О, для полярного угла 0 выбирается направление распространения «порции» излучения с энергией Лалее рассматривается следующая случайная величина Л: если луч, проведенный в выбранном направлении из выбранной точки, попадает на поверхность 5ь то величина Л принимает значение ~, в противном случае — нулевое значение. Очевидно, что математическое ожидание введенной случайной величины равно Е (Л) =с~;~ Я, и тогда угловой коэффициент можно определить так: В! — — Е ( Л)(а.