Моделирование миграции и осаждения примеси на основе конвекционно-диффузионной модели

Выпускная квалификационная (дипломная) работа

Моделирование миграции и осаждения примеси на основе конвекционно-диффузионной модели

конвекция уравнение программа двумерный Современная цивилизация осуществляет невиданное давление на природу. Загрязнение природной среды промышленными выбросами оказывает вредное действие на людей, животных, растения, почву, здания и сооружения, снижает прозрачность атмосферы, повышает влажность воздуха, увеличивает число дней с туманами, уменьшает видимость, вызывает коррозию металлических изделий. Под загрязнением окружающей среды следует понимать изменение свойств среды (химических, механических, физических, биологических и связанных с ними информационных), происходящие в результате естественных или искусственных процессов и приводящие к ухудшению функций среды по отношению к любому биологическому или технологическому объекту. Используя различные элементы окружающей среды в своей деятельности, человек изменяет её качество. Часто эти изменения выражаются в неблагоприятной форме загрязнения. В связи с этим неоценима важность и ценность математического моделирования в экологии, позволяющего уже на ранних этапах, оценить влияние, оказываемое на внешнюю среду.

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

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

Этапы эволюции методов исследования отображены, например, в работе М. Е. Берлянда. Применение моделей атмосферной диффузии к проблемам охраны окружающей среды получило значительное развитие в работах Г. И. Марчука [2 — 4].

В результате выполнения дипломной работы рассмотрены статические двумерные модели распространения загрязняющих примесей на фоне атмосферных процессов. Написана программа расчета и визуализации численных результатов на основе модели переноса однокомпонентных газовых примесей и аэрозолей, выбрасываемых неподвижным точечным источником. Для написания программы использованы следующие программные средства: Compaq Visual Fortran 6.5 и Maple 9.0.

1. Постановка задачи

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

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

1.1 Постановка задачи миграции примеси

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

Рассматриваемое в работе уравнение, описывающее перенос и диффузию однокомпонентной примеси при постоянных метеорологических характеристиках среды и в декартовой системе координат имеет вид

(1.1)

где — величина, характеризующая концентрацию ЗВ,

— коэффициент концентрации,

— вертикальный коэффициент диффузии,

— горизонтальный коэффициент диффузии,

— функция, описывающая источник загрязняющего вещества.

Для описания действия источника можно использоватьфункции Дирака. Тогда функция, где — мощность источника, — координаты источника в данный момент времени.

Рассматривается конечная область рассеяния примеси с границей S при следующих начальных и граничных условиях:

;

,

.

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

.

Дискретные аппроксимации для рассматриваемой задачи строятся на базе вариационного принципа в сочетании с методом расщепления [2 — 4]. При таком подходе осуществление технологически простой реализации дискретной модели на ЭВМ.

1.2 Параметризация

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

Рассмотрим уравнение Уравнение содержит группу параметров, характеризующих процесс распространения примеси. К числу этих параметров относятся:

— скорость в направлении оси Ox;

— скорость в направлении оси Oz[w]=м/с;

— коэффициент горизонтальной диффузии ;

— коэффициент вертикальной диффузии ;

— коэффициент естественного распада ;

— мощность источника выброса ;

— координата сосредоточенного источника ;

— координата сосредоточенного источника ;

— концентрация примеси ;

— ;

—; ;

—; .

Приведению к безразмерному виду подвергаются независимые параметры, затем производные. Задача считается решенной, если все безразмерные параметры выражаются через входные.

Пусть некоторые эталоны длины L, [L]=м, времени T, [T]=C и концентрации .

Тогда, ,, ,, .

Преобразуем уравнение В этом уравнении все слагаемые безразмерные. Вводя обозначения, ,, ,, , уравнение диффузии можно представить Установлением связей между размерными параметрами является предметом теории размерностей.

2. Модельные примеры для одномерного уравнения

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

В качестве модельного примера рассмотрим одномерное линейное уравнение переноса где — концентрация вещества, k=const — скорость переноса, — функция, описывающая источник вещества.

Краевая задача имеет следующий вид:

(2.1)

Рассмотрим конечно-разностный метод. Обозначим как Сформируем на области сеточную область — совокупность узлов на прямоугольной сетке с шагами и .

Обозначим через и точное и приближенное значение функции в узле .

Обозначим как — внутренние узлы сетки, — пограничные узлы.

