Введение
К дифференциальным уравнениям с частными производнымиприходим при решении самых разнообразных задач. Например, при помощи дифференциальныхуравнений с частными производными можно решать задачи теплопроводности,диффузии, многих физических и химических процессов.
Как правило, найти точное решение этих уравнений неудается, поэтому наиболее широкое применение получили приближенные методы ихрешения. В данной работе ограничимся рассмотрением дифференциальных уравнений счастными производными второго порядка, а точнее дифференциальными уравнениями счастными производными второго порядка параболического типа, когда эти уравненияявляются линейными, а искомая функция зависит от двух переменных
Для решения дифференциальных уравнений параболическоготипа существует несколько методов их численного решения на ЭВМ, однако особоеположение занимает метод сеток, так как он обеспечивает наилучшие соотношенияскорости, точности полученного решения и простоты реализации вычислительногоалгоритма. Метод сеток еще называют методом конечных разностей.
1 Теоретическая часть
1.1 Постановка задач для уравнений параболическоготипа
Классическим примером уравнения параболического типаявляется уравнение теплопроводности (диффузии). В одномерном по пространствуслучае однородное (без источников энергии) уравнение теплопроводности имеет вид
/> /> /> (1)
Если на границах />и />заданы значения искомой функции />в виде
/> />, />, (2)
/> />, />, (3)
т.е. граничные условия первого рода, и, кроме тогозаданы начальные условия
/> />, />, (4)
то задачу (1)-(4) называют первой начально-краевойзадачей для уравнения теплопроводности (1).
В терминах теории теплообмена />— распределение температуры в пространственно-временной области
/> a2 — коэффициент температуропроводности, а (2), (3) спомощью функций />,/>задают температуру на границах /> и />.
Если на границах /> и /> заданы значения производныхискомой функции по пространственной переменной:
/> />, />, (5)
/> />, />, (6)
т.е. граничные условия второго рода, то задачу (1),(5), (6), (4) называют второй начально-краевой задачей для уравнениятеплопроводности (1). В терминах теории теплообмена на границах в этом случаезаданы тепловые потоки.
Если на границах заданы линейные комбинации искомойфункции и ее производной по пространственной переменной:
/> />, />, (7)
/> />, />, (8)
т.е. граничные условия третьего рода, то задачу (1),(7), (8), (4) называют третьей начально-краевой задачей для уравнениятеплопроводности (1). В терминах теплообмена граничные условия (7), (8) задаюттеплообмен между газообразной или жидкой средой с известными температурами />на границе/>и /> на границе /> и границами расчетной области снеизвестными температурами />, />. Коэффициенты α, β –известные коэффициенты теплообмена между газообразной или жидкой средой исоответствующей границей.
Для пространственных задач теплопроводности в области /> перваяначально-краевая задача имеет вид
/> (9)
Аналогично ставится вторая и третья начально-краевыезадачи для пространственного уравнения (9). На практике часто ставятсяначально-краевые задачи теплопроводности со смешанными краевыми условиями,когда на границах задаются граничные условия различных родов.
1.2 Основные определения и конечно-разностные схемы
Основные определения, связанные с методом конечныхразностей, рассмотрим на примере конечно-разностного решения первойначально-краевой задачи для уравнения теплопроводности (1)-(4).
Согласно методу сеток в плоской области Dстроится сеточная область Dh, состоящая из одинаковых ячеек. При этом область Dh должна как можно лучше приближать область D.Сеточная область (то есть сетка) Dh состоит изизолированных точек, которые называются узлами сетки. Число узлов будетхарактеризоваться основными размерами сетки h: чемменьше h, тем больше узлов содержит сетка. Узел сеткиназывается внутренним, если он принадлежит области D, а всесоседние узлы принадлежат сетке Dh. В противномслучае он называется граничным. Совокупность граничных узлов образует границу сеточнойобласти Гh.
Сетка может состоять из клеток разной конфигурации:квадратных, прямоугольных, треугольных и других. После построения сеткиисходное дифференциальное уравнение заменяется разностным уравнением во всехвнутренних узлах сетки. Затем на основании граничных условий устанавливаютсязначения искомого решения в граничных узлах. Присоединяя граничные условия сеточнойзадачи к разностным уравнениям, записанных для внутренних узлов, получаемсистему уравнений, из которой определяем значения искомого решения во всехузлах сетки.
Нанесем на пространственно-временную область />, /> конечноразностную сетку ωh,τ:
/>(10)
с пространственным шагом h=l/N и шагом по времени τ=T/K.
/>
Рисунок 1 – Конечно-разностная сетка
Введем два временных слоя: нижний />, на которомраспределение искомой функции u(xj,tk), />, известно (при к = 0распределение определяется начальным условием (4)u(xj,tk)=ψ(xj)), иверхний временной слой tk+1=(k+1) τ, на которомраспределение искомой функции u(xj,tk+1), />.
Сеточной функцией задачи (1)-(4) называют однозначное отображение целых аргументов j,k в значения функции />.
На введенной сетке вводят сеточные функции />, /> первая из которыхизвестна, вторая подлежит определению. Для определения в задаче (1)-(4) заменяют(аппроксимируют) дифференциальные операторы отношением конечных разностей(более подробно это рассматривают в разделах численных методов «Численноедифференцирование»), получают
/>, (11)
/>, (12)
Подставляя (11), (12) в задачу (1)-(4), получим явнуюконечно-разностную схему для этой задачи в форме
/>(13)
В каждом уравнении этой задачи все значения сеточнойфункции известны, за исключением одного, />,которое может быть определено явно из соотношений (13). В соотношения (13)краевые условия входят при значениях j=1 и j=N-l,a начальное условие – при k = 0.
Если в (12) дифференциальный оператор попространственной переменной аппроксимировать отношением конечных разностей наверхнем временном слое:
/>, (14)
то после подстановки (11), (14) в задачу (1)-(4)получим неявную конечно-разностную схему для этой задачи:
/>(15)
Теперь сеточную функцию/>наверхнем временном слое можно получить из решения (15) с трехдиагональнойматрицей. Эта СЛАУ в форме, пригодной для использования метода прогонки, имеетвид
/>
где
/>/>;
/> />;
/>, />;
/>/>;
/>;
/>;
/>.
Шаблоном конечно-разностной схемы называют еегеометрическую интерпретацию на конечно-разностной сетке. На рисунке приведенышаблоны для явной и неявной конечно-разностных схем при аппроксимации задачи.
/>
Рисунок 2 — Шаблон явной конечно-разностной схемы для уравнения теплопроводности
/>
Рисунок 3 — Шаблон неявной конечно-разностной схемы для уравнения теплопроводности
В случае явных схем значения функции в узле очередногослоя можно найти, зная значения в узлах предыдущих слоев. В случае неявных схемдля нахождения значений решения в узлах очередного слоя приходится решать системууравнений. Для проведения вычислений самой простой схемой оказывается первая:достаточно на основании начального условия найти значения функции в узлах слоя />, чтобы вдальнейшем последовательно определять значения решения в узлах слоев /> и т.д. Вслучае второй схемы, которая является неявной, обязательно приходится решатьсистему уравнений для нахождения решения сеточной задачи. В любом случае согласнометоду сеток будем иметь столько уравнений, сколько имеется неизвестных(значения искомой функции в узлах). Число неизвестных равно числу всех узловсетки. Решая систему уравнений, получаем решение поставленной задачи.
Разрешимость этой системы для явных схем вопросов невызывает, так как все действия выполняются в явно определеннойпоследовательности. В случае неявных схем разрешимость системы следуетисследовать в каждом конкретном случае. Важным вопросом является вопрос о том,на сколько найденные решения хорошо (адекватно) отражают точные решения, иможно ли неограниченно сгущая сетку (уменьшая шаг по осям) получитьприближенные решения, сколь угодно близкие к точным решениям? Это вопрос осходимости метода сеток.
На практике следует применять сходящиеся разностныесхемы, причем только те из них, которые являются устойчивыми, то есть прииспользовании которых небольшие ошибки в начальных или промежуточныхрезультатах не приводят к большим отклонениям от точного решения. Всегдаследует использовать устойчивые разностные схемы, проводя соответствующие исследованияна устойчивость. Явные схемы просты для организации вычислительного процесса,но имеют один весьма весомый недостаток: для их устойчивости приходится накладыватьсильные ограничения на сетку. Неявные схемы свободны от этого недостатка, ноесть другая трудность – надо решать системы уравнений большой размерности, чтона практике при нахождении решения сложных уравнений в протяженной области свысокой степенью точности может потребовать больших объемов памяти ЭВМ ивремени на ожидание конечного результата. К счастью, прогресс не стоит на местеи уже сейчас мощности современных ЭВМ вполне достаточно для решенияпоставленных перед ними задач.
Вопрос устойчивости будет рассмотрен далее.
1.3 Аппроксимация
Из определения порядка аппроксимации ясно, что чемвыше порядок аппроксимации, тем лучше конечно-разностная схема приближается кдифференциальной задаче. Это не означает, что решение по разностной схеме можетбыть так же близко к решению дифференциальной задачи, так как разностная схемаможет быть условно устойчивой или абсолютно неустойчивой вовсе.
Для нахождения порядка аппроксимации используетсяаппарат разложения в ряды Тейлора точных (неизвестных, но дифференцируемых)решений дифференциальной задачи в узлах сетки (подчеркнем: значения сеточнойфункции uh дискретны, следовательно, не дифференцируемы ипоэтому не разлагаются в ряды Тейлора).
1.4 Устойчивость.Исследование устойчивости методом гармонического анализа
конечно-разностнаясхема устойчива, если для малых возмущений входных данных (начально-краевыхусловий и правых частей) конечно-разносная схема обеспечивает малые возмущениясеточной функции uh т.е. решение с помощьюконечно-разностной схемы находится под контролем входных данных.
Если вовходные данные fn входят тольконачальные условия или только краевые условия, или только правые части, тоговорят об устойчивости соответственно по начальным условиям, по краевымусловиям или по правым частям.
Изматематической физики известно, что решение начально-краевых задачпредставляется в виде следующего ряда:
/>, (16)
где λn – собственные значения
/> – собственныезначения функции, получаемые из решения соответствующей задачи Штурма-Лиувиля,т.е. решение может быть представлено в виде суперпозиции отдельных гармоник />, каждая изкоторых есть произведение функции времени и функции пространственнойпеременной, причем последняя по модулю ограничена сверху единицей при любыхзначениях переменной x.
В то же времяфункция времени />, называемая амплитудной частьюгармоники, никак не ограничена, и, по всей вероятности, именно амплитуднаячасть гармоник является источником неконтролируемого входными данными ростафункции и, следовательно, источником неустойчивости.
Такимобразом, если конечно-разностная схема устойчива, то отношение амплитуднойчасти гармоники на верхнем временном слое к амплитудной части на нижнемвременном слое по модулю должно быть меньше единицы.
Еслиразложить значение сеточной функции /> в ряд Фурье по собственнымфункциям:
/> (17)
гдеамплитудная часть/>может быть представлена в видепроизведения
/> (18)
где />– размерный ипостоянный сомножитель амплитудной части,
k – показатель степени(соответствующий номеру временного слоя) сомножителя, зависящего от времени.
Тогда подставив(17) в конечно-разностную схему, можно по модулю оценить отношение амплитудныхчастей на соседних временных слоях.
Однакопоскольку операция суммирования линейна и собственные функции ортогональны дляразличных индексов суммирования, то в конечно-разностную схему вместо сеточныхзначений достаточно подставить одну гармонику разложения (17) (при этом уамплитудной части убрать индекс n), т.е.
/> (19)
Таким образом,если конечно-разностная схема устойчива по начальным данным, то
/>, (20)
т. е. условие(20) является необходимым условием устойчивости.
1.5 СхемаКранка-Николсона
параболическое дифференциальное уравнение конечная разность
Явная конечноразностная схема, записанная в форме
/> (21)
обладает темдостоинством, что решение на верхнем временном слое tk+l получаетсясразу (без решения СЛАУ) по значениям сеточной функции на нижнем временном слоеtk, где решение известно (при k = 0 значения сеточной функции формируются из начальногоусловия). Но эта же схема обладает существенным недостатком, поскольку онаявляется условно устойчивой. С другой стороны, неявная конечно-разностнаясхема, записанная форме
/> (22)
приводит кнеобходимости решать СЛАУ, но зато эта схема абсолютно устойчива.
Проанализируемсхемы (21) и (22). Пусть точное решение, которое неизвестно, возрастает повремени, т.е. />. Тогда, в соответствии с явнойсхемой (21), разностное решение будет заниженным по сравнению с точным, так как/>определяетсяпо меньшим значениям сеточной функции на предыдущем временном слое, посколькурешение является возрастающим по времени.
Для неявнойсхемы (22) на возрастающем решении, наоборот, решение завышено по сравнению сточным, поскольку оно определяется по значениям сеточной функции на верхнемвременном слое.
На убывающемрешении картина изменяется противоположным образом: явная конечно-разностнаясхема завышает решения, а неявная — занижает (Рисунок 4).
На основеэтого анализа возникла идея о построении более точной неявно-явнойконечно-разностной схемы с весами при пространственных конечно-разностныхоператорах, причем при измельчении шагов тик точное (неизвестное) решение можетбыть взято в «вилку» сколь угодно узкую, так как если явная и неявная схемыаппроксимируют дифференциальную задачу и эти схемы устойчивы, то при стремлениисеточных характеристик τ и h к нулю решения по явной и неявнойсхемам стремятся к точному решению с разных сторон.
/>
Рисунок 4 –Двусторонний метод аппроксимации
Проведенныйанализ дал блестящий пример так называемых двусторонних методов, исследованныхВ. К. Саульевым
Рассмотримнеявно-явную схему с весами для простейшего уравнения теплопроводности:
/> (23)
где θ– вес неявной части конечно-разностной схемы,
θ-1– вес для явной части
Причем />. При θ=1имеем полностью неявную схему, при θ=0 – полностью явную схему, апри θ=1/2 – схему Кранка-Николсона.
Всоответствии с гармоническим анализом для схемы (23) получаем неравенство
/>,
откуда
/> (24)
причем правоенеравенство выполнено всегда.
Левое неравенствоимеет место для любых значений σ, если />. Если же вес θ лежитв пределах />,то между σ и θ из левого неравенства устанавливаетсясвязь
/> /> (25)
являющаясяусловием устойчивости неявно-явной схемы с весами (23), когда вес находится впределах />.
Такимобразом, неявно-явная схема с весами абсолютно устойчива при />и условно устойчива сусловием (25) при />.
Рассмотримпорядок аппроксимации неявно-явной схемы с весами, для чего разложим в рядТейлора в окрестности узла (xj,tk) на точномрешении значения сеточных функций /> по переменной t, />, /> по переменной хи полученные разложения подставим в (23):
/>
В этомвыражении дифференциальный оператор /> от квадратной скобки всоответствии с дифференциальным уравнением равен дифференциальному оператору />, всоответствии с чем вышеприведенное равенство приобретает вид
/>
Послеупрощения получаем
/>,
откуда видно,что для схемы Кранка-Николсона (θ = 1/2) порядок аппроксимациисхемы (23) составляет />, т.е. на один порядок по временивыше, чем для обычных явных или неявных схем. Таким образом, схемаКранка-Николсона при θ = 1/2 абсолютно устойчива и имеет второйпорядок аппроксимации по времени и пространственной переменной х.
Используем вуравнение (23) подстановку r=a2k/h2. Но в то же время его нужно решить для трех«еще не вычисленных» значений />, />, и />. Это возможно, если все значенияперенести в левую часть уравнения. Затем упорядочим члены уравнения (23) и врезультате получим неявную разностную формулу
/> (26)
для i=2,3,…,n-1. Члены вправой части формулы (26) известны. Таким образом, формула (26) имеет видлинейной трехдиагональной системы АХ=В. Шесть точек, используемых в формулеКранка-Николсона (26), вместе с промежуточной точкой решетки, на которой основанычисленные приближения, показаны на рисунке 5.
/>
Рисунок 5 –Шаблон (схема) метода Кранка-Николсона
Иногда вформуле (26) используется значение r=1. В этомслучае приращение по оси t равно />, формула (26)упрощается и принимает вид
/>, (27)
для i=2,3,…,n-1. Граничныеусловия используются в первом и последнем уравнениях (т. е. в /> и /> соответственно).
Уравнения(27) особенно привлекательны при записи в форме трехдиагональной матрицы АХ =В.
Если методКранка-Николсона реализуется на компьютере, то линейную систему АХ = В можнорешить либо прямым методом, либо итерационным.
/>
2. Практическая часть
2.1 Постановказадачи
Используемметод Кранка-Николсона, чтобы решить уравнение
/>,
где xϵ(0;1),
tϵ(0;0.1),
с начальнымусловием
/>,
где t=0,
xϵ(0;1),
и граничнымиусловиями
u(0,t) = g1(t) ≡ 0,
u(1,t) = g2(t) ≡ 0.
2.2Решение в ППП MatLab
Решение будемискать в ППП MatLab 7. Создадим четыре выполняемых m-фала: crnich.m – файл-функция с реализацией методаКранка-Николсона; trisys.m – файл-функция метода прогонки; f.m – файл-функциязадающая начальное условие задачи; fе.m – файл-функция задающая функциюопределяющую точное решение задачи(найдена аналитическим путем). Листинги программпредставлены в приложении А.
Для простотывозьмем шаг Δх = h = 0,1 и Δt = к = 0,01. Таким образом, соотношение r =1. Пусть решетка имеет n=11столбцов в ширину и m=11 рядов в высоту.
2.3 Анализрезультатов
Решения дляданных параметров отразим в таблице 1. Трехмерное изображение данных из таблицыпокажем на рисунке 5.
Таблица 1 – Значенияu(хi, ti),полученные методом Кранка-Николсона
xi 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ti 1.1180 1.5388 1.1180 0.3633 0.3633 1.1180 1.5388 1.1180 0.01 0.6169 0.9288 0.8621 0.6177 0.4905 0.6177 0.8621 0.9288 0.6169 0.02 0.3942 0.6480 0.7186 0.6800 0.6488 0.6800 0.7186 0.6480 0.3942 0.03 0.2887 0.5067 0.6253 0.6665 0.6733 0.6665 0.6253 0.5067 0.2887 0.04 0.2331 0.4258 0.5560 0.6251 0.6458 0.6251 0.5560 0.4258 0.2331 0.05 0.1995 0.3720 0.4996 0.5754 0.6002 0.5754 0.4996 0.3720 0.1995 0.06 0.1759 0.3315 0.4511 0.5253 0.5504 0.5253 0.4511 0.3315 0.1759 0.07 0.1574 0.2981 0.4082 0.4778 0.5015 0.4778 0.4082 0.2981 0.1574 0.08 0.1419 0.2693 0.3698 0.4338 0.4558 0.4338 0.3698 0.2697 0.1419 0.09 0.183 0.2437 0.3351 0.3936 0.4137 0.3936 0.3351 0.2437 0.1283 0.1 0.1161 0.2208 0.3038 0.3570 0.3753 0.3570 0.3038 0.2208 0.1161
Величины,полученные методом Кранка-Николсона, достаточно близки к
аналитическомурешению u(x,t) = sin(πx)e-π2t+ sin(3πx)e-9π2t, истинные значения для последнего представлены втаблице 2
Максимальнаяпогрешность для данных параметров равна 0,005
Таблица 2 –точные значения u(хi,ti), при t=0.1
xi 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t11 0.1 0.1153 0.2192 0.3016 0.3544 0.3726 0.3544 0.3016 0.2192 0.1153
/>
Рисунок 5 –Решениеu=u(хi, ti), для метода Кранка-Николсона
ЗАКЛЮЧЕНИЕ
В зависимости от формы области, краевых условий,коэффициентов исходного уравнения метод конечных разностей имеет погрешностиаппроксимации от первого до четвертого порядка относительно шага. В силу этогоони успешно используются для разработки программных комплексов автоматизированногопроектирования технических объектов.
В МКР строятся, как правило, регулярные сетки,особенности геометрии области учитываются только в около граничных узлах. Всвязи с этим МКР чаще применяется для анализа задач с прямолинейными границамиобластей определения функций.
Проблемой методов конечных разностей является высокаяразмерность результирующей системы алгебраических уравнений (несколько десятковтысяч в реальных задачах. Поэтому реализация методов конечных разностей в составеСАПР требует разработки специальных способов хранения матрицы коэффициентовсистемы и методов решения последней.
Библиографическийсписок
1 Березин И.С.,Жидков Н.П. Методы вычислений. Т.2. – М.: Физматгиз, 1962.
2 Мэтьюз,Джон, Г., Финк, Куртис, Д. Численные методы. Использование MATLAB, 3-еиздание.— М.: Вильяме, 2001. — 720 с
3 Тихонов А.Н.,Самарский А.А. Уравнения математической физики. – М.: Наука, 1972.
4 ФормалевВ.Ф., Ревизников Д.Л. Численные методы. – М.: ФИЗМАТЛИТ, 2004. — 400 с.
5 Пирумов У.Г.Численные методы. – М.: Издательство МАИ, 1998.
6 Калиткин Н.Н.Численные методы. – М.: Наука, 1976.
ПРИЛОЖЕНИЕ А
Листингпрограммы для расчета по методу Кранка-Николсона
f.m
function F=f(x)
F=sin(pi*x)+sin(3*pi*x);
ft.m
function F=f(x,t)
F=sin(pi*x)*exp(-pi^2*t)+sin(3*pi*x)*exp(-9*pi^2*t);
tisys.m
function X=trisys(A,D,C,B)
N=length(B);
for k=2:N
mult=A(k-1)/D(k-1);
D(k)=D(k)-mult*C(k-1);
B(k)=B(k)-mult*B(k-1);
end
X(N)=B(N)/D(N);
for k= N-1:-1:1
X(k)=(B(k)-C(k)*X(k+1))/D(k);
end
crnich.m
function[U,Y]=crnich(c1,c2,a,b,c,n,m)
clc;
% — c1=u(0,t) иc2=u(a,t)
% — а и b — правые точки интервалов [0, а]и [0, Ь]
% — с — постоянная уравнениятеплопроводности
% — n и m — число точек решетки наинтервалах [0, а] и [0, Ь]
%Выход — U — матрица решений
%Инициализация параметров и матрицы U
h=a/(n-1);
k=b/(m-1);
r=c^2*k/h^2;
s1=2+2/r;
s2=2/r-2;
U=zeros(n,m);
%Граничные условия
U(1,1:m)=c1;
Продолжение ПРИЛОЖЕНИЯ А
U(n,1:m)=c2;
%Генерирование первого ряда
U(2:n-1,1)=f(h:h:(n-2)*h)';
%Формирование диагональных и не лежащих надиагонали
%элементов А, вектора постоянных В '
%и решение трехдиагональной системы АХ=В
Vd(1,1:n)=s1*ones(1,n);
Vd(1)=1;
Vd(n)=1;
Va=-ones(1,n-1);
Va(n-1)=0;
Vc=-ones(1,n-1);
Vc(1)=0;
Vb(1)=c1;
Vb(n)=c2;
for j=2:m
for i=2:n-1
Vb(i)=U(i-1,j-1)+U(i+1,j-1)+s2*U(i,j-1);
end
X=trisys(Va,Vd,Vc,Vb);
U(1:n,j)=X';
end
U=U';
%точное решение и определение погрешности
x=0:h:h*(n-1);
t=0:k:k*(m-1);
for i=1:1:n
for j=1:1:m
Y(i,j)=ft(x(i),t(j));
end
end
Y=Y';
pogr=(abs(Y-U));
max(max(pogr))
surf(U)
xlabel('X');
ylabel('t');
zlabel('U(x,t)');
colorbar;
axis([0 n 0 m 0 max(max(U))]);