Цели и задачи математического моделирования процессов и систем

1. ЦЕЛИ И ЗАДАЧИ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ПРОЦЕССОВ И СИСТЕМ 1.1. ПОНЯТИЕ «МАТЕМАТИЧЕСКАЯ МОДЕЛЬ» Математическое моделирование позволяет до создания реальной системы (объекта) или возникновения реальной ситуации рассмотреть возможные режимы работы, выбрать оптимальные управляющие воздействия, составить объективный прогноз будущих состояний системы. Вычислительные эксперименты, проводимые на основе математических моделей, помогают увидеть за частным общее, развить универсальные методы анализа объектов различной физической природы, познать свойства изучаемых процессов и систем. Наконец, математическое моделирование является основой интенсивно разрабатываемых автоматизированных систем проектирования, управления и обработки данных. Основная задача математического моделирования – выделение законов в природе, обществе и технике и запись их на языке математики Например: 1) Зависимость между массой тела m, действующей на него силой F и ускорением его движения а записывается в форме 2-го закона Ньютона: F = m⋅ a; 2) Зависимость между напряжением в электрической цепи U, ее сопротивлением R и силой тока I записывается в виде закона Ома: I = U/R. Существует множество определений математической модели. Приведем одно из них: Математической моделью некоторого объекта, процесса или явления будем называть запись его свойств на формальном языке с целью получения нового знания (свойств) об изучаемом процессе путем применения формальных методов Альтернативой формальному (математическому) подходу является экспериментальный подход. К его недостаткам можно отнести: 1) высокая стоимость подготовки и проведения экспериментов; 2) получение частного знания (знания о конкретном объекте исследования, а не о классе объектов). Например, пусть требуется определить воздействие х на некоторый процесс или объект, при котором его результирующая характеристика у имеет максимально возможное значение (Рис. 1.1). у На ча ло уmax x1 x2 x3 х4 x5 x6 x7 x8 у = f(x) х хmax а x б Рис. 1.1 На рис. 1.1. а) показан эмпирический (экспериментальный) подход к решению поставленной задачи, который состоит в экспериментальном определении значения параметра у для нескольких значений входного воздействия х. Среди них найдено наибольшее, и оно принимается за максимум. Как видим из этого рисунка, возможно несколько значений воздействия х (х4 и х5), при которых у имеет наибольшее значение, но ни одно из них не является настоящим максимумом, который, возможно, лежит между ними. Математический подход (рис. 1.1. б) предполагает наличие математической модели процесса типа y = f(x). Взяв производную df и приравняв ее к нулю, получим dx уравнение, решением которого является точное значение xmax , доставляющее максимум функции у. Схема применения математической модели при решении реальных задач имеет вид, показанный на рис. 1.2. Реальный объект Постановка задач Разработка математическо Получение решения Выбор (доработка) методов Рис. 1.2 Модель сложного объекта (процесса, системы) не может быть простой. Из чего следует, что процесс использования математических моделей реальных систем является итерационным процессом, когда последовательно уточняется (дорабатывается) математическая модель и методы решения стоящих задач Важнейшей характеристикой моделей является их точность, адекватность действительности. При этом важно иметь в виду, что все модели представляют собой приближенное описание реальных объектов (процессов) и поэтому принципиально неточны. Интегральная оценка модели может быть получена путем сравнения результатов моделирования и экспериментальных данных для конкретных объектов или режимов. Для оценки значимости совпадения или несовпадения модельных и экспериментальных результатов широко используются методы математической статистики. Вместе с тем не следует переоценивать результаты такой проверки. Хорошее совпадение модельных и экспериментальных данных, вообще говоря, не доказывает точности модели, а лишь подтверждают ее функциональную пригодность для моделирования. Всегда может быть предложена модель, обеспечивающая лучшее совпадение с экспериментом, но не лучшее описание моделируемого объекта или процесса. 1.2. КЛАССИФИКАЦИЯ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ Существует несколько схем классификации математических моделей. Все они достаточно условны. Одна из таких схем приведена на рис. 1.3. Все математические модели по использованному формальному языку можно разбить на аналитические и имитационные. Аналитические – модели, в которых используется стандартный математический язык. Имитационные – модели, в которых использован специальный язык моделирования или универсальный язык программирования. Аналитические модели могут быть записаны в виде формул или уравнений. Если какой-либо процесс не может быть описан в виде аналитической модели, его описывают с помощью специального алгоритма или программы. Такая модель является имитационной. Аналитические модели в свою очередь разбиваются на теоретические и эмпирические модели. Теоретические модели отражают реальные структуры и процессы в исследуемых объектах, то есть, опираются на теорию их работы. Эмпирические модели строятся на основе изучения реакций объекта на изменение условий окружающей среды. При этом теория работы объекта не рассматривается, сам объект представляет собой так называемый «черный ящик», а модель – некоторую интерполяционную зависимость. Эмпирические модели могут быть построены на основе экспериментальных данных. Эти данные получают непосредственно на исследуемых объектах или с помощью их физических моделей. Математические модели Аналитические Имитационные Теоретические Эмпирические Теоретические Линейные Нелинейные Нелинейные Статические Динамические Динамические Детерминированные Стохастические Детерминированные Аналитически разрешимые Численно разрешимые Численно разрешимые Рис. 1.3 По форме описания аналитические модели подразделяются на линейные и нелинейные. Если все входящие в модель величины не зависят от времени, то имеем статическую модель объекта или процесса, в противном случае получаем динамическую модель. В детерминированных моделях все взаимосвязи, переменные и константы заданы точно, что приводит к однозначному определению результирующей функции. Если часть или все параметры, входящие в модель по своей природе являются случайными величинами или случайными функциями, то модель относят к классу стохастических моделей. В стохастических моделях задаются законы распределения случайных величин, что приводит к вероятностной оценке результирующей функции. Если аналитическое исследование может быть доведено до конца, модели называются аналитически разрешимыми. В противном случае говорят о численно разрешимых аналитических моделях. Математическое моделирование занимает особое место в различных отраслях науки и техники. Под математическим моделированием понимают процесс установления соответствия данному реальному объекту некоторого математического объекта, называемого математической моделью, и исследование этой модели, позволяющее получать характеристики рассматриваемого реального объекта. Вид математической модели зависит как от природы реального объекта, так и от задач исследования объекта и требуемой достоверности и точности решения этой задачи. Любая математическая модель, как и всякая другая, описывает реальный объект лишь с некоторой степенью приближения к действительности. Математическое моделирование для исследования характеристик процесса функционирования систем можно разделить на аналитическое, имитационное и комбинированное. Для аналитического моделирования характерно то, что процессы функционирования функциональных элементов соотношений системы записываются (алгебраических, в виде некоторых интегродифференциальных, конечно-разностных и т.п.) или логических условий. Аналитическая модель может быть исследована следующими методами: а) аналитическим, когда стремятся получить в общем виде явные зависимости для искомых характеристик; б) численным, когда, не умея решать уравнений в общем виде, стремятся получить числовые результаты при конкретных начальных данных; в) качественным, когда, не имея решения в явном виде, можно найти некоторые свойства решения (например, оценить устойчивость решения). Наиболее полное исследование процесса функционирования системы можно провести, если известны явные зависимости, связывающие искомые характеристики с начальными условиями, параметрами и переменными системы S . Однако такие зависимости удается получить только для сравнительно простых систем. При имитационном воспроизводит процесс моделировании функционирования реализующий системы S во модель времени, алгоритм причем имитируются элементарные явления, составляющие процесс, с сохранением их логической структуры и последовательности протекания во времени, что позволяет по исходным данным получить сведения о состояниях процесса в определённые моменты времени, дающие возможность оценить характеристики системы S. Основным преимуществом имитационного моделирования по сравнению с аналитическим является возможность решения более сложных задач. Имитационные модели позволяют достаточно просто учитывать такие факторы, как наличие дискретных и непрерывных элементов, нелинейные характеристики элементов системы, многочисленные случайные воздействия и др., которые часто создают трудности при аналитических исследованиях. В настоящее время имитационное моделирование – наиболее эффективный метод исследования больших систем, а часто и единственный практически доступный метод получения информации о поведении системы, особенно на этапе ее проектирования. Метод имитационного моделирования позволяет решать задачи анализа больших систем S , включая задачи оценки: вариантов структуры системы, эффективности различных алгоритмов управления системой, влияния изменения различных параметров системы. Имитационное моделирование может быть также положено в основу структурного, алгоритмического, параметрического синтеза больших систем, когда требуется создать систему с заданными характеристиками при определенных ограничениях, которая является оптимальной по некоторым критериям оценки эффективности. Метод статистического моделирования (метод статистических испытаний, или метод Монте-Карло) – многократное воспроизведение процесса путём "прогонов" имитационной модели на ЭВМ с последующей статистической обработкой информации для нахождения характеристик исследуемого процесса. Метод особенно эффективен, когда параметры модели и полученные результаты моделирования являются случайными величинами или реализациями случайных функций. В последнее время метод применяется для машинной имитации технологических систем, подверженных случайным воздействиям, с целью оптимизации обработки по производительности и другим параметрам при обязательном условии обеспечения регламентируемых параметров качества обработки с заданной надёжностью. Комбинированное (аналитико-имитационное) моделирование при анализе и синтезе систем позволяет объединить достоинства аналитического и имитационного моделирования. При предварительная построении декомпозиция комбинированных процесса моделей функционирования проводится объекта на составляющие подпроцессы и для тех из них, где это возможно, используются аналитические модели, а для остальных подпроцессов строятся имитационные модели. Такой комбинированный подход позволяет охватить качественно новые классы систем, которые не могут быть исследованы с использованием только аналитического и имитационного моделирования в отдельности. Кибернетическое моделирование характеризуется тем, что в нём отсутствует непосредственное подобие физических процессов, происходящих в моделях, реальным процессам. В этом случае стремятся отобразить лишь некоторую функцию и рассматривают реальный объект как "чёрный ящик", имеющий ряд входов и выходов, и моделируют некоторые связи между ними. Чаще всего при использовании кибернетических моделей проводят анализ поведенческой стороны объекта при различных воздействиях внешней среды. Таким образом, в основе кибернетических моделей лежит отражение некоторых информационных процессов управления, что позволяет оценить поведение реального объекта. Для построения имитационной модели в этом случае необходимо выделить исследуемую функцию реального объекта, попытаться формализовать эту функцию в виде некоторых операторов связи между входом и выходом и воспроизвести на имитационной модели данную функцию, причём на базе совершенно иных математических соотношений и, естественно, иной физической реализации процесса. 2. ТЕОРЕТИЧЕСКИЕ МАТЕМАТИЧЕСКИЕ МОДЕЛИ АНАЛИТИЧЕСКОГО ТИПА Простейшие аналитические модели могут быть заданы явно в виде функции одной или нескольких переменных. Обычно в виде функций задаются общие законы природы или общие закономерности, полученные в результате интегрирования дифференциальных уравнений. Примером такой модели может служить знаменитая формула К.Э. Циолковского: ∆vла = v ln M0 Mк , определяющая приращение скорости ракеты ∆v ла при импульсном сжигании топлива через скорость истечения рабочего тела v и отношение начальной М0 и конечной Mк масс ракеты. Модель, заданная в явном виде, дает исчерпывающее описание исследуемого объекта. Она позволяет построить зависимость его характеристик от управляющих факторов, взять производные и найти экстремумы модели, определить характеристики модели в окрестности экстремумов и т.д. Очень удобна графическая интерпретация таких моделей. Однако модели в виде формул могут быть разработаны только для очень простых объектов. 2.1. ПОСТРОЕНИЕ МАТЕМАТИЧЕСКОЙ МОДЕЛИ СВЕРЛЕНИЯ ЛАЗЕРОМ Примером аналитической теоретической модели может служить модель, описывающая глубину отверстия при лазерном сверлении. Резание и сверление металлов весьма важно для многих областей техники. Значительный интерес представляет создание новых устройств, предназначенных для специальных материалов, а также для тех случаев, когда желательно обеспечить некоторую степень автоматизации указанных процессов. В последнее время для этого были предприняты попытки использования мощных лазеров. Основная идея состоит в том, чтобы сфокусировать значительную мощность на малой площади поверхности материала, создавая таким образом интенсивный нагрев и испарение с последующим образованием отверстия. При сверлении необходимо постараться обеспечить такие условия процесса, чтобы проделанное отверстие прямо проходило сквозь материал, и избежать, таким образом, затекания расплавленного металла обратно в отверстие и застывания его там. Построим математическую модель, главная применимость которой – глубокое сверление. При помощи модели попытаемся ответить на вопрос, как быстро можно проделать отверстие, используя пучок излучения высокой мощности, и на какую глубину. Рассмотрим высокоэнергетический пучок лазерного W излучения, сфокусированный на малом участке поверхности металла S (Рис. 2.1). Определенная доля энергии поглощается, а остальная часть отражается. Поглощение энергии происходит внутри слоя, толщина которого много меньше миллиметра, вызывает поверхностный нагрев материала и рост температуры поверхности. Температура растет не безгранично. Существует два процесса, A Рис. 2.1 ограничивающие рост температуры: − перенос тепла в глубь материала от нагретых к холодным участкам, обусловленный теплопроводностью; − испарение. Когда температура материала достигает точки кипения, скрытое тепло поглощается без дальнейшего увеличения температуры в процессе испарения материала. При удалении пара от поверхности материала в металле образуется выемка. Задача количественного описания этого процесса и вызывает необходимость математического моделирования. Будем рассматривать модель, описывающую процесс разрушения материала, при котором вся энергия лазерного излучения используется только для испарения материала. Этот предельный режим испарения может возникать двумя путями: − когда энергия поступает на поверхность слишком быстро, так что тепло не успевает распространиться в глубь металла; − плотность мощности пучка постоянна, а распределение температуры впереди границы области испарения приближается к стационарному. Предположим, что мощность W распределена по некоторой площади А поверхности; излучение приложено по нормали к поверхности (см. рис. 2.1). За интервал времени ∆t поступает энергия W⋅∆t. Пусть глубина возникающей выемки равна ∆S, тогда объем испарившегося материала равен A⋅∆S. Используя закон сохранения энергии, получим h⋅ρ⋅A⋅∆S = W⋅∆t, где h – количество тепла, требуемое для испарения единицы массы материала; ρ – плотность материала. Преобразуем это выражение и положим ∆t → 0, получим скорость роста глубины выемки: dS W 1 = ⋅ . dt A hρ Это уравнение показывает, что для любого материала предельная скорость пропорциональна плотности энерговыделения W/A. Интегрируя это уравнение и полагая S = 0 при t = 0, найдем глубину выемки в произвольный момент времени t: S (t ) = 1 t ∫ Wdt hρ A 0 (1) или S (t ) = E (t ) , hρ A где E(t) – полная энергия, выделенная источником за промежуток времени (0, t). Таким образом, в предельном режиме испарения глубина выемки зависит только от полной энергии, поступившей на поверхность. Формула (2.1) представляет собой теоретическую аналитически-разрешимую динамическую детерминированную модель. На практике всегда существует перенос некоторого количества тепла в материал за счет теплопроводности. Общая задача движения границы раздела фаз с учетом теплопроводности известна как задача Стефана. 2.2. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НА ОСНОВЕ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С ЧАСТНЫМИ ПРОИЗВОДНЫМИ Ряд задач, связанных с использованием физических полей, приводит к моделям в виде дифференциальных уравнений в частных производных. Особенностью таких задач является то, что изучаемые параметры изменяются не только во времени, но и зависят от координат x, y, z рассматриваемого пространства. Такие модели называются нестационарными. Модели, в которых параметры не зависят от времени, называются стационарными. К таким моделям сводятся описания полей температур в элементах конструкции двигателя и полей скоростей при течении жидкости (газа). Уравнениями в частных производных описываются колебания элементов конструкции и поля напряжений, возникающих при работе этих элементов. Линейное дифференциальное уравнение в частных производных имеет вид a0 ∂Φ (t ) ∂Φ(t ) ∂Φ (t ) ∂Φ(t ) + a1 + a2 + ... + ak = f ( x1 , x2 ,..., xk , t ) . ∂t ∂x1 ∂x2 ∂xk Математическая модель, описанная дифференциальными уравнениями в частных производных, должна включать в себя необходимые для решения задачи краевые условия: 1) Должна быть задана область D, ограниченная поверхностью (на плоскости – кривой) Γ, в которой определяется решение. 2) Должны быть заданы условия на границе Γ этой области. В случае нестационарного поля эти граничные условия, так же как и сама область могут меняться во времени. Граничные условия могут быть 1-го, 2-го, 3-го и 4-го рода: а) Граничные условия 1-го рода предусматривают задание на границе величины искомой функции: Φ |Γ = f1 (Γ) – для стационарного поля; Φ(t ) |Γ = f1 (Γ, t ) – для нестационарного поля. б) Граничные условия 2-го рода – предусматривают задание производной искомой функции: ∂Φ = f 2 (Γ) – для стационарного поля; ∂xi Γ ∂Φ (t ) = f 2 (Γ, t ) – для нестационарного поля. ∂xi Γ в) Граничные условия 3-го рода – предусматривают комбинации функции и ее производной:  ∂Φ   aΦ + b  = f 3 (Γ) – для стационарного поля; ∂xi   Γ  ∂Φ   a(t )Φ (t ) + b(t ) (t )  = f 3 (Γ, t ) – для нестационарного поля. ∂xi   Γ г) Граничные условия 4-го рода (условия сопряжения). Ставятся на границе контакта двух сред с различными характеристиками. Для нестационарных полей должны быть заданы одно или два начальных условия, характеризующих состояние поля в начальный момент времени: Φ( xi ) t = 0 = f 4 ( xi ); (i = 1, 2, 3). ∂Φ ( xi ) = f 5 ( xi ), ∂t t =0 Здесь xi – координаты пространства. Совокупность уравнений и краевых (и начальных) условий полностью определяет модель и позволяет провести ее исследование. Решение часто задается в x2 Γ Φ = const x1 Рис. 2.2 виде семейств изолиний Φ = const (Рис. 2.2). 2.3. ПОСТАНОВКА ЗАДАЧ ДЛЯ НЕСТАЦИОНАРНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ В теории теплопроводности для отыскания температурного поля в твердом теле, а также в неподвижной жидкости или газообразной среде используется уравнение теплопроводности ρcp ∂T = div ( λ grad T ) + qv , ∂τ где T ( x, y, z, τ ) − температура тела в точке ( x, (1) y, z ) в момент времени τ , К ; ρ − плотность вещества, кг м3 ; c p − удельная теплоемкость вещества, Дж ( кг ⋅ К ) ; λ – коэффициент внутренней теплопроводности, Вт ( м ⋅ К ) ; qv − мощность источников теплоты, Вт м3 . Рассматривая постоянными ρ, cp , λ температуропроводности a = и вводя коэффициент λ , м 2 с , уравнение (1) принимает вид: ρcp q ∂T = a∇ 2T + v . ρcp ∂τ (2) При этом скалярное произведение оператора Гамильтона самого на себя дает оператор Лапласа: div ( grad T ) = ∇ 2T = ∆T . В декартовых координатах ∆T = координатах ∂ 2T 1 ∂T 1 ∂ 2T ∂ 2T ∆T = 2 + + + ∂r r ∂r r 2 ∂ϕ 2 ∂z 2 сферических координатах ∆T = ∆T = (3) ∂ 2T ∂ 2T ∂ 2T + + , ∂x 2 ∂y 2 ∂z 2 или ∆T = в цилиндрических 1 1 ( rTr′)′r + 2 Tϕϕ′′ + Tzz′′ , r r ∂ 2T 2 ∂T 1 ∂ 2T ctgθ ∂T 1 ∂ 2T + + + + r 2 ∂θ r 2 sin 2 θ ∂ϕ 2 ∂r 2 r ∂r r 2 ∂θ 2 в или 1 2 1 1 r ⋅ Tr′)′ + 2 ( sin θ ⋅ Tθ′ )′θ + 2 2 Tϕϕ′′ . 2 ( r r r sin θ r sin θ Без наличия различного рода источников или поглотителей тепла ( qv = 0 ) уравнение теплопроводности является однородным и имеет вид: ∂T = a∆T . ∂τ (4) При решении конкретных задач теплопроводности физические свойства тела ( ρ , c p , λ ) считаются известными, определены геометрическая форма и размеры тела или области пространства, в которой протекает процесс теплопроводности, заданы начальные и граничные условия. Начальным условием задается распределение температуры внутри тела и на его границах в момент времени τ = 0 : T ( x, y, z , τ ) τ =0 = f ( x, y, z ), (5) где f ( x, y, z ) − заданная функция. Граничным условием первого рода задается температура поверхности тела (при τ > 0 ) как функция координат точек его границ: (6) T ( x, y, z , τ ) Г = g ( x, y, z ,τ ), т. е. это условие определяет неизвестную функцию температуры T через известные функции. Например, T x =0 = Tc = const. Граничным условием второго рода задается плотность теплового потока qc , Вт м 2 , на границах тела: −λ ∂T = qc , ∂n Г (7) где λ − теплопроводность в данной точке границы тела, n − внешняя нормаль к элементу поверхности. При этом qc > 0 , если тепловой поток направлен от тела в окружающую среду, и qc < 0 , если наоборот. Например, на границе x = 0 получаем −λ ∂T ∂x = qc . x =0 Граничное условие третьего рода записывается в виде: −λ ∂T = α (T Г − Tж ), ∂n Г (8) где Tж − температура среды, α − коэффициент теплоотдачи, Вт ( м 2 ⋅ К ) , который характеризует интенсивность процесса конвективного теплообмена между поверхностью тела и окружающей средой. Граничным условием третьего рода задается зависимость теплового потока qc от температуры поверхности тела Tс и температуры среды Tж , окружающей тело, в соответствии с законом НьютонаРихмана: qc = α (Tc − Tж ). (9) В случае, когда помимо конвективного теплообмена происходит теплообмен излучением между данным и окружающими телами, в правую часть соотношения (1.8) надо включить qизл , Вт м 2 , − величину, равную разности плотностей потоков собственного и поглощенного излучений. Если поверхность окружающих тел 4 представить в виде оболочки с температурой Tокр , то qизл = ε σ 0 (Tс4 − Tокр ), где ε − степень черноты данного тела ( 0 ≤ ε ≤ 1 ), σ 0 = 5, 67 ⋅10−8 Вт ( м 2 ⋅ K 4 ) − постоянная Стефана–Больцмана. Таким образом, с учетом теплообмена излучением граничное условие третьего рода можно записать в виде: −λ ∂T ∂n Г 4 = α (T Г − Tокр ) + ε σ 0 (T 4 − Tокр ). Г (10) Уравнение теплопроводности является дифференциальным уравнением с частными производными. Его решение в конкретных случаях находиться с привлечением начальных и граничных условий. Аналитические методы решения задач теплопроводности достаточно сложны и подробно рассматриваются в специальных курсах математики, включающих тематику уравнений математической физики. Универсальность численных методов (конечных разностей, конечных элементов, контрольного объема) позволяет получать решение практически для любой задачи теплопроводности. Граничные условия четвертого рода задаются на границе двух сред, обладающих различными теплофизическими характеристиками. Условия идеального контакта имеют вид:  T1 г = T2 г ;   ∂T1 ∂T = λ2 2 , λ1 ∂n г  ∂n г (11) где T1 , T2 − температуры соприкасающихся сред. В случае неидеального контакта температура при переходе границы терпит разрыв, пропорциональный тепловому потоку и условия сопряжения имеют вид:  −1 ∂T  R1,2 (T2 г − T1 г ) = λ1 ∂n ;  г  λ ∂T1 = λ ∂T2 , 2  1 ∂n г ∂n г где R1,2 − тепловое сопротивление контакта. 2.4. НЕКОТОРЫЕ РАЗНОСТНЫЕ СХЕМЫ ДЛЯ НЕСТАЦИОНАРНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ (12) Рассмотрим применение метода сеток с использованием конечно-разностных аппроксимаций (метод конечных разностей) на примере однородного нестационарного уравнения теплопроводности ∂T ∂ 2T = a 2 , x ∈ ( 0, L ) , τ ∈ ( 0, τ m ) ; ∂τ ∂x (1) T ( x, 0 ) = γ = const при x ∈ [ 0, L ] ; (2) T ( 0, τ ) = α (τ ) , T ( L, τ ) = β (τ ) при τ ∈ [ 0, τ m ] . (3) Решение уравнения (1) зависит от двух переменных: τ − времени и x – пространства. Система координат выбирается так, чтобы в ней переменная x менялась вдоль оси абсцисс, а переменная τ − вдоль оси ординат (рис. 2.4). τ τm τj xi L x Рис. 2.4 Для решения уравнения (1) конечно-разностным методом построится сетка, покрывающая прямоугольник [0, L] ∪ [0,τ m ] . Координаты узлов сетки, образованные пересечением вертикальных и горизонтальных отрезков, определяются по формулам xi = i ⋅ h , τ j = j ⋅ k , где h − шаг по пространству, h = k= τm M L , 0 ≤ i ≤ N ; k − шаг по времени, N , 0≤ j≤M . Далее на построенной сетке вводят сеточные функции – это функции определенные в узлах сетки. T ( xi , τ j ) = Ti j , α (τ j ) = α j , β (τ j ) = β j . Для аппроксимации ∂T ∂ 2T и 2 можно использовать следующие формулы: ∂τ ∂x (4) ∂T Ti j +1 − Ti j = + O (k ), ∂τ k (5) ∂T Ti j − Ti j −1 = + O (k ), ∂τ k (6) ∂T Ti j +1 − Ti j −1 = + O (k 2 ), ∂τ 2k (7) ∂ 2T Ti +j1 − 2Ti j + Ti −j1 = + O ( h2 ) , 2 2 ∂x h (2.8) где O ( k ) , O ( k 2 ) и O ( h 2 ) соответствующие погрешности аппроксимации (невязки). Приведенные формулы разностных аналогов частных производных по времени (правая, левая и центральная разностные производные) и пространству (центральная разностная производная) получают, используя разложение функции двух переменных T ( x, τ ) в ряд Тейлора (ограничиваясь двумя первыми членами для производных первого порядка и тремя членами для производных второго порядка). Подставляя приведенные разностные аналоги производных по времени и пространству в соотношение (1), получим сеточные (разностные) уравнения: Ti j +1 − Ti j T j − 2Ti j + Ti −j1 = a i +1 , 1 ≤ i ≤ N − 1, 0 ≤ j ≤ M − 1; k h2 (9) Ti j − Ti j −1 T j − 2Ti j + Ti −j1 = a i +1 , 1 ≤ i ≤ N − 1, 1 ≤ j ≤ M ; k h2 (10) Ti j +1 − Ti j −1 Ti +j1 − 2Ti j + Ti −j1 =a , 1 ≤ i ≤ N − 1, 1 ≤ j ≤ M − 1; 2k h2 (11) Аппроксимация начальных (2) и граничных (3) условий имеет вид: Ti 0 = γ , 0 ≤ i ≤ N ; (12) T0 j = α j , Tn j = β j , 0 ≤ j ≤ M . (13) Определение 1. Сеточные уравнения в совокупности с аппроксимациями начальных и граничных условий для определения значений сеточной функции в области исследования называются разностными схемами. Конфигурация узлов расчетной области, используемых на каждом элементарном шаге вычисления, называют шаблоном разностной схемы. Шаблон разностной схемы в соответствии с сеточным уравнением (9) приведен на рисунке 2.5. i, j+1 i-1, j i, j i+1, j Рис. 2.5 На рисунке 2.6 представлен шаблон разностной схемы в соответствии с сеточным уравнением (10). i-1, j i, j i+1, j i, j-1 Рис. 2.6 Шаблон разностной схемы для сеточного уравнения (11) представлен на рисунке 2.7. i, j+1 i-1, j i, j i+1, j i, j-1 Рис. 2.7 Определение 2. Если в разностной схеме присутствуют значения искомой функции на двух соседних слоях, то схема называется двухслойной (например, рис. 2, 5), соответственно, на трех слоях – трехслойной (например, рис. 2.7). Определение 3. Явными схемами называются разностные схемы, для которых искомые значения на следующем слое по времени находятся непосредственно из данных на предыдущем слое. Например, формулы (9) и (11), шаблоны на рисунках 2.5 и 2.7 соответственно. Если для определения значений сеточной функции при переходе от слоя к слою необходимо решать систему линейных алгебраических уравнений, то схема называется неявной. Например, формула (10), шаблон на рисунке 2.6. 2.5. УСТОЙЧИВОСТЬ И СХОДИМОСТЬ РАЗНОСТНЫХ СХЕМ ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ Замена частных производных соответствующими конечными разностями приводит к появлению погрешности аппроксимации (невязки), которая представляет собой разность между точным решением и приближенным, т.е. ∂T ∂τ i, j Ti j +1 − Ti j = + O(k ). k (1) Определение 4. Разностная схема называется устойчивой, если малые изменения правой части дифференциального уравнения приводят к незначительным отклонениям решения. Рассмотрим анализ устойчивости на примере явной разностной схемы для уравнения теплопроводности. Обозначая сеточную функцию на верхнем временном слое Tˆi , запишем уравнение (9, п. 2.4) в следующем виде: Tˆi − Ti T − 2Ti + Ti −1 . = a i +1 k h2 (2) Выражая Tˆi , получим формулу для вычисления значений сеточной функции: ak Tˆi = Ti + 2 (Ti +1 − 2Ti + Ti −1 ) . h (3) Выясним, как изменяется погрешность решения при переходе от слоя к слою. Предположим, что на нижнем временном слое сеточная функция известна с некоторой погрешностью относительно точного решения T ( x ) : Ti = T ( x ) + ε i . (4) Используя формулу (4) для подстановки в разностное уравнение (3.3) с учетом того, что точные T ( x ) удовлетворяют исходному уравнению теплопроводности, получаем разностное уравнение для динамики погрешности ε i : εˆi = ε i + ak (ε i +1 − 2ε i + ε i −1 ) . h2 (5) Разделив выражение (5) на ε i , получим относительное изменение погрешности при переходе на новый слой по времени: εˆi ε +ε  ak  = 1 − 2  2 − i +1 i −1  . εi εi h   (6) Если правая часть равенства (6) меньше 1, то схема устойчива (погрешности не нарастают), иначе схема неустойчива. Предположим, что погрешность ε i имеет гармоническую зависимость (настоящая погрешность может быть выражена через спектр гармоник с различными ε и α ): ε i = ε exp (α lxi ) = ε exp (α lih ) , (7) где l 2 = −1 . Подставляя выражение (3.7) в оценку (3.6), получим: εˆi 2a k  exp (α lh ) + exp ( −α lh )  = 1 − 2 1 − . εi h  2  (8) Так как в формуле (3.8) каждая из экспонент от мнимого аргумента может принимать значения от -1 до 1, то выражение в скобках изменяется в диапазоне от 0 до 2. Следовательно, абсолютное значение относительной погрешности превысить 1, если εˆi может εi 2a k >1. h2 Таким образом, явная схема (3.2) для уравнения теплопроводности является условно устойчивой и при решении необходимо учитывать следующее соотношение сеточных параметров: 2a k ≤1. h2 Проводя аналогичные рассуждения для неявной схемы, получим оценку относительной погрешности всегда меньшую 1. Поэтому неявная схема является безусловно устойчивой, т.е. позволяет проводить расчеты при любых значениях h и k. Определение 5. Говорят, что решение разностной задачи y ( h ) сходится к решению дифференциальной задачи при h → 0, если проекция точного решения на разностную сетку. y( ) − Y ( h h) → 0, где Y ( h ) − В случае выполнения условия y( ) − Y ( h h) ≤ c ⋅ h p , c ≠ c ( h ) , сходимость имеет прядок p. Чтобы исследовать схему на сходимость, необходимо знать точное решение дифференциальной задачи. Если получение точного решения затруднительно, то разностные схемы исследуют на сходимость, используя теорему П. Лакса и В.С. Рябенького. Теорема (П. Лакса - В.С.Рябенького). Решение линейной разностной задачи сходиться к решению дифференциальной, если разностная задача устойчива и аппроксимирует дифференциальную задачу. При этом порядок аппроксимации совпадает с порядком сходимости. Таким образом, явная и неявная аппроксимации O ( k + h2 ) , двухслойные схемы имеют погрешность обладают свойством устойчивости, соответственно, решение разностной задачи сходится к решению дифференциальной задачи. 3. ТЕМПЕРАТУРНОЕ ПОЛЕ В ПРОЦЕССЕ ОХЛАЖДЕНИЯ ИЛИ НАГРЕВАНИЯ НЕОГРАНИЧЕННОЙ ПЛАСТИНЫ Рассмотрим задачу о симметричном нестационарном охлаждении пластины в среде с постоянной температурой и при постоянном коэффициенте теплоотдачи. Дана пластина толщиной 2 R. Если толщина пластины мала по сравнению с длиной и шириной, то такую пластину обычно считают неограниченной. В начале процесса охлаждения температура по всему объему пластины была одинаковой и равной tн . В течение всего процесса на её поверхности поддерживается постоянная температура tc . Требуется определить температуру t ( x, τ ) в любой точке в любой момент времени τ . В данной задаче используем обозначение: t ( x, τ ) − температура, °C. При заданных условиях охлаждения задача становится симметричной и начало координат удобно поместить на оси пластины. Таким образом, нужно найти функцию t ( x, τ ) , удовлетворяющую следующим условиям: ∂t ∂ 2t = a 2 ; x ∈ ( − R, R ) , τ ∈ ( 0, τ m ) ; уравнению теплопроводности: ∂τ ∂x (1) начальному условию: t ( x, 0 ) = tн , x ∈ [ − R, R ] ; (2) граничным условиям: t ( R, τ ) = t ( − R, τ ) = tc , τ ∈ [0, τ m ]. (3) Условие t ( − R, τ ) = tc можно заменить следующим отсутствие теплового потока в плоскости ∂t ∂x = 0, что означает x =0 x = 0 и следует из симметрии температурного поля пластины при заданных условиях. Используя двухслойную явную схему (2.9), имеем разностное уравнение: tij +1 − tij t j − 2ti j + tij−1 = a i +1 , 1 ≤ i ≤ N − 1, 0 ≤ j ≤ M − 1. k h2 (4) Выражая из соотношения (4.4) ti j +1 , получим формулу для вычисления значений сеточной функции: tij +1 = C ( tij−1 + tij+1 ) + (1 − 2C ) ti j , где C = (5) ak − число Куранта. h2 При реализации граничных условий первого рода значение функции на границе находится непосредственно из граничного условия: t Nj +1 = tc . (6) Для граничных условий второго рода производную в точке x0 ( x = 0 ) заменяют конечно-разностным аналогом со вторым порядком точности: t1j − t−j1 = 0. 2h (7) Следовательно, значение температуры, записанное в фиктивном узле i = −1, равно: t−j1 = t1j . (8) Из предположения, что решение непрерывно записываем уравнение (5) для i =0: t0j +1 = C ( t−j1 + t1j ) + (1 − 2C ) t0j . (9) Из полученного соотношения исключаем t−j1 с помощью (4.8): t0j +1 = 2C t1j + (1 − 2C ) t0j . (10) Формула (10) служит для определения значений температуры на левой границе (в точке x0 ). Из формул (5), (6) и (10) следует, что для получения значений сеточной функции на верхнем j+1 временном слое в i-ом узле необходимо знать три значения функции t ( x,τ ) на нижнем j-ом временном слое, а именно, в узлах i – 1, i, i + 1. Рассмотрим случай, когда тонкая пластина толщиной 2R , имевшая в начальный момент времени постоянную температуру tн , охлаждается в среде, температура которой tж = 0. На обеих сторонах пластины происходит конвективный теплообмен с постоянной интенсивностью α . Необходимо рассчитать температурное поле в пластине как функцию t ( x, τ ) . Используя симметрию температурного поля, поместим начало координат на ось пластины. Необходимо найти функцию t ( x, τ ) , удовлетворяющую следующим условиям: уравнению теплопроводности: tτ′ = a t xx′′ ; x ∈ ( − R, R ) , τ ∈ ( 0, τ m ) ; (11) начальному условию: t ( x, 0 ) = tн , x ∈ [ − R, R ] ; (12) граничным условиям: t x′ ( R, τ ) = − Условие t x′ ( − R, τ ) = α α t ( R, τ ) − tж  ; t x′ ( − R, τ ) = t ( − R, τ ) − tж  . (13) λ λ α t ( − R, τ ) − tж  можно заменить условием t x′ ( 0, τ ) = 0 , что λ означает отсутствие теплового потока в плоскости x = 0 и следует из симметрии температурного поля платины при заданных условиях. Применение при решении задачи двухслойной явной схемы, приводит к использованию для вычисления значений сеточной функции формулы (5). Определение значений температуры на левой границе (в точке x0 ) происходит с помощью формулу (10). Для граничных условий третьего рода производную в точке xN ( x = R ) заменяют конечно-разностным аналогом со вторым порядком точности: t Nj +1 − t Nj −1 α = − ( t Nj − tж ) . 2h λ (14) Следовательно, значение температуры, записанное в фиктивном узле i = N + 1, равно: t Nj +1 = t Nj −1 − 2hα λ Записываем уравнение (4.5) для i = N : ( tNj − t ) . ж (15) t Nj +1 = C ( t Nj −1 + t Nj +1 ) + (1 − 2C ) t Nj . (16) Из полученного соотношения исключаем t Nj +1 с помощью (15): 2Chα  j 2Chα  t Nj +1 = 2C t Nj −1 + 1 − 2C − tN + tж . λ  λ  (17) Формула (17) служит для определения значений температуры на правой границе (в точке xN ). Таким образом, для вычисления значений сеточной функции на верхнем временном слое необходимо использовать формулы (5), (10) и (17). В математическом пакете MathCAD численное решение смешанной задачи для уравнений теплопроводности с одной пространственной переменной осуществляется с помощью встроенной функции Pdesolve. Обращение к встроенному интегратору в этом случае выглядит следующим    0    образом: u:=Pdesolve  u, x,   , t,  tMax  , [ xpts ] , [ tpts]  . xMax       Дифференциальное уравнение для u(x,t) и краевые условия вводятся в привычной математической форме между служебным словом Given («дано») и обращением к интегратору pdesolve («решить уравнение в частных производных»). При записи уравнения начальных и граничных условий для указания частной производной используется буквенный подстрочный индекс. Буквенный индекс создают нажатием клавиши «десятичная точка» (точка лат. алф.). В перечне аргументов функции pdesolve указывают: 1) имя искомой функции u; 2) пространственную переменную x; 3) вектор-столбец (0, xMax), содержащий граничные значения координаты x; 4) переменную t; 5) вектор-столбец (0, tMax), содержащий граничные значения времени t; 6) xpts необязательный параметр, задающий число точек пространственной дискретизации; 7) tpts дискретизации. необязательный параметр, задающий число точек временной 4. ТЕМПЕРАТУРНОЕ ПОЛЕ В ПРОЦЕССЕ ОХЛАЖДЕНИЯ ШАРА ИЛИ БЕСКОНЕЧНО ДЛИННОГО ЦИЛИНДРА Рассмотрим задачи о симметричном нестационарном охлаждении шара или бесконечно длинного цилиндра в среде с постоянной температурой и при постоянном коэффициенте теплоотдачи. Пусть бесконечно длинный цилиндр радиусов R с постоянной начальной температурой tн охлаждается в среде с постоянной температурой tж при условии, что на поверхности тела происходит конвективный теплообмен (задача симметричная). Требуется найти температуру в любой точке через время τ после начала охлаждения. Если в начальный момент времени температура в каждой точке зависит только от ее расстояния r до оси цилиндра, то и в последующем температура t будет зависеть от r и времени τ . Тепловой поток при этом всегда направлен по радиусам цилиндра. Таким образом, t = t ( r , τ ) . Преобразуя общее уравнение теплопроводности к цилиндрическим координатам (учитывая, что функция t не зависит от ϕ и z ), получим уравнение: 1   tτ′ = a  trr′′ + tr′  . r   (1) Составление математической модели в условиях, когда на поверхности цилиндра происходит конвективный теплообмен, приводит к смешанной задаче для нестационарного уравнения теплопроводности (1) с начальным и граничными условиями: t ( r , 0 ) = tн , r ∈ [ 0, R ] ; (2) tr′ ( 0, τ ) = 0, τ ∈ [ 0, τ m ] ; (3) −λ tr′ ( R, τ ) = α ( t ( R, τ ) − tж ) , τ ∈ [ 0, τ m ] . Подставляя в уравнение (1) соотношения r = ih , trr′′ ≈ (4) ti j+1 − 2ti j + tij−1 tij+1 − ti j−1 ′ , t ≈ и r h2 2h ti j +1 − ti j tτ′ ≈ , получаем разностное уравнение: k  ti j+1 − 2tij + ti j−1 ti j+1 − ti j−1  tij +1 − ti j = a +  , 1 ≤ i ≤ N − 1, 0 ≤ j ≤ M − 1. 2ih2  k h2  (5) Выражая из уравнения (5) tij +1 , получим формулу для вычисления значений сеточной функции: 1 1   tij +1 = C 1 +  tij−1 + (1 − 2C ) ti j + C 1 −  ti j+1 ,  2i   2i  где C = (6) ak , 1 ≤ i ≤ N − 1, 0 ≤ j ≤ M − 1. h2 Для реализации граничных условий третьего рода в точке rN (r = R) записываем уравнение (6) для i = N : 1  j 1  j   j t Nj +1 = C 1 +  t N +1 + (1 − 2C ) t N + C 1 −  t N −1 .  2N   2N  (7) t Nj +1 − t Nj −1 Из аппроксимации граничных условий третьего рода −λ = α ( t Nj − tж ) 2h получим значение температуры, записанное в фиктивном узле i = N + 1 , t Nj +1 = t Nj −1 − 2hα λ (t j N − tж ) (8). Из соотношения (7) исключаем t Nj +1 с помощью формулы (8):  2Chα  1   j 2Chα  1  t Nj +1 = 2C t Nj −1 + 1 − 2C − t . 1 +  tN + 1 + λ  2N  λ  2 N  ж  (9) Формула (9) служит для определения значений температуры на правой границе. Для нахождения решения на верхнем временном слое в узле i = 0 нельзя применять формулу (6) поскольку существует особенность реализации граничных условий второго рода при r = 0. В этом случае устремляя r к нулю и раскрывая по правилу Лопиталя неопределенность tr′ , получим для уравнения (1) в точке r = 0 r следующий вид: (10) tτ′ = 2atrr′′ . Используя соотношения tτ′ ≈ ti j +1 − ti j k и trr′′ ≈ ti j+1 − 2tij + ti j−1 h2 для подстановки в уравнение (10), получаем разностное уравнение: tij +1 − ti j t j − 2tij + ti j−1 = 2a i +1 . k h2 Записываем уравнение (11) для i = 0 : (11) t0j +1 = 2C ( t−j1 + t1j ) + (1 − 4C ) t0j . (12) t1j − t−j1 = 0 получим 2h Из аппроксимации граничных условий второго рода значение температуры, записанное в фиктивном узле t−j1 = t1j , Из соотношения (12) исключаем t−j1 t0j +1 = 4C t1j + (1 − 4C ) t0j . (13) Формула (13) служит для определения значений температуры на левой границе (в точке r = 0 ). Таким образом, для вычисления значений сеточной функции на верхнем временном слое в данном случае необходимо использовать формулы (6), (9) и (13). 1 4 Условие устойчивости для предложенной разностной схемы C < . С учетом особенности реализации граничных условий на левой границе в точке r = 0, получить решение поставленной задачи в системе MathCAD с применение встроенной функции pdesolve невозможно. 5. ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧ ТЕПЛОПРОВОДНОСТИ МЕТОДОМ ПРОГОНКИ 5.1. ДВУХСЛОЙНАЯ НЕЯВНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ Рассмотрим теплопередачу через плоскую бесконечную пластину (или изолированный стержень). Дана пластина толщиной L, м . В начале процесса нагревания температура по всему объему пластины была одинаковой и равной t0 , °C. В течение всего процесса на одной границе (стороне поверхности) поддерживается постоянная температура tl , °C, на другой t p , °C. Требуется рассчитать температурное поле пластины через τ m , c . Пластина изготовлена из стали: λ , Вт ( м ⋅ о С ) , с, Дж ( кг ⋅ о С ) , ρ , кг м3 . В данной задаче используем обозначение: t ( x, τ ) − температура, °C. Таким образом, нужно найти функцию t ( x, τ ) , удовлетворяющую следующим условиям: уравнению теплопроводности: ρ c ∂t ∂ 2t = λ 2 ; x ∈ ( 0, L ) , τ ∈ ( 0, τ m ) ; ∂τ ∂x (1) начальному условию: t ( x, 0 ) = t0 , x ∈ [ 0, L] ; (2) граничным условиям: t ( 0, τ ) = tl , t ( L, τ ) = t p , τ ∈ [0, τ m ]. (3) Используя двухслойную неявную схему, имеем разностное уравнение: ρc ti j +1 − tij t j +1 − 2tij +1 + ti j−+11 = λ i +1 , 1 ≤ i ≤ N − 1, 0 ≤ j ≤ M − 1. k h2 (4) Поле температуры на новом временном слое определяется из системы (5) Ai ⋅ tij++11 − Bi ⋅ ti j +1 + Ci ⋅ ti j−+11 = Fi , где Ai = Ci = λ k , Bi = 2λ k + ρ ch 2 Fi = − ρ ch 2tij . Подставляя основное соотношение прогонки с уменьшенным на единицу нижним индексом tij−+11 = α i −1 ⋅ ti j +1 + βi −1 в уравнение (5), преобразуем трехточечное уравнение в двухточечное: Ai ⋅ tij++11 − Bi ⋅ tij +1 + Ci ⋅ (α i −1tij +1 + β i −1 ) = Fi . Следовательно, tij +1 = Ai Cβ −F ti j++11 + i i −1 i . Bi − Ciα i −1 Bi − Ciα i −1 Используя основное соотношение прогонки tij−+11 = α i −1 ⋅ ti j +1 + β i −1 , получаем формулы для расчета прогоночных коэффициентов: αi = Ai C ⋅β − F , β i = i i −1 i . Bi − Ci ⋅ α i −1 Bi − Ci ⋅ α i −1 (6) Формулы для определения начальных прогоночных коэффициентов α1 , β1 находим из левого граничного условия t1j +1 = α1t2j +1 + β1 = tl . Следовательно, α1j +1 = 0, β1j +1 = tl . (7) Значение температуры на правой границе определяем, исходя из правого граничного условия: t Nj +1 = t pj +1 (8) При решении задачи методом прогонки с применением двухслойной неявной схемы используем формулы (5), (6), (7) для расчета прогоночных коэффициентов, формулы для расчета значений температуры на правой границе (8) и основное соотношение прогонки. Замечание. Для успешного применения метода прогонки необходимо выполнение условия исключающего быстрый рост погрешности округления (α i ) < 1, i = 1, N − 1 и необращение знаменателей прогоночных коэффициентов в нуль. Это позволяет реализовать алгоритм решения задачи с применением данного метода, учитывая безусловную устойчивость предложенной неявной разностной схемы. 5.2. КОНЕЧНО-РАЗНОСТНЫЕ АППРОКСИМАЦИИ ГРАНИЧНЫХ УСЛОВИЙ ВТОРОГО И ТРЕТЬЕГО РОДА Граничные условия второго рода : −λ ∂T ( x,τ ) = q1 , τ > 0; ; ∂x x =0 −λ ∂T ( x,τ ) = q2 , τ > 0 . ∂x x=L (9) В случае нагревания: q1 > 0, q2 < 0 . В случае охлаждения: q1 < 0, q2 > 0 . При использовании неявной разностной схемы из левого граничного условия определяют первые прогоночные коэффициенты из соотношения: t1j +1 = α1 ⋅ t2j +1 + β1 . Проведем дискретизацию граничных условий второго рода с погрешностью аппроксимации O ( h ) . Для определения начальных прогоночных коэффициентов α1 , β1 находим из аппроксимации левого граничного условия −λ t1j +1 − t0j +1 qh = q1 , t0j +1 = t1j +1 + 1 h λ и разностного уравнение (4) при i = 1 , получаем соотношение: ρ ch 2 ( t1j +1 − t1j ) = λ kt2j +1 − λ kt1j +1 + khq1 , (10) Выражая из равенства (10) t1j +1 , получаем формулы для коэффициентов: α1j +1 = Значение аппроксимации температуры правого λk ρ ch 2 + λ k на , β1j +1 = правой граничного khq1 + ρ ch 2t1j . ρ ch 2 + λ k границе условия −λ определяем, (11) исходя t Nj ++11 − t Nj +1 qh = q2 , t Nj ++11 = t Nj +1 − 2 h λ из и разностного уравнение (4) при i = N , получаем соотношение ρ ch 2 ( t Nj +1 − t Nj ) = −λ kt Nj +1 + λ kt Nj +−11 − khq2 . (12) Подставляя основное соотношение t Nj +−11 = α N −1 ⋅ t Nj +1 + β N −1 прогонки в уравнение (12), получаем формулы для расчета значений температуры на правой границе: t j +1 N λ k β N −1 + t Nj ρ ch 2 + kh 2 q2 = , ρ ch 2 + λ k (1 − α N −1 ) (13) Проведем дискретизацию граничных условий второго рода с погрешностью аппроксимации O ( h 2 ) . ∂t ∂ 2t Пусть на границе выполняется уравнение: ρ c = λ 2 ∂τ ∂x Разложим функцию t ( x,τ ) в ряд Тейлора в окрестности точки x = 0 до членов второго порядка относительно h : t j +1 2 j +1 1 ≈t ∂t +h ∂x j +1 h 2 ∂ 2t + 2 ∂x 2 x=0 j +1 x=0 Используя уравнения теплопроводности, получаем: t2j +1 ≈ t1j +1 + h ∂t Тогда ∂x j +1 x =0 t2j +1 − t1j +1 ρ ch ∂t = − h 2λ ∂τ ∂t ∂x j +1 + x =0 ρ ch 2 ∂t 2λ ∂τ j +1 . x =0 j +1 . Из граничного условия (9) следует: x =0 t2j +1 − t1j +1 h ρ c t1j +1 − t1j q − =− 1. h 2λ k λ Выражая из полученного соотношения t1j +1 , получаем формулы для коэффициентов: α1j +1 = 2λ k 2khq1 + ρ ch 2t1j j +1 β , = . 1 ρ ch 2 + 2λ k ρ ch 2 + 2λ k (14) Определяем t Nj +1 , используя правое граничное условие. Разложим функцию t ( x,τ ) в ряд Тейлора в окрестности точки x = L до членов второго порядка относительно h : t j +1 N −1 ≈t j +1 N ∂t −h ∂x j +1 h 2 ∂ 2t + 2 ∂x 2 x=L j +1 x=L Используя уравнение теплопроводности, получаем: t Nj +−11 ≈ t Nj +1 − h ∂t Тогда ∂x j +1 x=L t Nj +1 − t Nj +−11 ρ ch ∂t = + h 2λ ∂τ j +1 ∂t ∂x + x=L ρ ch 2 ∂t 2λ ∂τ j +1 . x=L j +1 . Из граничного условия (9) следует: x= L t Nj +1 − t Nj +−11 h ρ c t Nj +1 − t Nj q + =− 2 . h 2λ k λ Подставляя основное соотношение прогонки t Nj +−11 = α N −1 ⋅ t Nj +1 + β N −1 в уравнение (12), получаем формулы для расчета значениий температуры на правой границе: 2k β N −1 + t Nj ρ ch 2 − 2kh 2 q2 , ρ ch 2 + 2λ k (1 − α N −1 ) t Nj +1 = (15) Граничные условия третьего рода: λ ∂t ( x,τ ) = γ 1 (t − tg ) , τ > 0 ; ∂x x =0 −λ ∂t ( x,τ ) = γ 2 (t − tg ) , τ > 0 , ∂x x = L (16) где t g − температура среды. Проведем дискретизацию граничных условий (16) с погрешностью аппроксимации O ( h 2 ) . Разложим функцию t ( x,τ ) в ряд Тейлора в окрестности точки x = 0 до членов второго порядка относительно h : t j +1 2 j +1 1 ≈t ∂t +h ∂x j +1 h 2 ∂ 2t + 2 ∂x 2 x=0 j +1 x=0 Используя уравнения теплопроводности, получаем: t ∂t Тогда ∂x j +1 x =0 j +1 2 j +1 1 ≈t t2j +1 − t1j +1 ρ ch ∂t = − h 2λ ∂τ ∂t +h ∂x j +1 ρ ch 2 ∂t + 2λ ∂τ x =0 j +1 . x =0 j +1 . Из граничного условия (16) следует: x =0 t2j +1 − t1j +1 h ρ c t1j +1 − t1j γ 1 j +1 − = ( t1 − t g ) . h 2λ k λ Выражая из полученного соотношения t1j +1 , получаем формулы для коэффициентов: α j +1 1 2khγ 1t g + ρ ch2t1j 2λ k j +1 = , β1 = . ρ ch 2 + 2k ( λ + hγ 1 ) ρ ch 2 + 2k ( λ + hγ 1 ) (18) Определяем t Nj +1 , используя правое граничное условие. Разложим функцию t ( x,τ ) в ряд Тейлора в окрестности точки x = L до членов второго порядка относительно h : t j +1 N −1 ≈t j +1 N ∂t −h ∂x j +1 h 2 ∂ 2t + 2 ∂x 2 x=L j +1 x=L Используя уравнение теплопроводности, получаем: t Тогда ∂t ∂x j +1 = x=L j +1 N −1 ≈t j +1 N t Nj +1 − t Nj +−11 ρ ch ∂t + h 2λ ∂τ ∂t −h ∂x j +1 ρ ch 2 ∂t + 2λ ∂τ x=L j +1 . x=L j +1 . Из граничного условия (16) следует: x= L t Nj +1 − t Nj +−11 h ρ c t Nj +1 − t Nj γ + = − 2 ( t Nj +1 − t g ) . h 2λ k λ Подставляя основное соотношение прогонки t Nj +−11 = α N −1 ⋅ t Nj +1 + β N −1 в уравнение (12), получаем формулы для расчета значений температуры на правой границе: t j +1 N = 2k λβ N −1 + t Nj ρ ch 2 + 2khγ 2t g ρ ch 2 + 2λ k (1 − α N −1 ) + 2khγ 2 , (19) 5.3. ТЕМПЕРАТУРНОЕ ПОЛЕ В МНОГОСЛОЙНОЙ ПЛАСТИНЕ Во многих практических приложениях проводиться оценка температурных полей в многослойных деталях конструкций. Рассмотрим процесс теплопроводности в теле, представляющем собой совокупность двух пластин с различными теплофизическими характеристиками. Математическая модель процесса теплопроводности в двухслойной пластине включает в себя линейные однородные дифференциальные уравнения для каждого из слоев: ∂T1 ( x,τ ) λ ∂ 2T ( x,τ ) = 1 ⋅ 1 2 , 0 < x < x* ; ∂τ c1 ρ1 ∂x (20) ∂T2 ( x,τ ) λ ∂ 2T2 ( x,τ ) * = 2 ⋅ , x < x < L; ∂τ c2 ρ 2 ∂x 2 (21) Tv ( x, 0 ) = 0, x ∈ [ 0, L ] ; v = 1, 2; (22) начальные условия: условия на внешних границах: T1 ( 0,τ ) = Tl (23) T2 ( L,τ ) = Tp (24) граничные условия на стыке слоев: ( ) ( ) T1 x* ,τ = T2 x* ,τ , λ1 (25) ∂T1 ( x,τ ) ∂T ( x,τ ) = λ2 2 . ∂x ∂x x = x* x = x* (26) Для решения задач (20)–(26) применяем метод сеток с использованием конечно-разностных аппроксимаций прямоугольную (метод конечных разностей). пространственно-временную { Вводим сетку } Gxτ = xi = ih, i = 1, N ; τ j = jk , j = 0, M , при этом число узлов по пространственной N = N1 + N 2 + 1 , координате определяется соотношением где Nv − количество промежутков, на которые разбивается первый или второй слои. При получении трехслойной неявной разностной схемы подставляем формулы соответствующих разностных аналогов частных производных по времени и пространству в гиперболическое уравнение теплопроводности для каждого слоя. В этом случае разностное уравнение имеет вид: Ti j +1 − Ti j Ti +j1+1 − 2Ti j +1 + Ti −j1+1 = av , k h2 где av = (27) λv , v = 1, 2 , i = 2, N − 1, j = 0, M − 1 . cv ρ v Полученную систему приводим к системе линейных алгебраических уравнений, имеющей трехдиагональную структуру и решаемой методом прогонки: Ai ⋅ Ti +j1+1 − Bi ⋅ Ti j +1 + Ci ⋅ Ti −j1+1 = Fi , (28) где Ai = Ci = av k , Bi = h2 + 2av k , Fi = −Ti j h 2 . Подставляя основное соотношение прогонки Ti j +1 = α i ⋅ Ti +j1+1 + βi с уменьшенным на единицу нижним индексом в уравнение (13), получаем формулы для расчета прогоночных коэффициентов: αi = Ai C ⋅β − F , β i = i i −1 i . Bi − Ci ⋅ α i −1 Bi − Ci ⋅ α i −1 (29) Формулы для определения начальных прогоночных коэффициентов α1 , β1 находим из левого граничного условия α1j +1 = 0, β1j +1 = Tl . (30) Значение температуры на правой границе определяем, исходя из правого граничного условия TNj +1 = Tp . (31) Прогоночные коэффициенты в точке контакта двух сред выводим из аппроксимации граничного условия четвертого рода, используя разложение сеточной функции в ряд Тейлора в окрестности точки сопряжения x = x* до членов второго порядка относительно h : j +1 i* −1 T j +1 i* +1 T j +1 ≈ Ti* j +1 ≈ Ti* ∂t −h ∂x ∂t +h ∂x j +1 j +1 x = x* h 2 ∂ 2t + 2 ∂x 2 x=x j +1 2 j +1 x = x* 2 h ∂t + 2 ∂x 2 ; * . x = x* Используя уравнения теплопроводности, получаем: Ti*j−+11 Ti*j++11 j +1 ≈ Ti* j +1 ≈ Ti* ∂t −h ∂x ∂t +h ∂x j +1 x = x* j +1 x = x* h 2 ∂t + 2a1 ∂τ h 2 ∂t + 2a2 ∂τ j +1 ; x = x* . j +1 . x = x* j +1 j ∂T Ti* − Ti* Тогда, учитывая , имеем: = ∂τ k ∂T ∂x ∂T ∂x j +1 = Ti*j +1 − Ti*j−+11 h x = x* j +1 = x = x* j +1 i* +1 T − Ti*j +1 h + j +1 j h Ti* − Ti* ; 2a1 k j +1 j h Ti* − Ti* − . 2 a2 k . Подставляя полученные соотношения в условия сопряжения, имеем: j +1 j j +1 j  Ti *j +1 − Ti *j−+11  Ti *j++11 − Ti *j +1 h Ti * − Ti *  h Ti * − Ti *  λ1  +  = λ2  −  h 2a1 k h 2a2 k     Используя прогоночное соотношение Ti j−+11 = α i −1 ⋅ Ti j +1 + βi −1 и выражая * * * * (32) из соотношения (32) Ti j +1 , получаем формулы для прогоночных коэффициентов в точке * x = x* : α i j +1 = * βi j +1 * = 2a1a2 k λ2 , 2a1a2 k λ1 1 − α i *−1 + λ2  + h 2 ( λ1a2 + λ2 a1 )   ( (33) ) 2a1a2 k λ1β i *−1 + h 2 ( λ1a2 + λ2 a1 ) Ti *j ( ) 2a1a2 k λ1 1 − α i *−1 + λ2  + h 2 ( λ1a2 + λ2 a1 )   . (34) Полученные формулы используются при реализации прямой и обратной прогонки для описания процедуры расчета значений функции на верхнем (промежуточном) временном слое. При этом алгоритм расчета предусматривает сначала определение прогоночных коэффициентов для первой среды, затем на границе i * , используя соотношения (33), (34), далее задаются прогоночные коэффициенты для второй среды. Определение температуры на правой границе происходит по формуле (31) в соответствии с временным слоем. Расчет температурного поля Ti j +1 проводиться с применением основного соотношения прогонки. Выполнение округления (α i условия ) < 1, i = 1, N − 1 исключающего и необращение быстрый рост знаменателей погрешности прогоночных коэффициентов в нуль обеспечивает успешное применение метода прогонки и позволяет реализовать алгоритм решения задачи с применением данного метода, учитывая безусловную устойчивость предложенной неявной разностной схемы с погрешностью аппроксимации O (τ + h 2 ) . 6. ТЕМПЕРАТУРНОЕ ПОЛЕ В ПРОЦЕССЕ ОХЛАЖДЕНИЯ НЕОГРАНИЧЕННОЙ ПЛАСТИНЫ С УЧЕТОМ ИЗЛУЧЕНИЯ Предположим, что тонкая пластина толщиной 2R с начальной температурой Tн охлаждается жидкостью, имеющей постоянную температуру Tж , при коэффициенте теплоотдачи α. Пластина участвует в радиационном взаимодействии с окружающими телами, имеющими температуру, равную Tж . Плотность теплового ( 4 ) потока результирующего излучения можно записать как qрез = ε eqσ T ( R, τ ) − Tж4 , где ε eq − степень черноты обменивающейся излучением системы, σ − постоянная Стефана–Больцмана. Необходимо рассчитать температурное поле в пластине как функцию T ( x, τ ) . В математической постановке этой задачи нужно найти функцию T ( x, τ ) , удовлетворяющую уравнению теплопроводности: ∂T ∂ 2T = a 2 ; x ∈ ( − R, R ) , τ ∈ ( 0, τ m ) ; ∂τ ∂x (1) начальному условию: T ( x, 0 ) = T0 , x ∈ [ − R, R ] ; граничному условию второго рода: −λ нелинейному граничному условию: −λ (2) ∂T ( x,τ ) = 0; ∂x x =0 (3) ∂T 4 ( R, τ ) = γ (T ( R, τ ) − Tж ) + ε eqσ T ( R, τ ) − Tж4 (4) ∂x ( ) При реализации двухслойной явной схемы значения сеточной функции определяются с помощью формулы (5, п 3). Значения температуры на левой границе (в точке x0 = 0 ) находятся с помощью формулу (10, п. 3). Для граничных условий (5.2) производную в точке xN ( x = R ) заменяют конечно-разностным аналогом со вторым порядком точности: ) ( 4 4 TNj+1 − TNj−1 −λ = α (TNj − Tж ) + εσ (TNj ) − (Tж ) . 2h (5) Выражая отсюда TNj+1 , получаем: TNj+1 = TNj−1 − 2hα λ (TNj − Tж ) − 2hεσ λ ( (T ) − (T ) ) . j 4 N 4 ж (6) Используя соотношение (6) для исключения TNj+1 в разностном уравнении (16, п.3), получим формулу для нахождения значений функции на правой границе: TNj +1 = 2CTNj−1 − 2Chα λ (TNj − Tж ) − 2Сhεσ λ ((T ) − (T ) ) + T (1 − 2C ). j 4 N 4 ж j N (7) Следовательно, для вычисления значений сеточной функции на верхнем временном слое необходимо использовать формулы (5, п.3), (10, п.3) и (7). Условие 1 2 устойчивости для явной разностной схемы C ≤ . В математическом пакете MathCAD численное решение смешанной задачи для уравнений теплопроводности с одной пространственной переменной и нелинейными граничными условиями можно реализовывать с помощью встроенной функции Pdesolve. Дифференциальное уравнение для искомой функции и краевые условия вводятся в привычной математической форме после служебного слова Given («дано»): Given Tτ ( x, τ ) = a ⋅ Txx ( x, τ ) T ( x, 0 ) = To Tx ( 0, τ ) = 0 ( -λ ⋅ Tx ( R, τ ) = α ⋅ ( T ( R, τ ) − Tc ) + ε ⋅ σ ⋅ T ( R, τ ) − Tc 4 4 ) Обращение к встроенному интегратору в этом случае выглядит следующим образом:    0   0  T:=Pdesolve  T, x,  ,τ ,  , [ xpts ] , [τ pts ]  .    xMax   τ Max    В перечне аргументов функции Pdesolve указывают: 1) имя искомой функции; 2) пространственную переменную; 3) вектор-столбец (0, xMax), содержащий граничные значения координаты x; 4) переменную τ ; 5) вектор-столбец (0, τMax ), содержащий граничные значения времени τ ; 6) xpts необязательный параметр, задающий число точек пространственной дискретизации; 7) τpts необязательный параметр, задающий число точек временной дискретизации. Для получения численного решения поставленной задачи с реализацией неявной разностной схемы для левой границы и внутренних узлов, а на правой границе, учитывая нелинейность граничных условий, используем явную разностную схему (явно-неявная разностная схема). Используя двухслойную неявную схему, имеем разностное уравнение: Ti j +1 − Ti j T j +1 − 2Ti j +1 + Ti −j1+1 = a i +1 , 1 ≤ i ≤ N − 1, 0 ≤ j ≤ M − 1. k h2 Умножая обе части равенства на множитель kh 2 , получаем соотношение: ak ⋅ Ti +j1+1 − ( h 2 + 2ak ) ⋅ Ti j +1 + ak ⋅ Ti −j1+1 = − h 2Ti j , (8) Поле температуры на новом временном слое определяется из системы Ai ⋅ Ti +j1+1 − Bi ⋅ Ti j +1 + Ci ⋅ Ti −j1+1 = Fi , где Ai = Ci = ak , Bi = 2ak + h 2 Fi = −h 2Ti j . Подставляя основное соотношение прогонки с уменьшенным на единицу нижним индексом Ti −j1+1 = α i −1 ⋅ Ti j +1 + βi −1 , преобразуем трехточечное уравнение в двухточечное: Ai ⋅ Ti+j1+1 − Bi ⋅ Ti j+1 + Ci ⋅ (α i−1Ti j+1 + βi−1 ) = Fi . Из соотношения Ti j +1 = Ai Cβ −F Ti +j1+1 + i i −1 i , получаем формулы для расчета Bi − Ciα i −1 Bi − Ciα i −1 прогоночных коэффициентов: αi = Ai , Bi − Ci ⋅α i−1 βi = Ci ⋅ βi −1 − Fi Bi − Ci ⋅ α i −1 (9) Формулы для определения начальных прогоночных коэффициентов α1 , β1 находим из левого граничного условия. Из предположения, что решение непрерывно записываем уравнение (8) для i = 1: ak ⋅ T2 j +1 − ( h 2 + 2ak ) ⋅ T1 j +1 + ak ⋅ T0 j +1 = − h 2T1 j . (10) Для исключения из полученного уравнения фиктивного узла T0j +1 используем граничные условия второго рода, заменяя производную в точке x1 ( x = 0 ) конечноразностным аналогом со вторым порядком точности: T2j +1 − T0 j +1 = 0. 2h Следовательно, подставляя значение температуры, записанное в фиктивном узле T0 j +1 = T2j +1 , получаем ak ⋅ T2 j +1 − ( h 2 + 2ak ) ⋅ T1 j +1 + ak ⋅ T2 j +1 = − h 2T1 j Из соотношения T1 j +1 = 2ak h 2T1 j j +1 T + , получаем формулы для 2 h 2 + 2ak h 2 + 2ak определения начальных прогоночных коэффициентов α1 , β1 : (11) α1 = 2ak h 2T1 j , β = 1 h 2 + 2ak h 2 + 2ak (12) Значение температуры на правой границе определяем, исходя из правого граничного условия. Используя двухслойную явную схему, имеем разностное уравнение: TNj +1 − TNj T j − 2TNj + TNj−1 = a N +1 , 1 ≤ i ≤ N − 1, 0 ≤ j ≤ M − 1. k h2 Тогда TNj +1 = C (TNj−1 + TNj+1 ) + (1 − 2C ) TNj , где C = (13) ak − число Куранта. h2 Для исключения из полученного уравнения фиктивного узла TNj+1 используем нелинейное граничное условие, заменяя производную в точке xN ( x = R ) конечноразностным аналогом со вторым порядком точности: −λ ) ( 4 4 TNj+1 − TNj−1 = γ (TNj − Tж ) + εσ (TNj ) − (Tж ) . 2h (14) Выражая отсюда TNj+1 , получаем: TNj+1 = TNj−1 − 2hγ λ (TNj − Tж ) − 2hεσ λ ( (T ) − (T ) ) . j 4 N 4 (15) ж Исключая из соотношения (13) фиктивный узел, используя формулу (15), получаем формулу для нахождения значений сеточной функции на правой границе: TNj +1 = 2CTNj−1 − 2Chγ λ (TNj − Tж ) − 2Сhεσ λ ( (T ) − (T ) ) + T j 4 N 4 ж j N (1 − 2C ) . (16) 1 2 Условие устойчивости для предложенной разностной схемы C ≤ . 7. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ ДЛЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С ЧАСТНЫМИ ПРОИЗВОДНЫМИ ГИПЕРБОЛИЧЕСКОГО ТИПА К уравнениям гиперболического типа приводят задачи связанные с процессами колебаний, например, задача о колебаниях струны, мембраны, газа, электромагнитных колебаниях и др. Среди задач теплопроводности в рамках данной тематики выделяют моделирование высокоинтенсивных процессов нагрева тел (плазменной, контактных лазерной соединений обработке в материалов, электрических высокоинтенсивном установках и др.). нагреве Характерной особенностью процессов, описываемых такими уравнениями, является конечная скорость их распространения. 7.1. ТЕЛЕГРАФНЫЕ УРАВНЕНИЯ Телеграфные уравнения ̶ пара линейных дифференциальных уравнений, описывающих распределение напряжения и тока по времени и расстоянию в линии электрической связи. Уравнения получены Оливером Хевисайдом, разработавшим в 1880 г. модель линии электрической связи. Теория применима к линиям передачи электрического тока всех частот (телеграфные, телефонные, силовые линии электропередачи, линии передачи постоянного тока). Рассмотрим двухпроводную линию. Напряжение между проводами u ( x, t ) и ток I ( x, t ) в некоторой точке линии зависит от расстояния x этой точки до начала линии и времени t . Рис.1. Схематическое изображение элементарных компонентов линии электрической связи Для составления дифференциальных уравнений, которым должны удовлетворять функции u ( x, t ) и I ( x, t ) , выделим бесконечно короткий участок линии от точки x до точки x + ∆x ( x + dx ) со следующими параметрами: R − погонное сопротивление проводников, Ом м ; L − погонная индуктивность, Гн м ; С − погонная ёмкость между проводниками, Ф м ; G − погонная проводимость диэлектрического материала (изоляции) разделяющего два проводника, сименс м ; В соответствии с законом Кирхгофа, разность напряжений в начале и конце рассматриваемого участка линии равна сумме падения напряжения на активном сопротивлении и индуктивного падения напряжения: u − (u + ∂u ∂I dx) = R ⋅ Idx + L dx . ∂x ∂t (1) Изменение тока на этом участке обусловлено током утечки и током смещения: I − (I + ∂I ∂u dx) = G ⋅ udx + C dx . ∂x ∂t (2) Приводим уравнения (1), (2) к виду: ∂u ∂I + R⋅I + L = 0, ∂x ∂t (3) ∂I ∂u + G ⋅u + C =0. ∂x ∂t (4) Дифференцируя уравнение (3) по x , а уравнение (4) по t , получаем: ∂ 2u ∂I ∂2 I R L + + =0, ∂x 2 ∂x ∂x∂t (5) ∂2I ∂u ∂ 2u +G +C 2 = 0 . ∂x∂t ∂t ∂t (6) ∂2I ∂u ∂ 2u = −G − C 2 . ∂x∂t ∂t ∂t (7) ∂I ∂u = −G ⋅ u − C . ∂x ∂t (8) Из соотношения (6) имеем: Из уравнения (4) получаем: Подставляя соотношения (7) и (8) в уравнение (5), получаем телеграфное уравнение для u ( x, t ) : ∂ 2u RC + LG ∂u RG 1 ∂ 2u + ⋅ + u − =0, ∂t 2 LC ∂t LC LC ∂x 2 (9) Аналогично можно получить телеграфное уравнение для функции I ( x, t ) : ∂ 2 I RC + LG ∂I RG 1 ∂2I + ⋅ + I − = 0, ∂t 2 LC ∂t LC LC ∂x 2 Рассмотрим некоторые частные случаи: (10) 1) Бесконечная линия (граничные условия отсутствуют). Предполагается, что в начальный момент вдоль линии задано распределение напряжения и тока: u ( x, t ) t =0 = f ( x) , I ( x, t ) t =0 = ϕ ( x) Используя уравнения (3) и (4), получаем начальные условия для функции u ( x, t ) и I ( x, t ) : u ( x, t ) t =0 = f ( x),   ∂u G 1 = − f ( x) − ϕ ′( x);  C C  ∂t t =0 (11)  I ( x, t ) t =0 = ϕ ( x),   ∂I R 1 = − ϕ ( x) − f ′( x);  L L  ∂t t =0 (12) 2) Линия без потерь (в случае, когда сопротивление провода мало и провод изолирован, можно принять R = G = 0 ). Телеграфное уравнение обращается в одномерное волновое уравнение: ∂ 2u 1 ∂ 2u = . ∂t 2 LC ∂x 2 (13) Начальные условия для функции u ( x, t ) : u ( x, t ) t =0 = f ( x),   ∂u 1 = − ϕ ′( x);  C  ∂t t = 0 (14) 3) Линия без искажений ( RС = LG ). Проводя замену u ( x, t ) = v( x, t )e − Rt L для уравнения (9), получаем одномерное волновое уравнение: ∂ 2v 1 ∂ 2v = . ∂t 2 LC ∂x 2 (15) v( x, t ) t =0 = f ( x),   ∂v 1 = − ϕ ( x);  C  ∂t t =0 (16) Начальные условия: 4) Линия конечной длины. Если к началу линии подключено синусоидальное напряжение, то u ( x, t ) x = 0 = E sin ωt . Если в конце линия коротко замкнута, то u ( x, t ) x =l = 0 . Если в конец линии изолирован, то ∂u ( x, t ) = 0. ∂x x =l 7.2. ВОЛНОВОЕ УРАВНЕНИЕ 1) Уравнение колебаний струны (одномерное волновое уравнение): 2 ∂ 2u ( x, t ) P ( x, t ) 2 ∂ u ( x, t ) = a + В , 2 2 ∂t ∂x ρ (17) где u ( x, t ) − смещение точек струны от положения равновесия, PВ ( x, t ) − внешняя сила, действующая на струну перпендикулярно оси ОX и рассчитанную на единицу длинны, a = T0 ρ , T0 − величина силы натяжения, ρ − линейная плотность. Начальные условия: u ( x, t ) t =0 = ϕ0 ( x),   ∂u = ϕ1 ( x).   ∂t t = 0 (18) Граничные условия первого рода: u ( x, t ) x =0 = 0, u ( x, t ) x =l = 0 . 2) Уравнение колебаний мембраны(двумерное волновое уравнение): 2 ∂ 2u ( x, y , t ) ∂ 2u ( x, y, t )  PВ ( x, y, t ) 2  ∂ u ( x, y, t ) =a  + , + ∂t 2 ∂x 2 ∂y 2 ρ   (19) Начальные условия: u ( x, y, t ) t =0 = ϕ0 ( x, y ),   ∂u = ϕ1 ( x, y ).   ∂t t =0 (20) Граничные условия в случае закрепления мембраны на контуре L : u ( x, y , t ) x = L = 0 . 3) Уравнение распространения звуковых волн в идеальной жидкости (трехмерное волновое уравнение): 2 ∂ 2u ( x, y, z , t ) ∂ 2u ( x, y, z, t ) ∂ 2u ( x, y, z, t )  2  ∂ u ( x, y , z, t ) = a + +  , ∂t 2 ∂x 2 ∂y 2 ∂z 2   где u ( x, y, z , t ) − потенциал скорости в момент времени t , a = (21) P0 ρ − скорость распространения звука в изотермической жидкости, P0 − равновесное давление, ρ − равновесная плотность. Начальные условия: u ( x, y, z, t ) t =0 = ϕ0 ( x, y, z ),   ∂u = ϕ1 ( x, y, z ).   ∂t t = 0 (22) Граничные условия (если жидкость занимает объем, ограниченный твердой непроницаемой поверхностью S ): ∂u ∂n = 0. x=S 7.3. ТРЕХСЛОЙНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С ЧАСТНЫМИ ПРОИЗВОДНЫМИ ГИПЕРБОЛИЧЕСКОГО ТИПА В математическом пакете MathCAD численное решение смешанных задач для уравнений гиперболического осуществляется с помощью типа с одной встроенной пространственной функцией Pdesolve. переменной При этом дифференциальное уравнение преобразуется в систему двух уравнений первого порядка по времени введением новой зависимой переменной v(x,t). Обращение к встроенному интегратору в этом случае выглядит следующим образом:  u   u  0   0   v  :=Pdesolve  v  , x,  xMax  , t,  tMax  , [ xpts] , [ tpts]          Система дифференциальных уравнений для u(x, t), и v(x, t), краевые условия вводятся в привычной математической форме между служебным словом Given («дано») и обращением к интегратору pdesolve («решить уравнение в частных производных»). При записи уравнения, граничных и начальных условий для указания частной производной используется буквенный подстрочный индекс. Буквенный индекс создают нажатием клавиши «десятичная точка». В перечне аргументов функции pdesolve указывают: 1) вектор-столбец с именами искомых функций u и v в качестве компонентов; 2) пространственную переменную x; 3) вектор-столбец (0, xMax), содержащий граничные значения координаты x; 4) переменную t; 5) вектор-столбец (0, tMax), содержащий граничные значения времени t; 6) xpts – необязательный параметр, задающий число точек пространствен-ной дискретизации; 7) tpts – необязательный параметр, задающий число точек временной дискретизации. Для задачи вида:  ∂ 2u ( x, t ) ∂ 2u ( x, t ) ∂u ( x, t ) = − + f ( x, t ), 0 < x < L, t > 0,  2 ∂x 2 ∂t  ∂t u ( x, t ) t =0 = β ( x), 0 ≤ x ≤ L,   ∂u = 0, 0 ≤ x ≤ L,  ∂t t =0   ∂u ( x, t ) = ϕ ( x), t ≥ 0  ∂x x =0   ∂u ( x, t ) + u ( x, t ) x = L = ψ (t ), t ≥ 0,  ∂x x = L система дифференциальных уравнений и краевые условия вводятся после служебного слова Given («дано»): Given u t ( x, t ) = v ( x, t ) v t ( x, t ) = u xx ( x, τ ) − v ( x, τ ) − f(x,t) u ( x, 0 ) = β(x) v ( x, 0 ) = 0 u x ( 0, t ) = φ ( t ) u x ( L, t ) + u ( L, t ) = ψ ( t ) Обращение к встроенному интегратору в этом случае выглядит следующим образом:  u   u  0   0   v  :=Pdesolve  v  , x,  xMax  , t,  tMax  , [ xpts] , [ tpts]          Для решения уравнения методом прогонки необходимо построить сетку, покрывающую прямоугольник [0, L] ∪ [0, T ] . Координаты узлов сетки, образованные пересечением вертикальных и горизонтальных отрезков, определяются по формулам xi = ( i − 1) ⋅ h , t j = j ⋅ k , где h − шаг по пространству, h = времени, k = L , 1 ≤ i ≤ N ; k − шаг по N −1 T , 0≤ j≤M . M При получении трехслойной неявной разностной схемы подставляем формулы соответствующих разностных аналогов частных производных по времени и пространству в гиперболическое уравнение для каждого слоя: ∂ 2u uij +1 − 2uij + uij −1 = , ∂t 2 k2 ∂ 2u uij++11 − 2uij +1 + uij−+11 ∂u uij +1 − uij = = , ∂x 2 ∂t h2 k где f i j +1 = − h(i − 1) ⋅ ( cos ( t ) + sin ( t ) ) , i = 2, N , j = 0, M − 1 . В этом случае разностное уравнение имеет вид: uij +1 − 2uij + uij −1 uij++11 − 2uij +1 + uij−+11 uij +1 − uij = − + fi j +1 , k2 h2 k (23) Обе части уравнения (23) умножаем на k 2 h 2 : h 2 ( uij +1 − 2uij + uij −1 ) = k 2 ( uij++11 − 2uij +1 + uij−+11 ) − kh 2 ( uij +1 − uij ) + fi j +1k 2 h 2 , (24) Следовательно, k 2uij++11 − uij +1 [ 2k 2 + h2 + kh 2 ] + k 2uij−+11 = − fi j +1k 2 h 2 − uij [ 2h 2 + kh 2 ] + uij −1h 2 Полученную систему приводим к системе линейных (25) алгебраических уравнений, имеющей трехдиагональную структуру и решаемую методом прогонки: Ai ⋅ Ti +j1+1 − Bi ⋅ Ti j +1 + Ci ⋅ Ti −j1+1 = Fi , (26) где Ai = Ci = k 2 , Bi = 2k 2 + kh 2 + h 2 , Fi = − f i j +1k 2 h 2 − uij ( kh 2 + 2h 2 ) + h 2 ⋅ uij −1 . Подставляя основное соотношение прогонки с уменьшенным на единицу нижним индексом uij−+11 = α i −1 ⋅ uij +1 + β i −1 в уравнение (26), получаем формулы для расчета прогоночных коэффициентов: αi = Ai C ⋅β − F , β i = i i −1 i . Bi − Ci ⋅ α i −1 Bi − Ci ⋅ α i −1 (27) Для определения коэффициентов Fi на первом временном слое ( j = 0 ) вместо значения сеточной функции в фиктивном узле ui−1 используем соотношение ui−1 = ui0 , полученное из аппроксимации второго начального условия: Fi1 = − f i1k 2 h 2 − ui0 ( kh 2 + h 2 ) . (28) Формулы для определения начальных прогоночных коэффициентов α1 , β1 находим из аппроксимации левого граничного условия u2j +1 − u0j +1 = ϕ (t ), u0j +1 = u2j +1 − 2hϕ (t ) и разностного уравнение (25) при i = 1 , получаем 2h соотношение: u1j +1 ( kh 2 + h 2 + 2k 2 ) = 2k 2u2j +1 + u1j ( kh 2 + 2h 2 ) − h 2u1j −1 − 2hk 2ϕ (t ) + f1 j +1k 2 h 2 , (29) Выражая из равенства (29) u1j +1 , получаем формулы для коэффициентов: α1j +1 = f1 j +1h 2 k 2 − 2hk 2ϕ (t ) + u1j h 2 ( k + 2 ) − h 2u1j −1 2k 2 j +1 , = . β 1 (1 + k ) h 2 + 2k 2 (1 + k ) h2 + 2k 2 (30) Определяя коэффициент β1 на первом временном слое, используем формулу (30) при условии j = 0 и значение сеточной функции в фиктивном узле u1−1 = u10 : f11h 2 k 2 − 2hϕ (t ) + u10 h2 ( k + 1) . β = (1 + k ) h 2 + 2k 2 (31) 1 1 Значение температуры аппроксимации правого на правой граничного границе условия определяем, исходя u Nj ++11 − u Nj +−11 + u Nj +1 = ψ (t ), 2h u Nj ++11 = u Nj +−11 − 2hu Nj +1 + 2hψ (t ) , и разностного уравнение (25) при из т.е. i = N , получаем соотношение k 2 ( u Nj +−11 + 4h cos t − 2hu Nj +1 ) − u Nj +1 ( 2k 2 + h 2 + kh 2 ) + k 2u Nj +−11 = (32) = − f Nj +1k 2 h 2 − u Nj [ 2h 2 + kh 2 ] + u Nj −1h 2 Подставляя основное соотношение прогонки u Nj +−11 = α N −1 ⋅ u Nj +1 + β N −1 в уравнение (32), получаем формулы для расчета значений температуры на правой границе: u j +1 N = 2k 2 β N −1 + u Nj ( kh 2 + 2h 2 ) − h 2u Nj −1 + f Nj +1k 2 h 2 + 2hk 2ψ (t ) 2hk 2 + kh 2 + h 2 + 2k 2 (1 − α N −1 ) , (33) Определяя значение температуры в узле u1N , используем формулу (33) при условии j = 0 и значение сеточной функции в фиктивном узле u N−1 = u N0 : 1 N u = 2k 2 β N −1 + u N0 ( kh 2 + h 2 ) + f N1 k 2 h 2 + 2hk 2ψ (t ) 2hk 2 + kh 2 + h 2 + 2k 2 (1 − α N −1 ) ,. (34) 7.4. ГИПЕРБОЛИЧЕСКОЕ УРАВНЕНИЕ ТЕПЛОПРОВОДНОСТИ Ограниченность применения классического уравнения теплопроводности в случаях, когда необходимо учитывать инерционность процессов теплопереноса и конечную скорость распространения тепловых возмущений, обуславливает использование его гиперболической модификации. Гиперболическое уравнение теплопроводности непосредственно следует из обобщения гипотезы Фурье: q = −λ grad T − τ r ∂q , ∂τ (35) где q − вектор теплового потока, определяющий количество тепла, проходящего в единицу времени и отнесенного к единице площади изотермической поверхности; gradT − градиент температуры (вектор, направленный по нормали к изотермической поверхности в сторону возрастания температуры); теплопроводности; τ r − время релаксации теплового потока. λ− коэффициент Предполагая зависимость теплового потока от времени, отсутствие источников тепла и, используя классическое уравнение баланса тепла с применением обобщения гипотезы Фурье, после некоторых выкладок выводится соответствующее одномерное уравнение теплопроводности гиперболического типа: ∂T ( x,τ ) ∂ 2T ( x,τ ) ∂ 2T ( x,τ ) W ( x,τ ) τ r ∂W ( x,τ ) +τr = + + a , ∂τ ∂τ 2 ∂x 2 ∂τ cρ cρ (36) λ − коэффициент температуропроводности, м 2 с; T ( x, τ ) − температура тела cρ в точке x в момент времени τ , К , W ( x,τ ) − мощность внутренних источников тепла, λ − коэффициент теплопроводности, c − теплоемкость, плотность вещества, τ r − время релаксации теплового потока. Условия однозначности смешанной задачи для уравнения (36) имеют вид: где a = начальные условия: T ( x, 0 ) = ϕ ( x), x ∈ [ 0, L] ; ∂T ∂τ (37) = 0, x ∈ [ 0, L ] ; (38) граничные условия первого рода: T ( 0, τ ) = ψ 1 (τ ) , T ( L, τ ) = ψ 2 (τ ) , τ ∈ [0,τ m ]. (39) τ =0 Граничные условия второго рода для гиперболического теплопроводности (36) имеют два следующих вида: ∂T ( x,τ ) d   −λ = 1 + τ r  f (τ ), ∂x г  dτ  уравнения (40) если функция, характеризующая плотность теплового потока, f (τ ) − переменная величина и −λ если f (τ ) = f = const , ∂T ( x,τ ) = (1 + τ rδ (τ ) ) f , ∂x г где δ (τ ) − (41) дельта-функция Дирака, +∞ δ (τ ) = 0, τ ≠ 0; δ (0) = ∞; ∫ δ ( x)dx = 1. −∞ Граничные условия третьего рода записываются в виде: −λ ∂T ( x,τ ) d   = α 1 + τ r  T ( x,τ ) г − Tg (τ )  , ∂x г dτ    (42) если температура среды Tg (τ ) − переменная величина и −λ ∂T ( x,τ ) d   = α 1 + τ r  T ( x,τ ) г − α (1 + τ rδ (τ ) ) Tg , ∂x г dτ   (43) если температура среды Tg (τ ) = Tg = const. 7.5. РЕШЕНИЕ ЗАДАЧИ НЕСТАЦИОНАРНОЙ ТЕПЛОПРОВОДНОСТИ С УЧЕТОМ КОНЕЧНОЙ СКОРОСТИ РАСПРОСТРАНЕНИЯ ТЕПЛА Для создания технологических процессов с использованием методов обработки различных материалов концентрированными потоками энергии с получением более высоких физико-механических свойств широко используется лазерная и плазменная обработка различных материалов. В частности методы поверхностного термоупрочнения с использованием концентрированных пучков энергии (лазерных, плазменных, электронных) применяются в рамках совершенствования технологии ремонта подвижного состава железных дорог. Например, технологии лазерного упрочнения гребней бандажей колесных пар тягового подвижного состава. При исследовании высокоинтенсивных процессов нагрева тел (плазменной, лазерной обработке материалов, высокоинтенсивном нагреве контактных соединений в электрических установках и др.) используется гиперболическое уравнение теплопроводности. Для случая, когда поглощение энергии излучения моделируется объемным источником тепла, математическая модель процесса теплопроводности при воздействии концентрированных потоков энергии на тело включает в себя линейное неоднородное дифференциальное уравнение: ∂q (τ ) ∂ 2T ( x,τ ) ∂T ( x,τ ) ∂ 2T ( x,τ ) τ rγ +γ =λ + Aq0 (τ ) α e −α x + τ rα e −α x A 0 2 2 ∂τ ∂τ ∂x ∂τ (44) начальные условия: T ( x, 0 ) = T0 , x ∈ [ 0, L ] ; (45) ∂T ( x,τ ) = 0, x ∈ [ 0, L ] ; ∂τ τ =0 (46) ∂T ( x,τ ) = 0, τ ∈ [ 0, τ k ] ; ∂x x =0 (47) ∂T ( x,τ ) = 0, τ ∈ [ 0, τ k ] , ∂x x=L (48) граничные условия: где начальная температура тела T0 , K , L, мкм , τ k , нc , время релаксации теплового потока τ r , нс , коэффициент теплопроводности λ , Вт ( мкм ⋅ К ) , γ = с ρ , Вт ⋅ нс ( мкм3 ⋅ К ) , поглощательная способность материала A , коэффициент поглощения α , 1 м км , плотность потока энергии излучения q0 (τ ) , где максимальное значение плотности теплового потока qmax , Вт мкм 2 , длительность импульса τ d , нc . Для решения уравнения (44) методом конечных разностей необходимо построить сетку, покрывающую прямоугольник [0, L] ∪ [0,τ m ] . Координаты узлов сетки, образованные пересечением вертикальных и горизонтальных отрезков, определяются по формулам xi = i ⋅ h , τ j = j ⋅ k , где h − шаг по пространству, h = 0 ≤ i ≤ N ; k − шаг по времени, k = τm M L , N , 0 ≤ j ≤ M . Далее на построенной сетке вводим обозначение для сеточной функции: T ( xi , τ j ) = Ti , j . (49) Для реализации трехслойной явной разностной схемы подставляем формулы разностных аналогов частных производных по времени и пространству в гиперболическое уравнение теплопроводности, получаем разностное уравнение: Ti j +1 − Ti j Ti j +1 − 2Ti j + Ti j −1 λ Ti +j1 − 2Ti j + Ti +j1 +τ r = ⋅ + fi j +1 , 2 2 k k h γ где f i j +1 = (50) ∂q0 (τ )  1 −α x −α x  Aq0 (τ ) α e + τ rα e A . ∂τ  γ Выражая из соотношения (50) Ti j +1 , получим формулу для вычисления значений сеточной функции: −1 Ti j +1 = (τ r + k )  −τ rTi j −1 + C (Ti −j1 + Ti +j1 ) + ( k + 2τ r − 2C ) Ti j + f i j +1k 2  , где C = (51) λk 2 − число Куранта. Явная схема для гиперболического уравнения γ h2 теплопроводности является условно устойчивой и при решении необходимо учитывать следующее соотношение C ≤ 1 . Для вычисления значений функции на первом временном слое записываем формулу (51) для j = 0 : −1 Ti1 = (τ r + k )  −τ rTi −1 + C (Ti −01 + Ti +01 ) + ( k + 2τ r − 2C ) Ti 0 + f i1k 2  . (52) Заменим в начальном условии (47) производную по времени левой разностной производной с первым порядком точности: Ti 0 − Ti −1 = 0. k (53) Значение температуры, записанное в фиктивном узле j = −1, (54) Ti −1 = Ti 0 . Из полученного соотношения (52) исключаем Ti −1 с помощью (54): −1 Ti1 = (τ r + k ) C (Ti −01 + Ti +01 ) + ( k + τ r − 2C ) Ti 0 + f i1k 2  . (55) Формула (55) служит для определения значений температуры на первом временном слое. При описании процедуры расчета значений функции на верхнем (промежуточном и первом) временном слое применяем формулы (51) и (55). Для граничных условий второго рода производную в точке x0 ( x = 0 ) заменяют конечно-разностным аналогом со вторым порядком точности: T1 j − T−j1 = 0. 2h (56) Следовательно, значение температуры, записанное в фиктивном узле i = −1, равно: T−1j = T1 j . (57) Из предположения, что решение непрерывно записываем соотношение (51) для i =0: T0 j +1 = −τ rT0 j −1 + C (T−j1 + T1 j ) + ( k + 2τ r − 2C ) T0 j + f 0 j +1k 2 τr + k (58) , Из полученного соотношения исключаем T−j1 с помощью (57): T0 j +1 = 2CT1 j + ( k + 2τ r − 2C ) T0 j − τ rT0 j −1 + f 0 j +1k 2 , τr + k (59) Формула (59) служит для определения значений температуры на левой границе (в точке x0 ). Для первого слоя записываем соотношение (59) при j = 0 и исключаем T0−1 = T00 : T01 = 2CT10 + ( k + τ r − 2C ) T00 + f 01k 2 , τr + k (60) Для граничных условий в точке xN ( x = L ) производную заменяют конечноразностным аналогом со вторым порядком точности: TNj +1 − TNj −1 = 0. 2h (61) Следовательно, значение температуры, записанное в фиктивном узле i = N + 1, равно: TNj+1 = TNj −1 . (62) Записываем уравнение (51) для i = N : j +1 N T −τ rTNj −1 + C (TNj−1 + TNj +1 ) + ( k + 2τ r − 2C ) TNj + f Nj +1k 2 = τr + k , (63) Из полученного соотношения исключаем TNj+1 с помощью (62): j +1 N T 2CTNj−1 + ( k + 2τ r − 2C ) TNj − τ rTNj −1 + f Nj +1k 2 = , τr + k (64) Формула (64) служит для определения значений температуры на правой границе (в точке xN ). Для первого слоя записываем соотношение (64) при j = 0 и исключаем TN−1 = TN0 : 2CTN0−1 + ( k + τ r − 2C ) TN0 + f N1 k 2 T = , τr + k (65) 1 N Для решения уравнения методом прогонки необходимо построить сетку, покрывающую прямоугольник [0, L] ∪ [0,τ m ] . Координаты узлов сетки, образованные пересечением вертикальных и горизонтальных отрезков, определяются по формулам xi = ( i − 1) ⋅ h , τ j = j ⋅ k , где h − шаг по пространству, h = времени, k = τm M L , 1 ≤ i ≤ N ; k − шаг по N −1 , 0≤ j≤M . При получении трехслойной неявной разностной схемы подставляем формулы соответствующих разностных аналогов частных производных по времени и пространству в гиперболическое уравнение теплопроводности для каждого слоя. В этом случае разностное уравнение имеет вид: γ Ti j +1 − Ti j T j +1 − 2Ti j + Ti j −1 Ti +j1+1 − 2Ti j +1 + Ti −j1+1 λ + τ rγ i = + f i j +1 , k k2 h2 где f i j +1 = Aq0 (τ ) α e −α x + τ rα e−α x A ∂q0 (τ ) ∂τ (66) , i = 2, N , j = 0, M − 1 . Обе части уравнения (66) умножаем на k 2 h 2 : γ kh 2 (Ti j +1 − Ti j ) + τ r γ h 2 (Ti j +1 − 2Ti j + Ti j −1 ) = λ k 2 (Ti +j1+1 − 2Ti j +1 + Ti −j1+1 ) + f i j +1k 2 h 2 , Полученную систему приводим к системе линейных (67) алгебраических уравнений, имеющей трехдиагональную структуру и решаемую методом прогонки: Ai ⋅ Ti +j1+1 − Bi ⋅ Ti j +1 + Ci ⋅ Ti −j1+1 = Fi , (68) где Ai = Ci = λ k 2 , Bi = 2λ k 2 + γ kh 2 + τ r γ h 2 , Fi = − f i j +1k 2 h 2 − Ti j (γ kh 2 + 2τ r γ h 2 ) + τ r γ h 2 ⋅ Ti j −1 . Подставляя основное соотношение прогонки с уменьшенным на единицу нижним индексом Ti −j1+1 = α i −1 ⋅ Ti j +1 + β i −1 в уравнение (68), получаем формулы для расчета прогоночных коэффициентов: αi = Ai C ⋅β − F , β i = i i −1 i . Bi − Ci ⋅ α i −1 Bi − Ci ⋅ α i −1 (69) Для определения коэффициентов Fi на первом временном слое ( j = 0 ) вместо значения сеточной функции в фиктивном узле Ti −1 используем соотношение Ti −1 = Ti 0 , полученное из аппроксимации второго начального условия: Fi1 = − f i1k 2 h 2 − Ti 0 ( γ kh 2 + τ r γ h 2 ) . (70) Формулы для определения начальных прогоночных коэффициентов α1 , β1 находим из аппроксимации левого граничного условия T2j +1 − T0 j +1 = 0, T0 j +1 = T2 j +1 и 2h разностного уравнение (67) при i = 1 , получаем соотношение: T1 j +1 ( γ kh 2 + τ r γ h 2 + 2λ k 2 ) = 2λ k 2T2 j +1 + Ti j ( γ kh 2 + 2τ r γ h 2 ) − τ r γ h 2Ti j −1 + f1 j +1k 2 h 2 , (71) Выражая из равенства (71) T1 j +1 , получаем формулы для коэффициентов: α1j +1 = f1 j +1h 2 k 2 − τ r γ h 2T1 j −1 + T1 jγ h 2 ( k + 2τ r ) 2λ k 2 j +1 , β = . (τ r + k ) γ h 2 + 2λ k 2 1 (τ r + k ) γ h2 + 2λ k 2 (72) Определяя коэффициент β1 на первом временном слое, используем формулу (72) при условии j = 0 и значение сеточной функции в фиктивном узле T1−1 = T10 : f11h 2 k 2 + T10γ h 2 ( k + τ r ) β = . (τ r + k ) γ h2 + 2λ k 2 (73) 1 1 Значение температуры на правой аппроксимации правого граничного условия границе определяем, исходя из TNj ++11 − TNj−+11 = 0, TNj++11 = TNj −+11 и разностного 2h уравнение (67) при i = N , получаем соотношение TNj +1 ( γ kh 2 + τ r γ h 2 + 2λ k 2 ) = 2λ k 2TNj−+11 + TNj ( γ kh 2 + 2τ r γ h 2 ) − τ r γ h 2TNj −1 + f Nj +1k 2 h 2 . (74) Подставляя основное соотношение прогонки TNj−+11 = α N −1 ⋅ TNj +1 + β N −1 в уравнение (74), получаем формулы для расчета значениий температуры на правой границе: j +1 N T = 2λ k 2 β N −1 + TNj ( γ kh 2 + 2τ r γ h 2 ) − τ r γ h 2TNj −1 + f Nj +1k 2 h 2 γ kh 2 + τ r γ h 2 + 2λ k 2 (1 − α N −1 ) , (75) Определяя значение температуры в узле TN1 , используем формулу (75) при условии j = 0 и значение сеточной функции в фиктивном узле TN−1 = TN0 : 1 N T = 2λ k 2 β N −1 + TN0 ( γ kh 2 + τ r γ h 2 ) + f N1 k 2 h 2 γ kh 2 + τ r γ h 2 + 2λ k 2 (1 − α N −1 ) ,. (76) 8. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЭЛЛИПТПЧЕСКИХ ЗАДАЧ ПРИ ИССЛЕДОВАНИИ СТАЦИОНАРНЫХ ПРОЦЕССОВ 1. ПОСТАНОВКА КРАЕВЫХ ЭЛЛИПТИЧЕСКИХ ЗАДАЧ Процессы и явления, которые не изменяются с течением времени, в большинстве случаев описываются краевыми эллиптическими задачами для уравнений Лапласа и Пуассона. Например, в теории теплопроводности рассматривается задача о нахождении стационарного распределения температуры внутри области вокруг источника (или поглотителя) тепла, в электростатике задача о распределении электрического потенциала внутри области. Составление математической модели при описании стационарных процессов при постоянстве физических параметров приводит к решению уравнения Пуассона: (1.1) ∆u = f , где ∆u − оператор Лапласа. Для стационарных задач без наличия внутренних источников используется уравнение Лапласа: (1.2) ∆u = 0. Левая часть уравнения Лапласа или Пуассона задается в форме лапласиана ∆u имеющего в декартовых, цилиндрических и сферических координатах следующий вид: 1) ∆u = 2) ∂ 2u ∂ 2 u + (прямоугольная система координат, двумерное пространство); ∂x 2 ∂y 2 ∆u = ∂ 2u 1 ∂u 1 ∂ 2u + ⋅ + ⋅ ∂r 2 r ∂r r 2 ∂ϕ 2 (полярная система координат, двумерное (прямоугольная система координат, трехмерное пространство); 3) ∆u = ∂ 2u ∂ 2u ∂ 2 u + + ∂x 2 ∂y 2 ∂z 2 пространство); 4) ∆u = 1 1 ∂ 2u 1 ∂u 1 ∂ 2u ∂ 2u ′ ′′ ′′ ru + u + u или ∆ u = + ⋅ + ⋅ + (цилиндрическая ( r ) r 2 ϕϕ zz r r ∂r 2 r ∂r r 2 ∂ϕ 2 ∂z 2 система координат, трехмерное пространство); 1 2 1 1 ∂ 2u 2 ∂u ′ ′ 5) ∆u = 2 ( r ⋅ ur′ ) r + 2 ( sin θ ⋅ uθ′ ) θ + 2 2 uϕϕ′′ или ∆u = 2 + ⋅ + r r sin θ r sin θ ∂r r ∂r + 1 ∂ 2u ctgθ ∂u 1 ∂ 2u ⋅ + ⋅ + ⋅ r 2 ∂θ 2 r 2 ∂θ r 2 sin 2 θ ∂ϕ 2 (сферическая система координат, трехмерное пространство). Для одномерного случая решение задачи сводится к решению обыкновенного дифференциального уравнения, не вызывая затруднений. Среди двумерных задач рассматривается три основных вида: плоские (решаются в двумерной прямоугольной системе координат); плоскопараллельные (предполагается, что геометрия расчетной области, свойства сред и параметры, характеризующие источники поля, неизменны в направлении оси OZ); осесимметричные (решаются в цилиндрической системе координат, физические свойства и источники поля предполагаются не зависящими от угловой координаты). Для получения полного математического описания определенного стационарного процесса к условиям однозначности наряду с физическими и геометрическими добавляют граничные условия, определяющие особенности протекания процесса на поверхности тела. В отличие от нестационарного случая постановка краевых задач для стационарного уравнения не требует наличие начальных условий. Выделяют следующие виды краевых задач: 1) задача Дирихле (задача с граничными условиями первого рода); 2) задача Неймана (задача с граничными условиями второго рода); 3) задача Робэна (задача с граничными условиями третьего рода); 4) задача с условиями сопряжения (задача с граничными условиями четвертого рода). 2. ТРЕХСЛОЙНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ УРАВНЕНИЯ ПУАССОНА Рассмотрим применение метода сеток с использованием конечно-разностных аппроксимаций (метод конечных разностей) на примере уравнения Пуассона. Заменяем двумерную область сеточной областью, при этом шаг по обеим координатам берется равный. Для аппроксимации каждой из частных производных центральной разностной производной по соответствующей координате, можно использовать следующие формулы: ∂ 2u ui −1, j − 2ui , j + ui +1, j = + O ( h2 ) , 2 2 h ∂x (2.1) ∂ 2u ui , j −1 − 2ui , j + ui , j +1 = + O ( h2 ) , ∂y 2 h2 (2.2) где ui , j = u ( xi , y j ) − значение сеточной функции в точке ( i, j ) , O ( h 2 ) – погрешность аппроксимации (невязка), h − шаг по координатам. Подставляя приведенные разностные аналоги производных по пространственным переменным в соотношение (1.1), получим сеточное (разностное) уравнение: ui −1, j − 2ui , j + ui +1, j h 2 + ui , j −1 − 2ui , j + ui , j +1 h2 = fi, j . (2.3) Выражая ui , j , получим формулу для вычисления значений сеточной функции во внутренних узлах сетки: ui , j = ui −1, j + ui +1, j + ui , j −1 + ui , j +1 − Si , j 4 (2.4) . где Si , j = fi , j ⋅ h 2 . Значение сеточной функции в точке i, j+1 ( i, j ) находится как среднее арифметическое между значениями сеточной функции в четырех соседних узлах. i-1, j Трехслойная разностная схема такого типа имеет название «крест», шаблон предложенной схемы представлен на рис. 1. i, j i+1, j i, j-1 Рис. 1 Если заданы граничные условия по всему периметру расчетной области, то для всех внутренних точек имеем систему алгебраических уравнений, матрица коэффициентов которой содержит ненулевые элементы на пяти диагоналях (главной, двух верхних и нижних). К недостатку реализации данной схемы относится тот факт, что количество уравнений в системе определяется величиной выбранного шага. Уменьшение величины шага приводит к увеличению количества точек решетки, для каждой из которых выписывается соответствующее уравнение. Дальнейшее решение можно вести по двум направлениям: 1) модифицировать для пятидиагональной матрицы алгоритм прогонки, предложенный для матриц с тремя диагоналями; 2) использовать итерационные методы. Для реализации итерационного метода релаксации в правой части уравнения (2.4) прибавляют и вычитают ui , j , затем представляют разностное уравнение в виде: ui , j = ui , j + ui −1, j + ui +1, j + ui , j −1 + ui , j +1 − 4ui , j − Si , j 4 (2.5) Из формулы (2.5) следует, что итерационный метод релаксации основывается на последовательном вычислении для каждого внутреннего узла нового значения исходя из полученного на предыдущей итерации. Обязательным условием является предварительное задание начального приближения для всех внутренних точек с использованием граничных условий. Большое количество узлов сетки приводит к крайне медленной сходимости решения. Для ускорения сходимости в формулу (2.5) вводится дополнительный множитель ω ∈ (0; 2) : ui , j = ui , j + ω ui −1, j + ui +1, j + ui , j −1 + ui , j +1 − 4ui , j − Si , j 4 (2.6) Формула (2.6) последовательно применяется к узлам сетки на каждой итерации. Вычисления останавливаются, если для всех точек разность между значениями, полученными на данной и предыдущей итерации, окажется ниже заданной величины ошибки. Рассматривая в формуле (2.4) правую часть на предыдущей итерации, а левую – на текущей, получим каноническую форму записи сеточного метода: u k +1 i, j = uik−1, j + uik+1, j + uik, j −1 + uik, j +1 − Si , j 4 , (2.7) где uik, j − значение сеточной функции в точке ( i, j ) , k − номер итерации. Формула (2.7) используется в методе Якоби для нахождения следующего приближения к решению. Разностное уравнение (2.3) при реализации метода Зейделя, учитывающего результаты вычисления на k + 1 итерации, записывается в виде: uik−+1,1 j − 2uik, +j 1 + uik+1, j h2 + uik, +j 1−1 − 2uik, +j 1 + uik, j +1 h2 = fi, j . (2.8) Выражая uik, +j 1 , получим каноническую формулу для метода Зейделя: k +1 i, j u = uik−+1,1 j + uik+1, j + uik, +j −11 + uik, j +1 − Si , j 4 . (2.9) Дальнейшее обобщение метода Зейделя с использованием формулы (2.6) и введением в правую часть параметра релаксации 1 < ω < 2 , приводит к записи формулы для метода верхней релаксации: k +1 i, j u k i, j = (1 − ω )u + ω uik−+1,1 j + uik+1, j + uik, +j −11 + uik, j +1 − Si , j (2.10) 4 При ω = 1 формула (2.10) приводится к канонической формуле для метода Зейделя. Используя в соотношении (2.10) 0 < ω <1 получаем формулу для метода последовательной нижней релаксации. В случае неравных координатных шагов разностное уравнение (2.8) можно записать в виде: uik−+1,1 j − 2uik, +j 1 + uik+1, j hx2 + uik, +j −11 − 2uik, +j 1 + uik, j +1 hy2 (2.11) = fi , j . Значения сеточной функции в этом случае находятся по формуле: u k +1 i, j = Fi − Ai uik−+1,1 j − Bi uik, +j −11 Ci (2.12) , uik+1, j uik, j +1 1 2 2 1 где Ai = 2 , Bi = 2 , Ci = − 2 − 2 , Fi = fi , j − 2 − 2 . hx hy hx hy hx hy Правильная реализация вычислительного процесса позволяет получить в каждом последующем приближении более точный результат по сравнению с предыдущим. Итерационный процесс, при котором погрешность с каждой итерацией уменьшается, определенность называют матрицы сходящимся. коэффициентов Симметричность соотношения и (2.7) положительная обосновывает сходимость методов Якоби и Зейделя при любом выборе начального приближения. Среди всех рассмотренных итерационных процедур метод последовательной верхней релаксации сходится наиболее быстро. При этом в рамках самого метода наибольшая скорость сходимости достигается при некотором оптимальном значении параметра релаксации ω , которое зависит от элементов матрицы коэффициентов в уравнении (2.10). Считается, что при ω = 1, 5 достигается достаточно быстрая сходимость, что приводит к эффективной численной реализации метода для широкого круга проблем. 3. РЕШЕНИЕ ДВУМЕРНОЙ СТАЦИОНАРНОЙ ЗАДАЧИ С ИСПОЛЬЗОВАНИЕМ ВСТРОЕННЫХ ФУНКЦИЙ СИСТЕМЫ MATHCAD В системе MathCAD для решения эллиптических уравнений вида: A +L ∂ 2 u ( x, y ) ∂ 2 u ( x, y ) ∂ 2 u ( x, y ) + 2 B + C + ∂x 2 ∂x∂y ∂y 2 ∂u ( x, y ) ∂u ( x, y ) +M + N ⋅ u ( x, y ) = F ( x, y ), ∂x ∂y (3.1) где B 2 − AC < 0, предназначена функция relax, реализующая итерационный метод релаксации. Используя аппроксимации частных производных, уравнение (3.1) можно привести к следующему разностному уравнению: (3.2) ai , j uij+1 + bi , j uij−1 + ci , j uij +1 + di , j uij −1 + ei , j uij = fi , j h2 Для корректной работы функции relax необходимо чтобы по обеим координатам рассматривался равный шаг h и одинаковые промежутки определения i и j , в этом случае расчетная область представляет собой квадрат. В частности, для уравнения Пуассона и уравнения Лапласа в соответствии с формулой (2.3) коэффициенты ai , j = bi , j = ci , j = di , j = 1, ei , j = −4. Функция relax(a, b, c, d, e, f, v, rjac) возвращает квадратную матрицу значений сеточной функции, в случае, когда известны значения искомой функции на всех четырех сторонах квадратной области. В перечне аргументов функции relax указывают: 1) a, b, c, d, e – квадратные матрицы одного и того же размера, содержащие коэффициенты аппроксимирующей дифференциальное уравнение пятиточечной разностной схемы; 2) f – квадратная матрица, содержащая значения правой части уравнения в каждой точке внутри квадрата; 3) v – квадратная матрица, содержащая граничные значения функции на краях области, а также начальное приближение решения во внутренних точках области; 4) rjac – параметр, управляющий сходимостью процесса релаксации, изменяющийся в диапазоне от 0 до 1. С помощью встроенной функции relax можно решать задачи нестационарной теплопроводности, если разносная схема имеет шаблон, соответствующий рис. 1. В случае, когда по всему периметру расчетной области заданы нулевые граничные условия, можно использовать специальную функцию multigrid(F,ncycle). В перечне аргументов функции указывают: 1) F – матрица размерности ( M + 1) × ( M + 1) , определяющая правую часть уравнения Пуассона. При использовании этой функции M = 2n ( n ∈ N ) и сторона квадрата обязательно содержит 2n + 1 узлов. 2) ncycle – параметр численного алгоритма, определяющий количество циклов, используемых на каждой итерации (для большинства случаев используют значение 2). 4. ТЕМПЕРАТУРНОЕ ПОЛЕ В ПЛАСТИНЕ С НЕРАВНОМЕРНЫМ ПОДВОДОМ И ОТВОДОМ ТЕПЛОТЫ НА ГРАНИЦЕ Достаточно часто на практике поверхность теплообмена охлаждается (нагревается) неравномерно, в этом случае плотность теплового потока непостоянна и меняется вдоль поверхности теплообмена (например, экранные трубы в топке котла, обращенные наполовину к факелу и продуктам сгорания, а наполовину к стенке топки). Рассматривая математическую модель расчета температурного поля в пластине с неравномерным подводом и отводом теплоты на границе, получаем задачу о стационарном температурном поле в квадратной области с граничными условиями третьего рода. Численное решение задачи в этом случае можно получить на основе итерационного метода Зейделя. Математическая формулировка задачи имеет вид: Найти в области D : 0 ≤ x ≤ L, 0 ≤ y ≤ L решение Т ( x, y ) уравнения стационарной теплопроводности λ λ ∂T ∂y ∂T ∂x q ∂ 2T ∂ 2T + 2 = − v , удовлетворяющее граничные условия: 2 ∂x ∂y λ = α1 (T ( x, 0 ) − Tg1 ) , −λ y =0 x=0 = α 3 (T ( 0, y ) − Tg3 ) , −λ ∂T ∂y y=L ∂T ∂x x=L = α 2 (T ( x, L ) − Tg 2 ) , x ∈ [ 0, L ] ; (4.1) = α 4 (T ( L, y ) − Tg 4 ) , y ∈ [ 0, L ] , (4.2) где Tg i − температура среды по соответствующей стороне пластины. При реализации метода Зейделя для вычисления значений сеточной функции во внутренних узлах сетки используется формула (2.9). Для аппроксимации граничных условий третьего рода производную на границах y = 0 и x = 0 заменяют правой разностной производной, а на границах y = L и x = L − левой разностной производной: Ti ,1 − Ti ,0 h Ti ,n − Ti ,n −1 h Вводя сеточное T1, j − T0, j α 3 α1 Ti ,0 − Tg1 ) , = (T0, j − Tg 3 ) , ( λ λ h (4.3) Tn , j − Tn −1, j α2 α Ti ,n − Tg 2 ) , = − 4 (Tn , j − Tg 4 ) ( λ λ h (4.4) = =− число Био Bik = αk ⋅ h λ и выражая Ti ,0 , Ti ,n , T0, j , Tn , j в соотношениях (4.3) и (4.4), получим формулы для определения значений температуры в узлах на границе области через температуру окружающей среды и температуру ближайшего внутреннего узла: Ti ,0 = T0, j = Bi1 ⋅ Tg1 + Ti ,1 Bi ⋅ Tg 2 + Ti ,n −1 , Ti ,n = 2 , 1 + Bi1 1 + Bi2 Bi3 ⋅ Tg 3 + T1, j 1 + Bi3 , Tn , j = Bi4 ⋅ Tg 4 + Tn −1, j 1 + Bi4 (4.5) . (4.6) Для получения численного решения поставленной задачи в системе MathCAD с реализацией метода сеток используем трехслойную неявную схему вида «крест». При описании процедуры расчета значений сеточной функции применяем формулы для метода Зейделя (2.9), (4.5) и (4.6). Алгоритм решения поставленной задачи основан на применении итерационного метода Зейделя. Особенностью реализации данного метода является использование при вычислении значений сеточной функции на k + 1 итерации уже найденных значений на этой же итерации в узлах с меньшими номерами. Процедура расчета температурного поля пластины включает цикл по переменной iter ∈1 .. maxiter . Максимальное число итераций и точность tol задаются в процедуре. Расчет температурного поля с учетом формул для нахождения значения сеточной функции во внутренних узлах и на границе области реализуется через внутренние циклы по переменным i, j. Вычисления в цикле прекращаются, если максимальная разность в двух последних итерациях принимает значение меньше заданной точности. 5. РАСЧЕТ ПРОСТРАНСТВЕННОГО РАСПРЕДЕЛЕНИЯ ПОТЕНЦИАЛА ЭЛЕКТРИЧЕСКОГО ПОЛЯ В ПОЛОСКОВОЙ ЛИНИИ В настоящее время в технике СВЧ широко применяется особый класс линий передачи с волнами типа Т, называемых полосковыми волноводами или полосковыми линиями. В этих волноводах токонесущие проводники представляют собой тонкие полоски металла, между которыми находится подложка – плоский слой диэлектрика с малыми потерями. Полосковые линии бывают симметричными и несимметричными; поперечные сечения их изображены на рис. 2. В Рис. 2. Полосковые линии передачи: а) симметричная, б) несимметричная технической литературе несимметричные полосковые линии для сантиметрового и миллиметрового диапазонов часто называют микрополосковыми линиями, подчеркивая этим термином миниатюрность конструкций. Строгий электродинамический анализ полей в несимметричной полосковой линии является достаточно сложной задачей и проводится в основном численными методами. Это связано с тем, что параметры заполняющей среды в полосковой линии неоднородны по сечению. На практике обычно применяют линии, у которых толщина подложки существенно меньше ширины верхнего проводника. Поэтому электрическое поле в поперечном сечении линии распределено примерно так же, как и электростатическое поле в плоском конденсаторе. Достаточно высокое значение относительной диэлектрической проницаемости подложки снижает роль краевых эффектов, так что поле во внутренней области оказывается приблизительно однородным. Следовательно, в этом случае можно обоснованно пренебречь сравнительно малыми продольными проекциями электрического и магнитного полей. Конструкция полосковой линии (ПЛ) включает: металлический проводник (полоска) шириной w и толщиной t лежит на обеспечивающей жесткость конструкции подложке толщиной h, выполненной из однородного диэлектрика с относительной диэлектрической проницаемостью ε. В качестве диэлектрика чаще всего используется керамика на основе окиси алюминия с ε = 9, 6 ÷ 9,8 (поликор, сапфирит и др.). С внешней стороны подложка покрыта слоем металла (экраном). Поперечное сечение ПЛ приведено на рис. 3., где 1 – полоска, 2 – подложка, 3 – экран. Рис.3. Поперечное сечение полосковой линии При многослойном заполнении линия передачи содержит несколько сред с постоянными диэлектрическими проницаемостями. В пределах каждого из однородных диэлектриков распределение потенциала описывается уравнением Лапласа (дифференциальное уравнение в частных производных второго порядка относительно скалярного электрического потенциала) для двух- или трехмерного случая: ∂ 2ϕ ∂ 2ϕ + =0 ; ∂x 2 ∂y 2 (5.1) ∂ 2ϕ ∂ 2ϕ ∂ 2ϕ + + = 0. ∂x 2 ∂y 2 ∂z 2 (5.2) Для нахождения погонных емкостей линий передачи СВЧ используется краевая задача нахождения решения уравнения Лапласа для двух- или трехмерного потенциала, удовлетворяющего заданным граничным условиям. Затем, когда распределение потенциала найдено, по нему находится структура электрического поля и емкость линии. Далее определяется эффективная диэлектрическая проницаемость и волновое сопротивление линии. Для построения конечно-разностной схемы в приведенное выше двухмерное уравнение Лапласа подставляются аппроксимации вторых производных с учетом равномерности сетки и получают разностное уравнение вида: ϕi −1, j + ϕi +1, j + ϕi , j −1 + ϕi , j +1 − 4ϕi , j = 0 . (5.3) Для потенциалов точек, лежащих на границах раздела подобластей (рис. 4), необходимо использовать уравнение связи узловых потенциалов. В частности, если узловые точки расположены на параллельной оси x границе двух сред с диэлектрическими проницаемостями ε1 и ε 3 , а сетка является равномерной, то ε 4 = ε 3 , ε 2 = ε1 . В этом случае уравнение связи узловых потенциалов имеет вид: ε13ϕi −1, j + ε13ϕi +1, j + ε1ϕi , j −1 + ε 3ϕi , j +1 − 4ε13ϕi , j = 0 , где ε13 = ε1 + ε 3 2 (5.4) . Аналогично, для параллельной оси y границы раздела диэлектриков с проницаемостями ε 2 и ε 4 имеем: ε 4ϕi −1, j + ε 2ϕi +1, j + ε 24ϕi , j −1 + ε 24ϕi , j +1 − 4ε 24ϕi , j = 0 , где ε 24 = ε2 + ε4 2 (5.5) . Рис. 4. Сетка с узловыми точками, лежащими на границах диэлектриков В итерационную формулу метода последовательной релаксации вносится соответствующая корректировка для точек, принадлежащих границам раздела слоев диэлектрика. Например, исходя из конечно-разностного уравнения (2.10), нетрудно заключить, что итерации следует проводить по формуле ϕip, j+1 = (1 − ω ) ϕip, j + ω (ϕ 4 p +1 i −1, j + ϕip+1, j + ϕip, j+−11 + ϕip, j +1 ) , (5.6) Потенциалы лежащих на границе раздела двух диэлектриков узловых точек методом последовательной верхней релаксации уточняются по формуле: ϕip, j+1 = (1 − ω ) ϕip, j + ω ε13ϕip−+1,1j + ε13ϕip+1, j + ε 1ϕip, j+−11 + ε 3ϕip, j +1 ) , ( 4ε 13 (5.7) Для получения численного решения поставленной задачи в системе MathCAD с реализацией метода сеток используем трехслойную неявную схему вида «крест». При описании процедуры расчета значений сеточной функции применяем формулы для метода последовательной верхней релаксации (5.6), (5.7) в случае равномерной сетки. Процедура расчета потенциала электрического поля включает в себя цикл по переменной m . Максимальное число итераций задается в процедуре. Расчет поля с учетом формул для нахождения значения сеточной функции во внутренних узлах области реализуется через внутренние циклы по переменным i, j. Вычисления в цикле прекращаются, если максимальная разность в двух последних итерациях принимает значение меньше заданной точности δ . 9. ОПТИМИЗАЦИОННЫЕ МОДЕЛИ 1. ОСНОВНЫЕ ПОНЯТИЯ Рис.1 Рис.2 Содержание математического программирования составляют теория и методы решения задач, в математических моделях которых требуется определить экстремум некоторой функции нескольких переменных, связанных условиями в виде равенств и неравенств:  f ( x) → max (min);  g ( x) ≤ 0, i = 1, … , k ;  i   gi ( x) = 0, i = k + 1, … , m;  x ∈ Rn ,  где f ( x), g i ( x) (i = 1, … , m) − некоторые функции n переменных. (1) Если функции f ( x), g i ( x) (i = 1, … , m) линейные, то задача (1) называется задачей линейного программирования (ЛП). Отдельный класс составляют так называемые задачи ЛП транспортного типа и задачи сетевого планирования, для решения которых разработаны методы, учитывающие их специфику. 2. ЗАДАЧА ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ Определение 1. Задачей линейного программирования (ЛП) называется задача математического программирования вида: n f ( x) = ∑ c j x j → max (min); (1) j =1 n gi ( x) = ∑ aij x j ≤ bi , i = 1, … , k ; (2) j =1 n gi ( x) = ∑ aij x j ≥ bi , i = k + 1, … , l; (3) j =1 n gi ( x) = ∑ aij x j = bi , i = l + 1, … , m; (4) j =1 x j ≥ 0 для некоторых j ∈ { 1, … , n } . (5) Определение 2. Функция (1) называется целевой. Определение 3. Любой вектор x ∈ R n , удовлетворяющий условиям (2)−(5), называется допустимым решением, или планом задачи ЛП. Множество всех планов задачи обозначим D. Определение 4. Вектор x∗ ∈ D называется оптимальным решением (оптимальным планом) задачи ЛП, если f ( x∗ ) ≥ f ( x) для любого x ∈ D в задаче максимизации или f ( x∗ ) ≤ f ( x) для любого x ∈ D в задаче минимизации. Задача ЛП (1)−(5) называется разрешимой, если она имеет хотя бы один оптимальный план. Из определения задачи ЛП ясно, что она относится к классу задач на условный экстремум, однако частные производные целевой функции постоянны, т. е. не зависят от точек ∂f = c j ( j = 1, … , n) ∂x j . Следовательно, обычные методы математического анализа неприменимы для решения задач ЛП, поэтому для таких задач созданы специальные методы, учитывающие их специфику. Особый класс задач составляют так называемые канонические задачи. Определение 5. Канонической задачей ЛП (КЗЛП) называется задача вида: n  f ( x ) c j x j → max; = ∑  j =1   n  ∑ aij x j = bi , i = 1, … , m;  j =1  x ≥ 0, j = 1, … , n. j   (6) Заметим, что характерными особенностями КЗЛП являются следующие: целевая функция максимизируется; из ограничений присутствуют только ограничения, имеющие вид равенств; все переменные задачи неотрицательны. Теорема 1 (теорема об эквивалентности). Для всякой задачи ЛП существует эквивалентная ей каноническая задача ЛП. Правила приведения ЗЛП к каноническому виду: 1. Если в исходной задаче некоторое ограничение (например, первое) было неравенством, то оно преобразуется в равенство, введением в левую часть некоторой неотрицательной переменной, причем в неравенства «≤» вводится дополнительная неотрицательная переменная со знаком «+»; в случаи неравенства «≥» - со знаком «-» В каждое из неравенств вводится своя “уравнивающая” переменная, после чего система ограничений становится системой уравнений. 2. Если в исходной задаче некоторая переменная не подчинена условию неотрицательности, то ее заменяют (в целевой функции и во всех ограничениях) разностью неотрицательных переменных 3. Если в ограничениях правая часть отрицательна, то следует умножить это ограничение на (-1). 4. Если исходная задача была задачей на минимум, то введением новой целевой функции F1 = -F мы преобразуем нашу задачу на минимум функции F в задачу на максимум функции F1. Таким образом, всякую задачу сформулировать в канонической форме. линейного программирования можно В стандартной же форме задача линейного программирования является задачей на максимум (минимум) линейной целевой функции. Система ограничений ее состоит из одних линейных неравенств типа «<=» или «>=». Все переменные задачи неотрицательны. 3. ГРАФИЧЕСКИЙ МЕТОД РЕШЕНИЯ ЗАДАЧИ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ Задача ЛП, имеющая только две переменные, может быть решена графически. В этом случае допустимое множество D задачи ЛП (если оно не пусто) образует на плоскости некоторую область, границами которой являются прямые линии. Такое множество назовем многогранным. Задача ЛП сводится к поиску точки многогранного множества D, в которой целевая функция f ( x) = c1 x1 + c2 x2 принимает наибольшее (наименьшее) значение. Напомним, что линией уровня функции f(x) называется линия, во всех точках которой f(x) принимает одно и то же значение C (C ∈ R) , т. е. линия, заданная уравнением f ( x) = C . Градиентом функции f(x) называется вектор grad f = где ∂f ∂f , ∂ x1 ∂ x2 ∂f ∂f i + j, ∂ x1 ∂ x2 − частные производные функции f(x). Вектор (1) − grad f , т. е. противоположный вектору grad f , называется антиградиентом. Если f(x) − линейная функция, то grad f = (c1 , c2 ) − постоянный вектор, а линии уровня − прямые линии c1 x1 + c2 x2 = C , перпендикулярные градиенту. Из курса высшей математики известно, что градиент указывает направление наискорейшего возрастания функции f(x), а антиградиент − убывания, поэтому, чтобы найти наибольшее (наименьшее) значение целевой функции f(x), необходимо мысленно сдвигать линию уровня этой функции в направлении градиента (антиградиента). Точка, в которой линия уровня функции в последний раз касается множества D, будет оптимальной. З а м е ч а н и я. 1) Если линия уровня целевой функции f(x) касается множества D в последний раз по целому ребру, то задача ЛП имеет бесконечное множество оптимальных решений, они соответствуют каждой точке этого ребра. 2) Если допустимое множество D таково, что сколько бы ни сдвигали линию уровня целевой функции f(x) в направлении градиента для задачи на максимум (антиградиента − для задачи на минимум), она всегда будет пересекать множество D, то функция неограниченно возрастает (убывает) на множестве D. Это означает, что задача ЛП не имеет решения. Алгоритм графического метода решения задач ЛП: 1. Построить область допустимых решений. 2. Если область допустимых решений является пустым множеством, то задача не имеет решения ввиду несовместности системы ограничений. 3. Если область допустимых решений является непустым множеством, построить нормаль линий уровня n = (c1 ,c2 ) и одну из линий уровня, имеющую общие точки с этой областью. 4. Линию уровня переместить до опорной прямой в задаче на максимум в направлении нормали, в задаче на минимум − в противоположном направлении. 5. Если при перемещении линии уровня по области допустимых решений в направлении, соответствующем приближению к экстремуму целевой функции, линия уровня уходит в бесконечность, то задача не имеет решения в виду неограниченности целевой функции. 6. Если задача линейного программирования имеет оптимальное решение, то для его нахождения решить совместно уравнения прямых, ограничивающих область допустимых решений и имеющих общие точки c соответствующей опорной прямой. Если целевая функция задачи достигает экстремума в двух угловых точках, то задача имеет бесконечное множество решений. Оптимальным решением является любая выпуклая линейная комбинация этих точек. После нахождения оптимальных решений вычислить значение целевой функции на этих решениях. 10. СИМПЛЕКС-МЕТОД РЕШЕНИЯ ЗАДАЧ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ 1. СПЕЦИАЛЬНАЯ ЗАДАЧА ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ Рассмотрим систему линейных уравнений, входящую в ограничения канонической задачи ЛП         n f ( x) = ∑ c j x j = (c, x) → max; j =1 n ∑A x j j = B; j =1 (1) x j ≥ 0, j = 1, … , n, Матрицу коэффициентов левых частей уравнений этой системы обозначим А. Определение 1. Система линейных уравнений называется системой с базисом, если в каждом уравнении имеется переменная, которая входит в это уравнение с коэффициентом 1 и отсутствует в остальных уравнениях. Такие переменные называются базисными, остальные − небазисными. Матрицу А системы с базисом с помощью перенумерации переменных можно привести к виду:  1 0 ⋯ 0 a1, m +1  0 1 ⋯ 0 a2, m +1 A= ⋯ ⋯ ⋯ ⋯ ⋯   0 0 ⋯ 1 am, m +1 a1, m + 2 a2, m + 2 ⋯ am , m + 2 ⋯ a1n   ⋯ a2 n  . ⋯ ⋯  ⋯ am n  (2) Определение 2. Каноническая задача ЛП называется специальной, если 1) ее система линейных уравнений является системой с базисом; 2) правые части уравнений системы неотрицательны: bi ≥ 0, i = 1, 2, … , m ; 3) целевая функция f ( x) выражена только через небазисные переменные. Специальная задача ЛП имеет вид:  f ( x) = c0 − cm +1 xm +1 − … − cn xn → max;  x + a1, m +1 xm +1 + … + a1n xn = b1 ;  1 x2 + a2, m +1 xm +1 + … + a2 n xn = b2 ;   ………………………………………………  xm + am, m +1 xm +1 + … + am n xn = bm ,  где x1 , x2 , … , xm bi ≥ 0, i = 1, 2, … , m . − базисные, xm +1 , xm + 2 , … , xn − небазисные (3) переменные, 2. СИМПЛЕКС-МЕТОД Идея симплекс-метода состоит в переборе базисных решений специальной задачи ЛП, число которых конечно, но может быть настолько большим, что перебор всех базисных решений потребует слишком много времени даже при применении современной вычислительной техники. Кроме того, теоретически существует возможность зацикливания симплекс-метода, когда проходятся по кругу одни и те же базисные решения, не достигая оптимального. Однако несмотря на эти отрицательные моменты, симплекс-метод за годы своего существования доказал свою практическую пригодность. Дело в том, что вероятность зацикливания симплексметода ничтожно мала. К тому же появились модификации симплекс-метода, избавленные от этого недостатка. Описание симплекс-метода решения задачи ЛП дадим в терминах симплексных таблиц. Сначала изложим метод для специальной задачи ЛП, а затем покажем, как каноническую задачу ЛП привести к виду специальной. Определение 3. Симплексной таблицей специальной задачи ЛП (3) называется таблица вида табл. 1. Строка симплексной таблицы, соответствующая целевой функции f ( x) , называется индексной. Определение 4. Симплексная таблица называется допустимой, если все элементы столбца x0 , за исключением быть может c0 , неотрицательны. Определение 5. Две симплексные таблицы называются эквивалентными, если соответствующие им специальные задачи ЛП эквивалентны. Заметим, что вектор x 0 = (b1 , b2 , … , bm , 0, … , 0) является базисным решением задачи (3), причем f ( x 0 ) = c0 . Действительно, если bi = 0, i = 1, 2, … , m, то вектор x0 = 0 ; если не все bi равны нулю, то система векторов-столбцов Aj , соответствующих ненулевым компонентам x j 0 вектора x 0 , линейно независима, как подсистема системы столбцов единичной матрицы Em . Будем говорить, что базисное решение x 0 соответствует симплексной таблице. Таблица 1 Перейдем к изложению симплекс-метода. Шаг 0 (начальный шаг). Для специальной задачи (3) составляем начальную симплексную таблицу. Соответствующее ей базисное решение x 0 = (b1 , b2 , … , bm , 0, … , 0) назовем текущим. Переходим к шагу 1. Шаг 1 (проверка на оптимальность). Если c j ≥ 0 для всех j = 1, 2, … , n , то конец работы метода − текущее базисное решение x 0 оптимально. В противном случае переходим к шагу 2. Шаг 2 (проверка на разрешимость). Если существует такой номер q, что сq < 0 , и среди элементов q-го столбца нет положительных: aiq ≤ 0, i = 1, 2, … , m, то конец работы метода − задача неразрешима. В противном случае переходим к шагу 3. Шаг 3 (выбор ведущего столбца). Столбец, имеющий номер q, назовем ведущим, если сq < 0 и aiq > 0 для некоторого i ∈ { 1, 2, … , m } . Если таких столбцов несколько, то выбираем любой из них. Переходим к шагу 4. Шаг 4 (выбор ведущей строки). Строку номер p назовем ведущей, если bp a pq  b  = min  i : aiq > 0  .  aiq  (4) Если таких строк несколько, то выбираем любую из них. Переходим к шагу 5. ′ Шаг 5 (преобразование симплексной таблицы). Элементы aij , bi′ , c j′ новой симплексной таблицы вычисляем по следующим формулам: b p′ = bi′ = bi − bp a pq a pq a pq ; a pj′ = aiq ; aij′ = aij − c0′ = c0 − Элемент bp bp a pq a pj a pq a pj a pq , где j = 1, 2, … , n; (5) aiq , где i ≠ p; j = 1, 2, … , n; cq ; c j′ = c j − a pj a pq cq , где j = 1, 2, … , n. называется ведущим, а вычисления по формулам (5) − преобразованиями с ведущим элементом a pq . В столбце базисных переменных В новой симплексной таблицы заменяем x p на xq . Возвращаемся к шагу 1. Далее вся процедура повторяется до тех пор, пока либо на шаге 1 не обнаружится, что текущее базисное решение оптимально, либо на шаге 2 не обнаружится, что задача неразрешима. За м е ч а н и я. 1) Формулы (5) задают элементарные преобразования системы линейных ограничений специальной задачи (3) (как при решении системы линейных уравнений методом Гаусса), которые, как известно, не меняют множества решений исходной системы. Следовательно, новая симплексная таблица, элементы которой найдены по формулам (5), эквивалентна исходной таблице. 2) Правило выбора ведущей строки (4) гарантирует, что решение, соответствующее новой симплексной таблице, является допустимым и базисным. Следующие три утверждения дают теоретическое обоснование симплексметода. Теорема 1 (условия оптимальности решения). Если cj ≥ 0 для всех j = 1, 2, … , n , то текущее базисное решение задачи (3) оптимально. Теорема 2 (условие неразрешимости задачи). Если существует такой номер q, что сq < 0 и среди элементов q-го столбца нет положительных: aiq ≤ 0, i = 1, 2, … , m, то целевая функция f(x) неограниченно возрастает на множестве D допустимых решений задачи (3). Теорема 3 (о переходе к новому базисному решению). Если существует такой номер q, что сq < 0 и для некоторого i ∈ {1, 2, … , m} , то преобразование симплексной таблицы по формулам (19) приводит к получению допустимого базисного решения, причем f ( x′) ≥ f ( x 0 ) , где − x′, x 0 базисные решения, соответствующие старой и новой таблицам. При решении задач можно использовать следующий алгоритм перехода от таблицы к таблице: • просматривается последняя строка (индексная) таблицы и среди коэффициентов этой строки (исключая столбец свободных членов) выбирается наименьшее отрицательное число при отыскании max, либо наибольшее положительное при задачи на min. Если такового нет, то исходное базисное решение является оптимальным и данная таблица является последней; • просматривается столбец таблицы, отвечающий выбранному отрицательному (положительному) коэффициенту в последней строке - ключевой столбец, и в этом столбце выбираются положительные коэффициенты. Если таковых нет, то целевая функция неограниченна на области допустимых значений переменных и задача решений не имеет; • среди выбранных коэффициентов столбца выбирается тот, для которого абсолютная величина отношения соответствующего свободного члена (находящегося в столбце свободных членов) к этому элементу минимальна. Этот коэффициент называется разрешающим (ведущим), а строка в которой он находится ключевой (ведущей); • в дальнейшем базисная переменная, отвечающая строке разрешающего элемента, должна быть переведена в разряд свободных, а свободная переменная, отвечающая столбцу разрешающего элемента, вводится в число базисных. Строится новая таблица, содержащая новые названия базисных переменных: • разделим каждый элемент ключевой строки (исключая столбец свободных членов) на разрешающий элемент и полученные значения запишем в строку с измененной базисной переменной новой симплекс таблицы. • строка разрешающего элемента делится на этот элемент и полученная строка записывается в новую таблицу на то же место. • в новой таблице все элементы ключевого столбца = 0, кроме разрезающего, он всегда равен 1. • таким же. столбец, у которого в ключевой строке имеется 0,в новой таблице будет • строка, у которой в ключевом столбце имеется 0,в новой таблице будет такой же. • в остальные клетки новой таблицы записывается результат преобразования элементов старой таблицы: В результате получают новую симплекс-таблицу, отвечающую новому базисному решению. Теперь следует просмотреть строку целевой функции (индексную), если в ней нет отрицательных значений (в задачи на нахождение максимального значения), либо положительных (в задачи на нахождение минимального значения) кроме стоящего на месте свободного столбца, то значит, что оптимальное решение получено. В противном случае, переходим к новой симплекс таблице по выше описанному алгоритму. 11. НАХОЖДЕНИЕ НАЧАЛЬНОГО БАЗИСНОГО РЕШЕНИЯ МЕТОДОМ ИСКУССТВЕННОГО БАЗИСА Рассмотрим каноническую задачу ЛП:  f ( x) = c1 x1 + c2 x2 + … + cn xn → max;  a11 x1 + a12 x2 + … + a1n xn = b1 ;   ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯  a x + a x +… + a x = b ; m1 1 m2 2 mn n m  x j ≥ 0, j = 1, 2, … , n.  Можно считать, что все bi ≥ 0, i = 1, 2, … , m (1) (в противном случае решения только соответствующее уравнение домножим на −1). Поскольку симплекс-метод может применяться для специальных задач ЛП, то следует ответит на вопрос: существует ли специальная задача, эквивалентная канонической задаче (1)? И в случае положительного ответа указать способ ее построения. Рассмотрим вспомогательную задачу ЛП: h( x ) = −(t1 + t2 + … + tm ) → max;   a11 x1 + a12 x2 + … + a1n xn + t1 = b1 ;   ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯  a x + a x +…+ a x + tm = bm ; m1 1 m2 2 mn n   x j ≥ 0, j = 1, 2, … , n; ti ≥ 0, i = 1, 2, … , m, где x = ( x1 , x2 , … , xn , t1 , … , tm ) − (n+m)-мерный вектор. (2) Переменные ti , i = 1, 2, … , m, называются искусственными базисными переменными. З а м е ч а н и е. Задача (2) легко приводится к специальному виду. Для этого функцию h(x ) надо выразить через небазисные переменные x j , j = 1, 2,… , n, используя линейные ограничения задачи (2). Кроме того, отметим, что h( x ) ≤ 0 при любых допустимых значениях переменных. Теорема 1 (критерий существования допустимых решений канонической задачи ЛП). Для канонической задачи (1) множество D допустимых решений не пусто тогда и только тогда, когда оптимальное значение целевой функции h( x ) вспомогательной задачи (2) равно нулю. Теорема 2 (о преобразовании канонич. задачи ЛП в эквивалентную ей специальную задачу ЛП). Если множество D допустимых решений канонической задачи (1) не пусто, то существует эквивалентная ей специальная задача ЛП, которая соответствует симплексной таблице, полученной преобразованием завершающей симплекс. таблицы вспомогательной задачи (2). Метод построения специальной задачи ЛП, эквивалентной канонической задаче ЛП, с помощью вспомогательной задачи называется методом искусственного базиса. З а м е ч а н и е. Решение симплекс-таблицы • преобразования симплекс-таблицы, стараясь выводить из базиса дополнительные переменные. • если какая-то дополнительная переменная выведена из базиса, то соответствующий столбец симплекс-таблицы можно просто вычеркнуть и больше к нему не возвращаться. • В конце возможны два варианта. 1. Все векторы, соответствующие введенным дополнительным переменным, будут выведены из базиса. В этом случае мы просто вернемся к исходной задаче, попав в какую-то вершину допустимой области. Все столбцы симплекс-таблицы, соответствующие дополнительным переменным, тогда исчезнут и дальше будет решаться исходная задача. 2. Получающийся оптимальный план будет все-таки содержать какую-то из дополнительных переменных. Это означает, что допустимая область исходной задачи пуста, то есть ограничения исходной задачи противоречивы и поэтому исходная задача вообще не имеет решений. 12. ДВОЙСТВЕННОСТЬ В ЛИНЕЙНОМ ПРОГРАММИРОВАНИИ Рассмотрим пару задач ЛП вида: I II n m f ( x) = ∑ c j x j → max; g ( y ) = ∑ bi yi → min; j =1 n ∑a x ij j ≤ bi , i = 1, … , k ; j = bi , i = k + 1, … , m; i =1 ↔ yi ≥ 0, i = 1, … , k ; j =1 n ∑a x ij yi − любые, i = k + 1, … , m; j =1 x j ≥ 0, j = 1, … , l ; m ↔ ∑a y ≥ c , ij i j j = 1, … , l ; i =1 x j − любые, j = l + 1, … , n. m ∑a ij yi = c j , j = l + 1, … , n. i =1 Задачу I называют прямой задачей ЛП, а задачу II − двойственной. Однако задача, двойственная к задаче II, есть исходная прямая задача I, поэтому можно считать любую из такой пары задач прямой, а другую − двойственной. Неравенства задач I и II, соответствующие друг другу, называются сопряженными (отмечены стрелками). Множества допустимых решений задач I и II будем обозначать DI и DII соответственно. Теорема 1 (первая теорема двойственности). Если одна из пары двойственных задач I, II разрешима, то разрешима и другая задача, причем f ( x∗ ) = g ( y ∗ ) , где x∗ , y ∗ − оптимальные планы задач I, II соответственно. Если одна из пары задач I, II не имеет решения из-за неограниченности целевой функции, то у другой задачи множество допустимых решений пусто. Говорят, что планы x ∈ DI , y ∈ DII удовлетворяют условиям дополняющей нежесткости (УДН), если при подстановке этих векторов в ограничения задач I, II хотя бы одно из любой пары сопряженных неравенств обращается в равенство. Теорема 2 (вторая теорема двойственности). Планы x∗ ∈ DI , y ∗ ∈ DII оптимальны в задачах I, II тогда и только тогда, когда они удовлетворяют условиям дополняющей нежесткости. 13. ТЕОРИЯ ГРАФОВ 1. НЕКОТОРЫЕ ОПРЕДЕЛЕНИЯ ИЗ ТЕОРИИ ГРАФОВ На практике часто возникают задачи на конфигурациях, состоящих из точек и соединяющих их линий. При этом является несущественным, прямые это линии или непрерывные кривые, соединяющие две концевые точки; где расположены эти линии и длинные они или короткие. Существенно то, что линии соединяют две данные точки. Примерами могут служить электрическая схема, графическое изображение иерархического строения какой-либо компании или государства, схематическое изображение молекулы некоторого вещества и т. п. Определение 1. Графом G называется математический объект, состоящий из двух множеств: непустого конечного множества V элементов, называемых вершинами, и множества Е неупорядоченных и упорядоченных пар вершин. Обозначается граф G = (V, E). Неупорядоченная пара вершин называется ребром, упорядоченная − дугой. Определение 2. Граф, содержащий только ребра, называется неориентированным, а граф, содержащий только дуги, − ориентированным, или сетью. Обычно граф представляют на плоскости диаграммой, изображая его вершины точками, а ребра или дуги − линиями, соединяющими соответствующие пары точек. Например, две диаграммы, изображенные на рис. 1, представляют один и тот же неориентированный граф, имеющий четыре вершины и пять ребер. Ориентированный граф с пятью вершинами и пятью дугами изображен на рис. 2. v1 v2 v1 v2 v4 v3 v4 v3 а б Рис. 1 Определение 3. Если e = (u, v) − ребро (дуга) графа G, то вершины u и v называются смежными и инцидентными ребру (дуге) е. Ориентированный граф можно задать с помощью таблицы, устанавливающей инцидентность всех его дуг и вершин. v2 Определение 4. Степенью вершины v называется число инцидентных ей ребер (дуг). e1 e2 v3 e3 Обозначается степень вершины σ (v) . Определение 5. Граф называется простым, если в нем нет петель и кратных ребер, т. е. любая пара вершин инцидентна не v1 e4 v4 е5 v5 Рис. 2 более чем одному ребру (дуге) и любое ребро (дуга) соединяет две несовпадающие вершины. Определение 6. Подграфом графа G называется всякий граф G′ = (V ′, E ′) с множеством вершин V ′ ⊆ V и множеством ребер (дуг) E ′ ⊆ E , каждое из которых инцидентно только вершинам из V ′ . Если не учитывать ориентацию дуг графа, то все они становятся его ребрами. Определение 8. Маршрутом в графе G называется последовательность ребер e1 , e2 , … , er , в которой каждые два соседних ребра имеют одну общую вершину. Маршрут называется цепью, если в нем все ребра различны, и простой цепью, если все ее вершины различны. Определение 9. Граф G называется связным, если любая пара его вершин соединена простой цепью. Максимальный связный подграф графа G называется компонентой связности графа G. Таким образом, несвязный граф имеет по меньшей мере две компоненты связности. Определение 10. Путем из вершины i в вершину j ориентированного графа G называется такая последовательность дуг e1 , e2 , … , er , в которой вершина i является началом дуги e1 , вершина j − концом дуги er и конец дуги ek совпадает с началом дуги ek +1 , k = 1, … , r − 1 . Количество дуг r называется длиной пути. Если начало пути совпадает с его концом (i = j), то путь называется контуром. Примерами пути в графе, изображенном на рис. 2, являются последовательности дуг e4 , e1 , e2 ; e3 , e5 и др. Лишь один путь представляет собой контур: e1 , e3 , e4 . 2.ПОНЯТИЕ СЕТЕВОГО ГРАФИКА Выполнение комплексных исследований, проектирование и строительство крупных объектов и многие другие виды деятельности требуют календарной увязки большого числа взаимосвязанных работ, выполняемых различными организациями. Составление и анализ соответствующих календарных планов − весьма сложная задача, при решении которой используются так называемые методы сетевого планирования. Сеть из m вершин – ориентированный граф, в котором выделены две вершины – вход и выход. Зависимость между отдельными видами работы некоторого комплекса задается с помощью ориентированного графа Г(m, n, S), где m − количество вершин, n − количество дуг, S = { ik , jk }k =1, …, n − ориентированные дуги. Каждую отдельную дугу будем обозначать s = { is , js } , она соответствует определенной работе. Если дуги s1 , s2 , … , sl образуют путь из вершины is1 в вершину jsl , то это означает, что выполнение работ s, для которых is = jsl , может начаться только после завершения всех работ: s1 , s2 , … , sl . Следовательно, каждую вершину i графа Г можно интерпретировать как событие, означающее завершение всех работ s, для которых i = js , и возможность начала выполнения всех работ s, для которых i = is . Очевидно, что граф Г не имеет контуров. Обозначим множества дуг, входящих в вершину i и выходящих из нее, Ni + и Ni − соответственно. Так как граф Г не имеет контуров, то в нем найдется вершина i0 , в которую не входит ни одна дуга, и вершина i00 , из которой ни одна дуга не выходит: N i0 + = ∅, N i00 − = ∅. Будем считать, что i0 = 1 и это единственная начальная вершина, а i00 = m − единственная конечная вершина. Событие i0 − начало выполнения работ комплекса, i00 − завершение всех работ. Понятно, что для любой вершины i ≠ 1 в графе Г существует путь из вершины i0 = 1 в вершину i и для любой вершины i ≠ m − путь из этой вершины в вершину i00 = m. Пример графа Г(8, 14, S) приведен на рис. 3. 4 Пусть известно время τ s , 7 необходимое 1 для выполнения каждой работы s = 1, … , n . В этом 3 8 6 случае говорят, что имеется сетевой график Г ( m, n, S , T ) , где T = { τ s }s =1, …, n . Календарный план 5 2 выполнения сетевого определяется Рис. 3 графика вектором t = ( t1 , … , tm ) , компоненты которого означают планируемые сроки наступления соответствующих событий. План t является допустимым, если t js − tis ≥ τ s , s = 1, … , n . При этом время, затраченное на выполнение всего комплекса работ, λ (t ) = tm − t1 . При составлении календарного плана необходимо поставить цель минимизации величины λ (t ) . 3. ЗАДАЧА О КРАТЧАЙШЕМ СРОКЕ Обозначим через Θ множество всех допустимых векторов t. Задача о кратчайшем сроке состоит в том, чтобы при заданном сетевом графике Г ( m, n, S , T ) найти такой календарный план t 0 , для которого ( ) λ t 0 = min λ ( t ) . t∈Θ (1) При этом величина λ0 = λ ( t 0 ) = min λ ( t ) называется кратчайшим сроком. t∈Θ Обозначим через Pi + для i ≠ 1 множество всех путей в сети Г из вершины 1 в вершину i, а через Pi − для i ≠ m − множество всех путей из вершины i в вершину m. Рассмотрим вектор t 0 , имеющий следующие координаты: t10 = 0 ; ti 0 = max ∑ τ s , i = 2, … , m . + (2) p∈Pi s∈ p Очевидно, что t js 0 − tis 0 ≥ τ s , s = 1, … , n ; tm 0 = λ0 . Следовательно, такой вектор t 0 является решением задачи о кратчайшем сроке. Заметим, что для вычисления координат этого вектора не надо перебирать все пути из множеств Pi + , так как соотношения (2) эквивалентны следующим: ( ) t10 = 0 ; ti 0 = max tis 0 + τ s , i = 2, … , m . + s∈Ni (3) Алгоритм решения задачи о кратчайшем сроке Шаг 0. Полагаем ti = 0, i = 1, … , m . Шаг 1. Полагаем М = 0. Для s = 1, … , n делаем следующее: 1) определяем t js ′ = max { t js , tis + τ s } ; 2) если t js ′ ≠ t js , то полагаем М = 1; принимаем t js = t js ′ . Шаг 2. Если М = 0, т. е. значения координат вектора t не изменились, то полагаем ti 0 = ti , i = 1, … , m . Конец работы алгоритма. Если М ≠ 0, то переходим на шаг 1. Замечание. Пусть l0 − длина максимального (по числу дуг) пути из вершины 1 в вершину m. Очевидно, что не более чем за l0 просмотров списка дуг сетевого графика все координаты вектора t 0 ЗАДАЧА ЛП, ЭКВИВАЛЕНТНАЯ ЗАДАЧЕ О КРАТЧАЙШЕМ СРОКЕ Введем переменные ti , i = 1, … , m . Тогда задача о кратчайшем сроке как задача ЛП имеет следующую постановку: tm − t1 → min ; t − t ≥ τ , s = 1, … , n. s  js is (I) 4. ЗАДАЧА О КРИТИЧЕСКОМ ПУТИ Координаты рассмотренного в предыдущем пункте вектора t 0 определяют максимально ранние возможные сроки наступления всех событий. При анализе сетевого графика также важно установить имеющиеся резервы времени наступления событий от первого до (m-1)-го, т. е. определить некоторый интервал времени, в течение которого наступление данного события не повлияет на время завершения всего комплекса работ. Для этого найдем максимально поздние возможные сроки наступления событий, не изменяющие планируемый срок окончания комплекса работ. Обозначим их ti 00 , i = 1, … , m − 1 . Рассмотрим множество вершины i в {P} i − i =1, …, m −1 вершину , где Pi − − совокупность всех путей из m. Вычислим tm 00 = tm 0 ; величины: ti 00 = tm 00 − max ∑ τ s , i − 1, … , m − 1 . p∈Pi − s∈ p Очевидно, что Pm + = P1− , следовательно, t100 = tm 00 − max ∑ τ s = tm 00 − max ∑ τ s = tm 00 − tm 0 = 0 = t10 . p∈P1− s∈ p p∈Pm + s∈ p Заметим, что для i = 2, … , m − 1 ti 0 ≤ ti 00 , поэтому для тех же значений i ρi = ti 00 − ti 0 ≥ 0 представляют собой резервы времени для наступления соответствующих событий. Определение 11. Событие i сетевого графика Г(m, n, S, T) называется критическим, если резерв времени для выполнения этого события равен нулю: ρi = 0 , т. е. ti 00 = ti 0 . Заметим, что ti 00 = min ( t js 00 − τ s ) , i = 1, … , m − 1 , где − τs is js s∈Ni Ni − − множество дуг, выходящих из вершины i. τ s1 τ s2 Определение 12. Работа s сетевого графика Г(m, n, S, T) называется напряженной, если события is , js являются критическими и t js 0 = tis 0 + τ s (или, что то же самое, tis 00 = t js 00 − τ s ). Рис. 4 Заметим, что из критичности событий is , js еще не следуют приведенные в определении равенства. Например, для подграфа, изображенного на рис. 4, может быть τ s1 + τ s2 > τ s . Алгоритм вычисления резервов времени Шаг 0. Полагаем ti = tm 0 , i = 1, … , m . Шаг 1. Полагаем М = 0. Для s = n, … , 1 выполняем следующие действия: 1) определяем tis ′ = min { tis , t js − τ s } ; 2) если tis ′ ≠ tis , то полагаем М = 1; принимаем tis = tis ′ . Шаг 2. Если М = 0 (т. е. значения ti не изменялись), то полагаем ti 00 = ti , ρi = ti 00 − ti 0 , i = 1, … , m . Конец работы алгоритма. Если М ≠ 0, то переходим на шаг 1. Замечание. Пусть l0 − длина максимального (по числу дуг) пути из вершины 1 в вершину m. Очевидно, что не более чем за l0 просмотров списка дуг сетевого графика все координаты вектора ρ будут определены. Определение 13. Продолжительностью пути p = {s1 , s2 , … , sl } называется время, необходимое для выполнения всех работ, лежащих на этом пути: l τ ( p) = ∑ τ s = ∑τ sk . s∈ p (5) k =1 Определение 14. Путь, имеющий наибольшую продолжительность, называется критическим. Так как t10 = 0 , ti 0 = max ∑τ s ( i = 2, … , m ) и + p∈Pi s∈ p tm 0 = λ0 , то продолжительность любого критического пути совпадает с кратчайшим сроком выполнения всего комплекса работ: max+ τ ( p ) = min λ (t ) . p∈Pm (6) t∈Θ От продолжительности работ, лежащих на критическом пути (критических работ), зависит срок завершения всех работ. Сокращение или увеличение сроков выполнения критических работ соответственно продолжительность выполнения всего проекта. сокращает или увеличивает Теорема. Для критичности пути p = {s1 , s2 , … , sl } из вершины i0 = 1 в i00 = m необходимо и достаточно, чтобы все его дуги отвечали напряженным работам. Доказательство. Рассмотрим путь p = {s1 , s2 , … , sl } из вершины i0 = 1 в i00 = m . Пронумеруем его вершины: k0 = 1; sr = ( kr −1 , kr ) , r = 1, … , m . Напомним, что работа s называется напряженной, если: а) события is , js являются критическими, т. е. ρis = ρ js = 0 ; (7) б) выполняется равенство: t js 0 = tis 0 + τ s ( t js 00 = tis 00 + τ s ). (8) Доказательство теоремы проведем в два этапа. 1) Докажем, что путь р является критическим тогда и только тогда, когда для всех его дуг выполняются соотношения (8). Из допустимости векторов t 0 , t 00 следует, что при r = 1, … , l tk0r − tk0r −1 ≥ τ sr , tk00r − tk00r −1 ≥ τ sr . (9) Оценим время, затраченное на выполнение всего комплекса работ: l l r =1 r =1 ( ) λ ( p) = ∑ τ s = ∑τ sr ≤ ∑ tk0r − tk0r −1 = s∈ p = tk01 − tk00 + tk02 − tk01 + … + tk0l − tk0l −1 = tk0l − tk00 = tm0 − t10 = tm0 = λ0 . (10) Аналогично получаем для вектора t 00 : l l r =1 r =1 ( ) λ ( p) = ∑ τ s = ∑τ sr ≤ ∑ tk00r − tk00r −1 = tm00 − t100 = tm00 = λ0 . s∈ p (11). Путь р критичен в том и только в том случае, если λ ( p) = λ0 , следовательно, все неравенства в формулах (8)─(10) выполняются как равенства, т. е. tk0r − tk0r −1 = τ sr , tk00r − tk00r −1 = τ sr и соотношения (8) верны для всех дуг пути р. Таким образом установлено, что путь является критическим тогда и только тогда, когда все его дуги отвечают напряженным работам. 2) Докажем, что для пути, ведущего из вершины i0 = 1 в i00 = m , условия (7) и (8) эквивалентны. Действительно, так как k0 = 1 , то tk00 = t10 = t100 = tk000 = 0 . Тогда из соотношений (22) следует, что tk01 = tk00 + τ s1 = tk000 + τ s1 = tk001 и ρ k1 = 0 . Аналогично устанавливаем, что условия (7) выполняются для всех остальных вершин пути р. Таким образом, из соотношений (8) вытекает справедливость формул (7). Обратная связь также очевидна. Теорема доказана. ЗАДАЧА ЛП, ЭКВИВАЛЕНТНАЯ ЗАДАЧЕ О КРИТИЧЕСКОМ ПУТИ Исходные данные содержатся в сетевом графике Г(m, n, S, T). Введем переменные: 1, если дуга s входит в критический путь; xs =  0, если не входит. Тогда задача ЛП, эквивалентная задаче о критическом пути, имеет вид: n   ∑τ s xs → max ;  s =1 − ∑ xs = −1;  s∈N1−   ∑ xs − ∑ xs = 0, i = 2, … , m − 1;  s∈N + s∈Ni−  i  ∑ xs = 1;  s∈Nm+   xs ≥ 0, s = 1, … , n;  xs ∈ {0; 1} , s = 1, … , n. Заметим без доказательства, что отбрасывание (II) последнего условия ( xs ∈ {0; 1} ; s = 1, … , n ) не изменяет оптимального решения задачи (II), т. е. исключив это условие, ее можно решать как обычную задачу ЛП. Нетрудно установить, что задача (II) (без последнего условия) является двойственной к задаче (I).