Тогда дифференциальное уравнение можно переписать следующим образом:

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

Таким образом, разностная схема для задачи (2.1) примет вид

(2.2)

Задача (2.2) аппроксимирует на решениях задачи (2.1) и порядок аппроксимации — первый.

Рассмотрим случай одномерного диффузионного процесса где — концентрация вещества, k=const коэффициент диффузии — функция, описывающая источник вещества.

Краевая задача имеет следующий вид:

(2.3)

В области Соответствующую сеточную область обозначим .

Явная разностная схема имеет вид

(2.4)

Разностная схема (2.4) аппроксимирует задачу (2.3) с погрешностью первого порядка по и второго по. Устойчивость разностной схемы выполняется при. Исходя из порядка аппроксимации и устойчивости получаем, что явная разностная схема (2.4) сходится к решению задачи (2.3).

Неявная разностная схема имеет вид

(2.5)

3. Методы решения двумерной задачи

Численное моделирование проблем распространения загрязняющих веществ проводится чаще всего на основе использования разностных схем.

Запишем поставленную в 1.1 задачу в операторном виде [4], т. е.

в области, где .

Здесь

где в D при t=0.

Для решения этой задачи вводится сеточная область. Чтобы построить приближенное решение, необходимо заменить исходную дифференциальную задачу некоторой конечномерной дискретной задачей, обычно представляющей собой систему линейных алгебраических уравнений (СЛАУ) для компонентов вектора. Одной и той же дифференциальной задаче можно поставить в соответствие множество различных дискретных моделей, однако, далеко не все из них пригодны для практической реализации. Вопросы построения разностных схем подробно рассмотрены в [6, 7].

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

(3.1)

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

С учетом изложенного рассмотрим задачу Коши

(3.2)

где в D при t=0

Рассмотрим простейшие методы аппроксимации задачи (3.2) по времени, полагая, что не зависит от времени. Одна из них — явная схема первого порядка аппроксимации на сетке

(3.3)

где

Неявная схема первого порядка аппроксимации имеет вид

(3.4)

где

Это схемы первого порядка аппроксимации (в предположении, что существуют вторые производные по времени функции).

Схема Кранка — Николсона имеет следующий вид [6]:

(3.5)

где

Эта схема аппроксимирует исходную задачу со вторым порядком по времени.

При этом схема (3.3) будет устойчива при выполнении определенного условия, например, если — симметричная положительно определенная матрица с собственными числами из интервала [с, b], а удовлетворяет соотношениям .

Рассмотрим методы расщепления, применяемые при решении двумерных задач, описывающих процессы различной природы. Рассматриваемое эволюционное уравнение имеет вид (3.1) в. Оператор не зависит от времени и представим в виде при условии, что .

Будем полагать, что эта задача уже редуцирована к разностному виду.

3.1 Метод стабилизации

Рассмотрим так называемый метод стабилизации. Для этого рассмотрим разностную схему решения (3.1) в предположении f=0:

(3.6)

где

Эта схема аппроксимирует исходную задачу со вторым порядком аппроксимации по .

С помощью ряда преобразований приводится к виду

(3.7)

где .

Отсюда видно, что (3.7) при достаточной гладкости решения совпадает со схемой Кранка — Николсона, т. е. имеет второй порядок аппроксимации по .

Эта разностная схема допускает удобную компьютерную реализацию.

В случае неоднородной задачи

(3.8)

где при t=0.

В этом случае задача запишется следующим образом:

(3.9)

где ,.

При условии схема будет обладать вторым порядком аппроксимации по .

3.2 Метод покомпонентного расщепления

В случае зависимости операторов от времени применяется метод покомпонентного расщепления.

Пусть в (3.1) оператор не зависит от времени и представим в виде при условии, что и. Рассмотрим аппроксимацию этих матриц на интервале в форме в предположении, что их элементы имеют достаточную гладкость.

Построим систему разностных уравнений, состоящую из последовательности простейших схем Кранка — Николсона

(3.10)

Система разностных уравнений при исключении вспомогательных функций может быть приведена к одному уравнению

(3.11)

(3.12)

Если, , то при достаточной гладкости элементов этих матриц и решения задачи (3.1) разностная схема (3.10) абсолютно устойчива и аппроксимирует исходное уравнение со вторым порядком по в случае, если и коммутируют, т. е., и с первым, если не коммутируют.

Рассмотрим двуцикличный метод покомпонентного расщепления. Будем аппроксимировать операторы и на интервале. Положим

.

Построим следующие две системы разностных уравнений

(3.13)

(3.14)

Имеем

(3.15)

Двуцикличный метод абсолютно устойчив, а схема (3.15) аппроксимирует исходное уравнение (3.1) со вторым порядком по

Будем искать решение неоднородной задачи с помощью двуцикличного полного расщепления.

Рассмотрим систему разностных уравнений вида (3.13), (3.14), записанных в более удобной форме

(3.16)

где. Разрешая эти уравнений относительно, получим

(3.17)

(3.18)

(3.19)

С помощью разложения по степеням малого параметра придем к соотношению

(3.20)

которое преобразуем к виду

(3.21)

Исключим, используя разложение решения в ряд Тейлора в окрестности точки. С точностью до будем иметь

(3.22)

Производную исключим с помощью соотношения

(3.23)

Подставим (3.23) в (3.22). Тогда

(3.24)

(3.25)

Подставим соотношение (3.25) в (3.24). В результате будем иметь

(3.26)

Очевидно, что уравнение (3.26) аппроксимирует исходное уравнение (3.1) на интервале со вторым порядком по. Таким образом, найдена разностная аппроксимация неоднородного эволюционного уравнения второго порядка с помощью двуцикличного метода.

Если, , то при достаточной гладкости решения, функции, и элементов матриц, система разностных уравнений (3.16) абсолютно устойчива на интервале и аппроксимирует исходное уравнение со вторым порядком по .

4. Описание программной реализации решения двумерной задачи

Для численного решения поставленной задачи разработана программа. Алгоритм основан на использовании выше описанного метода покомпонентного расщепления решения дифференциальных уравнений с частными производными. Программа предназначена для расчёта концентрации загрязняющих веществ. Входными данными являются:

— параметры области решения задачи;

— компоненты вектора скорости воздушных масс в направлении осей соответственно ;

— коэффициенты диффузии в направлении осей соответственно ;

— мощность источника примеси;

— величина, характеризующая взаимодействие примесей с подстилающей поверхностью (л);

— координаты источников, в которых проводятся наблюдения за распространением примеси;

— количество итераций по времени исследования;

— количество шагов по временной переменной;

— количество шагов по пространственным переменным.

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

4.1 Выбор среды реализации

Для написания программы была использована среда Compaq Visual Fortran 6.5. Этот выбор обусловлен тем, что Fortran — один из самых мощных языков программирования, позволяющий работать с типами данных повышенной точности, что очень важно при выполнении математических расчётов. Fortran изначально был создан для научных и численных расчетов и всё его последующее развитие ориентировано прежде всего на подобные приложения.

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

4.2 Описание программы

Все входные параметры находятся в отдельном текстовом файле in.txt, их можно изменять непосредственно в этом файле. Ввод данных осуществляется под управлением именованного списка.

Имя списка ввода в данной программе: input. Таким образом, оператор NAMELIST объявления именованного списка ввода есть в разделе объявлений программной единицы и имеет вид: namelist /input/ L, h, tau, T, n, u, w, mu, nu, M, constX, constZ, lambda где после /input/ идёт перечисление заранее определённых в программе в принадлежности какому-либо типу переменных. При вводе именованного списка оператор ввода ищет в файле in.txt начало списка, которое имеет вид: $input. Перечень принадлежащих именованному списку данных завершается знаком доллара ($). Имена переменных во входном файле при использовании именованного списка ввода должны совпадать с соответствующими именами списка переменных оператора NAMELIST.

Написанная программа реализует схему расщепления по физическим процессам.

Рассмотрим первый этап, который соответствует переносу. Разобьем весь промежуток [0, T] на элементарные. Переносу соответствует следующий оператор:

.

Соответствующий разностный аналог имеют вид где — шаги сетки вдоль соответствующих осей.

Рассмотрим следующий этап, соответствующий диффузии и поглощения субстанции. А именно будем рассматривать оператор

.

Разностный аналог оператора будет иметь вид где

За при этом принимается решение задачи на предыдущем этапе (когда рассматривался только перенос).

В результате получается схема второго порядка аппроксимации по и по .

Для граничных значений получаем конечноразностные формулы следующего вида:

.

Основные вычисления происходят в цикле по временной переменной tc. Вначале описывается удовлетворение граничным условиям. Начальное распределение концентрации примеси записывается в двумерный массив fi. Сначала в этом цикле происходит решение системы уравнений, описывающих явление переноса загрязняющих веществ, затем уравнений, описывающих диффузию примеси. Это делается при помощи метода прогонки (так как матрица этой системы является блочно-диагональной с трехдиагональными блоками, первое и последнее уравнение в блоке аппроксимируют соответствующие граничные условия, поэтому в начале цикла они не учитывались). Промежуточный результат записывается в двумерный массив fipr затем, используя получившиеся данные, переходим к следующей системе и проделываем аналогичные действия. После того как, завершился шаг покомпонентного расщепления, результаты записываются в массив fi, выполняется присваивание массиву fi значений массива fipr, происходит вывод результатов на текущем шаге по времени в файл и управление передаётся на начало цикла. На каждой итерации по времени можно увидеть картину распространения примеси.

Вывод производится в текстовый файл out.txt, а из него нужные данные можно выбрать, например, в файл с именем out1.txt и визуализировать в Maple в готовом файле out1.mws.

5. Анализ численных результатов

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

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

В качестве примера использования метода расщепления было решено несколько модельных задач, в частности решалась задача с постоянным точечным источником мощности 20 в области, где D={0

Рисунок 1 — Поверхность распределения примеси на двадцатом моменте времени Рисунок 2 — Поверхность распределения примеси на двадцатом втором моменте времени Рисунок 3 — Поверхность распределения примеси на двадцатом пятом моменте времени

Заключение

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

На основе изученной модели распространения загрязняющих примесей разработана программа, предоставляющая возможность проведения вычислительных экспериментов и анализа их результатов. Программа написана в среде Compaq Visual Fortran 6.5. Численные результаты выводятся в текстовый файл out.txt, из которого можно построить визуализацию в среде Maple.

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

Список использованных источников

Берлянд, М. Е. Современные проблемы атмосферной диффузии / М. Е. Берлянд. — Ленинград: Гидрометеоиздат, 1975. — 448 с.

Марчук, Г. И. Математическое моделирование в проблеме охраны окружающей среды / Г. И. Марчук. — М.: Наука, 1982. — 320 с.

Марчук, Г. И. Численное решение задач динамики атмосферы и океана / Г. И. Марчук. — Ленинград: Гидрометеоиздат, 1974. — 303 с.

Марчук, Г. И. Методы вычислительной математики / Г. И. Марчук. — М.: Наука, 1989. — 608 с.

Алоян, А. Е. Моделирование динамики и кинетики газовых примесей и аэрозолей в атмосфере / А. Е. Алоян. — М.: Наука, 2008. — 415 с.

Бабешко, В. А. Математическое моделирование экологических процессов распространения загрязняющих веществ / В. А. Бабешко, А. В. Павлова, О. М. Бабешко, О. В. Евдокимова. — Издательско-полиграфический центр Кубанского государственного университета, 2009. — 138 с.

Самарский, А. А. Численные методы решения задач конвекции-диффузии / А. А. Самарский, П. Н. Вабищевич. — М.: URSS: Либроком, 2009. — 246 с.

Бахвалов Н. С. Численные методы / Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. — М.: Лаборатория Базовых Знаний, 2001. — 632 с.

Агошков, В. И. Методы решения задач математической физики / В. И. Агошков, П. Б. Дубовский, В. П. Шутяев. — М.: ФИЗМАТЛИТ, 2002. — 303 с.

Пененко, В.В., Модели и методы для задач охраны окружающей среды / В. В. Пененко, А. Е. Алоян. — Новосибирск: Наука, 1985. -254 с.

Бартеньев, О. В. Современный фортран — М.: Диалог-МИФИ, 2000. — 390 с.

Приложение

Листинг программы

program DiffTrans

implicit none

! Variables

integer n

real (8) L, T, M, constX, constZ, u, w, mu, nu, lambda, tau, tc

namelist /input/ L, T, M, constX, constZ, n, u, w, mu, nu, lambda, tau

integer i, j, k, z

real (8) fi (0:100,0:100), fipr (0:100,0:100), fi1 (0:100,0:100)

real (8) a1 (0:100), a2 (0:100), a3 (0:100), b (0:100), pr (0:100)

real (8) h, o, s, f

external f

! Body of DiffTrans

open (1, FILE = 'in.txt')

open (2, FILE = 'out.txt')

read (1, input)

h=L/n! 0.1

do i=1, n-1

do j=1, n-1

fi (i, j)=0.0

end do

end do

do j=1, n-1

fi (0, j)=0.0

fi (n, j)=0.0

end do

do i=1, n-1

fi (i, n)=0.0

fi (i, 0)=fi (i, 1)*nu/(h*lambda+nu)

end do

do k=0, n

write (2,*) k

write (2,' ((1x, f8.5))') (fi (i, k), i=0, n)

end do

tc=0.0

fi1=fi

z=0.0

do while ((abs (tc-T)>=tau).and. (z<=1 000 000))

z=z+1

write (2,*), 'tc=', tc

tc=tc+tau

fi=fi1

do j=1, n-1

a1 (0)=0.0

a2 (0)=1.0

a3 (0)=-1.0

b (0)=0.0

a1 (n)=1.0

a2 (n)=-1.0

a3 (n)=0.0

b (n)=0.0

do k=1, n-1

a1 (k)=(((-u*tau)/(4*h)) — (mu*tau)/(2*h*h))

a2 (k)=(1.0+(mu*tau)/(h*h))

a3 (k)=(((u*tau)/(4*h)) — (mu*tau)/(2*h*h))

b (k)=(fi (k, j) — ((tau*u*(fi (k+1, j) — fi (k-1, j)))/(4*h)+(tau*mu*(fi (k+1, j) — 2*fi (k, j)+fi (k-1, j)))/(2*h*h)))

end do

call progonka (a1, a2, a3, b, pr, n);

do k=0, n

fipr (k, j)=pr (k)

end do

end do

do k=1, n-1

fipr (k, n)=0.0

fipr (k, 0)= fipr (i, 1)*nu/(h*lambda+nu)

end do

fi=fipr

do i=1, n-1

a1 (0)=0.0

a2 (0)= - nu/h-lambda

a3 (0)=nu/h

b (0)=0.0

a1 (n)=0.0

a2 (n)=1.0

a3 (n)=0.0

b (n)=0.0

do k=1, n-1

a1 (k)=(((-w*tau)/(4*h)) — (nu*tau)/(2*h*h))

a2 (k)=(1.0+(nu*tau)/(h*h))

a3 (k)=(((w*tau)/(4*h)) — (nu*tau)/(2*h*h))

b (k)=(fi (i, k) — ((tau*w*(fi (i, k+1) — fi (i, k-1)))/(4*h)+(tau*nu*(fi (i, k+1) — 2*fi (i, k)+fi (i, k-1))/(2*h*h)))+tau*f (i*h, k*h, M, constX, constZ))

end do

call progonka (a1, a2, a3, b, pr, n)

do k=0, n

fipr (i, k)=pr (k)

end do

end do

fi=fipr

fi1=fipr

do k=0, n

write (2,*) k

write (2,' ((1x, f8.5))') (abs (fi (i, k)), i=0, n)

end do

end do

close (1)

close (2)

end program DiffTrans

function f (o, s, consf, cx, cy)

real (8) f, consf, cx, cy, buf, o, s

buf=0

if (o==cx.and. s==cy) then

buf=consf*consf

else

buf=0

endif

f=buf

end function f

subroutine progonka (a1, a2, a3, b, pr, n)

real (8) a1 (0:100), a2 (0:100), a3 (0:100), b (0:100), pr (0:100)

real (8) pr1 (0:100), pr2 (0:100)

integer i1, n

intent (out) pr

pr1 (n-1)=-a1 (n)/a2 (n)

pr2 (n-1)=b (n)/a2 (n)

do i1=n — 1,1, — 1

pr1 (i1−1)=-a1 (i1)/(pr1 (i1)*a3 (i1)+a2 (i1))

pr2 (i1−1)=(b (i1) — pr2 (i1)*a3 (i1))/(pr1 (i1)*a3 (i1)+a2 (i1))

end do

pr (0)=(b (0) — a3 (0)*pr2 (0))/(a2 (0)+a3 (0)*pr1 (0))

do i1=1, n

pr (i1)=pr1 (i1−1)*pr (i1−1)+pr2 (i1−1)

end do

end subroutine