close

Вход

Забыли?

вход по аккаунту

?

31.Решение краевых задач для дифференциальных уравнений в частных производных методом конечных ра.

код для вставкиСкачать
ТОМСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ СИСТЕМ
УПРАВЛЕНИЯ И РАДИОЭЛЕКТРОНИКИ
( ТУСУР )
Кафедра сверхвысокочастотной и квантовой радиотехники
( СВЧ и КР )
РЕШЕНИЕ КРАЕВЫХ ЗАДАЧ ДЛЯ
ДИФФЕРЕНЦИАЛЬНЫХ
УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ
МЕТОДОМ КОНЕЧНЫХ РАЗНОСТЕЙ
2011
Министерство образования и науки РФ
ТОМСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ СИСТЕМ
УПРАВЛЕНИЯ И РАДИОЭЛЕКТРОНИКИ
( ТУСУР )
Кафедра сверхвысокочастотной и квантовой радиотехники
( СВЧ и КР )
УТВЕРЖДАЮ
Заведующий кафедрой
__________С. Н. Шарангович
Методы математической физики
РЕШЕНИЕ КРАЕВЫХ ЗАДАЧ ДЛЯ ДИФФЕРЕНЦИАЛЬНЫХ
УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ
МЕТОДОМ КОНЕЧНЫХ РАЗНОСТЕЙ
Методические указания к компьютерной лабораторной работе
для студентов специальности
210401 – Физика и техника оптической связи
Разработчик
профессор кафедры СВЧ и КР
_____________ Г. Г. Гошин
Программное обеспечение
студент гр. 43-в Д. А. Заев
2011
Содержание
1. Введение…………………………………………………………………... 3
2. Цель работы………………………………………………………………... 3
3. Основные теоретические сведения………………………………………. 3
3.1. Решение задачи Дирихле для уравнения Лапласа…………………. 3
3.2. Решение задачи для уравнения параболического типа……………. 5
3.3. Решение задачи для уравнения гиперболического типа…………... 6
3.4. Вид аналитических решений………………………………………… 8
4. Описание интерфейса и методика работы с программой……………….. 8
5. Задания и рекомендации по выполнению работы……………………….. 9
5.1. Уравнение гиперболического типа…………………………………. 10
5.2. Уравнение параболического типа…………………………………... 10
5.3. Уравнение эллиптического типа……………………………………. 11
6. Рекомендации по оформлению отчета…………………………………… 11
7. Контрольные вопросы……………………………………………………… 12
8. Рекомендуемая литература………………………………………………… 12
1. Введение
Решение задач математической физики на ЭВМ обычно сводится к
составлению систем алгебраических уравнений той или иной структуры и их
последующему решению [1] – [6]. При построении алгебраических уравнений
необходимо знать, к какому классу функций принадлежит решение задачи,
свойства операторов задачи, свойства входных данных и т.д. Свойства
операторов дифференциальных уравнений должны быть, по возможности,
сохранены при переходе от функций непрерывных аргументов к дискретным.
Именно так, например, осуществляется переход от уравнений с частными
производными к системе алгебраических уравнений в методе конечных
разностей.
В лабораторной работе рассматриваются уравнения математической физики
в частных производных эллиптического, параболического и гиперболического
типов, которые решаются численно с помощью метода конечных разностей.
Полученные результаты в частных случаях сравниваются с расчетами,
выполненными по известным аналитическим формулам.
Впервые данная работа в ТУСУР была поставлена в 2000г. и тогда же были
изданы к ней методические указания. В новое издание методических указаний
внесены исправления и
уточнения, учитывающие накопленный опыт
проведения лабораторной работы.
2. Цель работы
Целью работы является: закрепление лекционного материала по краевым
задачам математической физики, ознакомление с методом конечных разностей
численного решения этих задач, а также развитие навыков физической
интерпретации полученных результатов.
3. Основные теоретические сведения
3.1. Решение задачи Дирихле для уравнения Лапласа
Необходимо найти непрерывную функцию U ( x, y ) , удовлетворяющую
внутри прямоугольной области Ω = {( x, y )  0 ≤ x ≤ a , 0 ≤ y ≤ b } уравнению
Лапласа
∂ 2U ∂ 2U
∆U = 2 + 2 = 0
(1)
∂x
∂y
и принимающую на границе области Ω заданные значения:
U (0, y ) = f1 ( y ) ≡ 0 , y ∈ [0, b] ; U ( a , y ) = f 2 ( y ) ≡ 0 , y ∈ [0; b] ;
x ∈[0; a] ,
(2)
U ( x,0) = f 3 ( x ) ≡ 0 , x ∈ [0, a] ; U ( x, b ) = f 4 ( x ) ,
где
U = const , x ∈ [ x1 , x2 ],
f 4 ( x) = 
0, x ∉ [ x1 , x2 ]⋅
Функция U ( x, y ) описывает распределение электрического потенциала в
прямоугольной области, являющейся поперечным сечением прямоугольной
полости, стенки которой кроме верхней заземлены, а на верхней стенке
расположен изолированный электрод с приложенным к нему постоянным
напряжением.
Выбрав по обоим направлениям одинаковый шаг h , строим сетку xi = i ⋅ h ,
i = 1, 2, ..., n ; y j = j ⋅ h , j = 1, 2, ..., m . Вводя обозначения U i , j = U ( xi , y j ) ,
аппроксимируем частные производные в каждом внутреннем узле сетки
центральными разностными производными второго порядка и заменим
уравнение Лапласа системой конечно-разностных уравнений
1
U i , j = (U i +1, j + U i −1, j + U i , j +1 + U i , j −1 ) .
(3)
4
Погрешность замены дифференциального уравнения разностным составляет
величину 0( h 2 ). При получении сеточных уравнений (3) была использована
схема узлов, изображенная на рис.1 (шаблон типа “крест”). Система называется
Рис. 1. Четырехточечная схема для уравнения Лапласа
неявной, так как не позволяет в явном виде выразить U i , j через значения
функции на предыдущих слоях.
Численное решение задачи Дирихле состоит в нахождении приближенных
значений U i , j искомой функции U ( x, y ) во внутренних узлах сетки. Для этого
нужно решить систему линейных алгебраических уравнений (3). В данном
случае система решается итерационным методом Гаусса-Зейделя [4], который
состоит в построении последовательности вида:
U is, +j 1 =
1 s +1
U i −1, j + U is+1, j + U is, j +1 + U is, j −1 ) ,
(
4
(4)
где
верхним
индексом
s
обозначен
номер
итерации.
При
s→∞
s
i, j
последовательность U
сходится к точному решению системы (3). В
качестве условия окончания итерационного процесса можно принять
max (Uis, j − Uis,+j 1 ) / Uis, j < ε , 1 ≤ i ≤ n − 1 , 1 ≤ j ≤ m − 1,
s
(5)
s +1
где ε – заданная величина, на которую могут отличаться U i , j и U i , j .
Таким образом, погрешность решения численным методом задачи (1)
складывается из двух погрешностей: погрешности аппроксимации
дифференциального
уравнения
разностным
и
погрешности
из-за
приближенного решения системы (3). Полученная разностная схема обладает
устойчивостью, т.е. малые отклонения начальных данных приводят к малым
изменениям решения задачи.
3.2. Решение задачи для уравнения параболического типа
Рассмотрим задачу для однородного уравнения теплопроводности. Задача
состоит в отыскании функции
U ( x, t ) , удовлетворяющей в области
D={ ( x, t ) │ 0 ≤ x ≤ L , 0 ≤ t ≤ T } уравнению
начальному условию
∂U
∂ 2U
− K 2 2 = 0 , K = const > 0 ,
∂t
∂x
U ( x,0) = f ( x)
и граничным условиям первого рода
U (0, t ) = µ1 (t ) , U ( L, t ) = µ 2 (t ) ,
(6)
(7)
(8)
где K = k cρ – коэффициент температуропроводности,
k – коэффициент внутренней теплопроводности,
c – удельная теплоёмкость,
ρ – плотность стержня.
К задаче (6) – (8) приводит, в частности, задача о распределении тепла в
однородном стержне длиной L, на торцах которого заданы некоторые значения
температуры; боковая поверхность стержня считается теплоизолированной.
Построим в области D равномерную прямоугольную сетку с шагом h в
направлении оси x и шагом τ в направлении оси t (рис. 2).
Рис. 2. Четырехточечная схема для уравнения теплопроводности
Обозначим узлы сетки через ( x1 , t1 ),...( xi , t j )... , а приближенные значения
функции U ( x, t ) в этих узлах как U i , j . Тогда
xi = i ⋅ h , i = 0, 1, 2,..., n , h = L / n ;
t j = j ⋅ τ , j = 0,1, 2,..., m , τ = T / m .
Аппроксимируем уравнение (6) на четырехточечном шаблоне, который
изображен на рис. 2 жирными линиями. В результате получим неявную схему
λU i +1, j − (1 + 2λ )U i , j + λU i −1, j = −U i , j −1 ,
(9)
где λ = τ K 2 / h 2 .
2
Разностная схема (9) аппроксимирует уравнение (6) с погрешностью 0(τ + h ) .
Число уравнений в (9) меньше числа неизвестных U i , j . Недостающие
уравнения получают из краевых условий
U 0, j = µ1 (t j ) , U k , j = µ 2 (t j ) .
(10)
Для решения системы линейных уравнений (9) рациональнее всего
применить метод прогонки [4]. Алгоритм численного решения задачи имеет
следующий вид: на нулевом слое ( j = 0) решение известно из начального
условия U i ,0 = f ( xi ) ; на каждом последующем слое искомая функция
определяется как решение системы (9), (10). Разностная схема (9) устойчива
при любых значениях параметра
λ = τ K 2 / h2 > 0 .
(11)
Схема обладает сходимостью. Это означает, что при h → 0 и τ → 0 решение
дискретной задачи стремится к точному решению непрерывной задачи (6) – (8).
3.3. Решение задачи для уравнения гиперболического типа
Рассмотрим задачу для однородного уравнения колебаний струны.
Необходимо найти функцию U ( x, t ) , описывающую смещение точек струны и
удовлетворяющую при τ > 0 уравнению
2
∂ 2U
2 ∂ U
−C
= 0, 0 ≤ x ≤ L , 0 ≤ t ≤T ,
∂t 2
∂x 2
начальным условиям
∂U
U ( x,0) = f ( x) ,
( x, t ) = V ( x) при t = 0 , 0 ≤ x ≤ L ,
∂t
и граничным условиям
U (0, t ) = µ1 (t ) , U ( L, t ) = µ2 (t ) , 0 ≤ t ≤ T .
где C =
T0
ρ
(12)
(13)
(14)
– коэффициент, характеризующий скорость распространения
колебательного процесса;
T0 =const – сила натяжения струны;
ρ =const – плотность струны.
В (13) функции f ( x) и V ( x) описывают начальное воздействие на струну: f ( x)
определяет начальное положение точек струны относительно оси x , а V ( x) –
начальную скорость смещения точек струны, например в результате удара.
Для построения разностной схемы задачи (12) – (14) зададим в области
D = { ( x, t ) │ 0 ≤ x ≤ L , 0 ≤ t ≤ T } сетку
xi = i ⋅ h , i = 0,1, 2,..., n + 1 ; L = h ⋅ ( n + 1) ;
t j = j ⋅ τ , j = 0,1,2,..., m ; T = m ⋅ τ
и аппроксимируем уравнение (12) в каждом внутреннем узле на шаблоне типа
“крест” (рис. 3).
Рис. 3. Четырехточечная схема для уравнения колебаний струны
Используя для аппроксимации частных производных
разностные производные, приходим к уравнению
центральные
U i , j +1 − 2U i , j + U i , j −1
= C2
U i +1, j − 2U i , j + U i −1, j
,
(15)
h2
где U i , j – приближенное значение функции U ( x, t ) в узле ( xi , t j ) . Полагая
λ = τ ⋅ C / h , получаем трехслойную разностную схему
U i , j +1 = 2(1 − λ2 )U i , j + λ2 (U i +1, j − U i −1, j ) − U i , j −1 , i = 1, 2,..., n .
(16)
В данной работе выбирем нулевые граничные условия, т.е. считаем, что концы
струны закреплены и тогда µ1 (t ) ≡ 0 , µ 2 (t ) ≡ 0 . Это означает, что в (16)
U 0, j = 0 , U k , j = 0 для всех j . Разностная схема (16) явная, т.е. позволяет в
явном виде выразить U i , j через значения функции на предыдущих слоях.
Описанная схема аппроксимирует непрерывную задачу (12) – (14) с точностью
0(τ + h 2 ). Схема устойчива, если выполнено условие Куранта [4]
τ ⋅C ≤ h.
(17)
При выполнении этого условия малые погрешности, возникающие,
например, при вычислении решения на первом слое, не будут неограниченно
возрастать при переходе к каждому новому слою. Если выполняется условие
Куранта, то схема обладает также равномерной сходимостью, т.е. при τ → 0 ,
h → 0 решение разностной задачи равномерно стремится к решению задачи
(12) – (14). Как видно из (12) и (17), C имеет размерность скорости, т.е.
описывает скорость распространения колебательного процесса.
τ2
3.4. Вид аналитических решений
Аналитические решения известны для весьма ограниченного круга задач.
Обычно они находятся методом разделения переменных [6].
Решение задачи для уравнения Лапласа в случае, когда электрод занимает
всю верхнюю стенку полости, имеет вид [5]
πx
πy
] ⋅ sh[(2n + 1) ]
∞ sin[(2n + 1)
4
a
a ,
U ( x, y ) = U 0
(18)
π
b
π
n =0
(2n + 1) sh[(2n + 1) ]
a
где U 0 = const – напряжение на электроде.
Решение задачи о распространении тепла в стержне длиной L, на концах
которого поддерживаются заданные постоянные значения температуры T1 и
T2 , а начальная температура самого стержня T0 = 0 , имеет вид [6]
∑
x 2
T ( x, t ) = T1 + (T2 − T1 ) +
L π
∞
( −1) n T2 − T1 − (
e
n
n =1
∑
nπ 2
) t
L
⋅ sin(
n πx
).
L
(19)
Формула записана при K = 1 .
Следует заметить, что выписанные формулы представляют собой ряды
Фурье с бесконечными верхними пределами, область сходимости которых
зависит от значений входящих в них параметров и переменных. Поэтому их
суммирование также выполняется с некоторой погрешностью.
4. Описание интерфейса и методика работы с программой
Программа имеет стандартный интерфейс – есть главное меню и “быстрые
кнопки” с подсказками. Размеры главного окна программы можно изменять,
т.е. минимизировать или раскрыть на весь экран. Главное окно содержит
“Рисунок”– компонент, в котором происходит построение всех графиков, и
четыре таблицы. Границу между рисунком и таблицами можно передвигать
влево или вправо, подведя курсор “мыши” к границе и нажав левую кнопку.
В первую таблицу заносятся результаты решения задач численным методом.
В некоторых частных случаях известны аналитические решения задач. Расчеты
по ним заносятся во вторую таблицу. Относительная погрешность,
вычисленная по формуле
(20)
( fa − fч ) / fa ,
заносится программой в третью таблицу, которая прежде должна быть
активизирована
соответствующей
кнопкой.
Четвертая
таблица
вспомогательная и служит для накопления любых данных с последующей
возможностью построения по ним графика. Эти данные могут браться из
предыдущих таблиц. В (20) индексы указывают на принадлежность к
аналитическому или численному решениям.
Если какая-либо таблица заполнена, то для того, чтобы построить график,
нужно выбрать соответствующую строку в таблице параметров задачи,
подвести курсор “мыши” к первой колонке и нажать левую кнопку. При этом в
список параметров графика (находится справа от координатной сетки)
помещается содержимое ячейки первой колонки для выбранной строки и
строится график. Для быстрого просмотра графиков предназначена крайняя
правая кнопка. Кнопка “Рисунок” дает возможность построения трехмерных
(3D) графиков с возможностью выбора нужного ракурса. Если курсор “мыши”
поставить на график, нажать левую кнопку и повести вниз, то можно
просмотреть выделенную часть графика в увеличенном виде. Возвращение в
исходное состояние – нажать левую кнопку и повести вверх.
Для того, чтобы очистить рисунок, нужно выбрать в главном меню
“Вид/Очистить рисунок” или воспользоваться соответствующей кнопкой.
Построенный график можно вывести на печать (команда главного меню
“Вид/Вывести рисунок на печать”) или скопировать в буфер (“Вид/Копировать
в буфер”) и затем вставить его в текст в редакторе Word. Для выполнения
команд также можно воспользоваться “быстрыми кнопками“.
5. Задания и рекомендации по выполнению работы
Варианты задания для каждого типа уравнений представлены в таблицах. В
них содержатся значения параметров, входящих в уравнения, в граничные и
начальные условия, а также размеры шагов сетки. Далее следуют сами задания.
Для того чтобы выполнить вариант задания, нужно:
• выбрать в главном меню ’’Уравнение’’ и соответствующий вид уравнения
или воспользоваться кнопками;
• в появившемся окне задать необходимые параметры, соответствующие
варианту задания; если указаний нет, то оставить по умолчанию;
• установить переключатели “Выполнить аналитическое решение”,
“Выполнить численное решение” в нужные положения и не требовать
“Выполнить аналитическое решение”, если оно отсутствует;
• нажать кнопку “ОК”, если вы хотите провести расчет для введенных
параметров. Нажать кнопку “Отмена”, если нужно закрыть окно и не
проводить вычисления;
• для построения графиков нажать в соответствующих таблицах кнопки с
выбранными значениями параметров. С целью отбора графиков и
получения информации о динамике процессов рекомендуется
пользоваться кнопкой быстрого просмотра. Окончательное формирование
графиков для представления в отчет следует проводить с использованием
3D-графики;
• выбрать в меню “Таблица/Вычислить относительную погрешность” или
нажать соответствующую кнопку, если необходимо вычислить
погрешность;
• повторить все предыдущие пункты для других типов уравнений.
Параметры
5.1. Уравнение гиперболического типа
Таблица 1.
Варианты
L, м
1
4,0
2
4,5
3
5,0
4
5,5
5
6,0
6
3,5
7
3,0
8
2,5
9
2,0
10
1,5
hx, м
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,05
h t, с
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,05
Решение уравнения выполнять численным методом.
1. Выбрать T = 5; C = 1; А = 20. Последовательно задавая n = 1, 2, 3, описать
процесс распространения волн отклонения вдоль струны.
2. Выбрать А = 20; C = 1; (x2 - x1) = 0,1)L. Последовательно задавая x1 = (0;
0,1; 0,5)L, описать процесс распространения вдоль струны образующейся в
результате удара волны импульса, в зависимости от трех указанных положений
интервала воздействия (места удара). Время T выбрать достаточным для
наблюдения и понимания физики процесса (T = 10…20).
3. В случае волн отклонения при n = 2 выявить зависимость амплитуды
колебательного процесса от параметра C, принимающего значения C = 0,2; 0,5;
1,0; 1,2. Определить критические значения Cкр > C , начиная с которых
наблюдаются “аномалии”. Объяснить их.
5.2. Уравнение параболического типа
Таблица 2.
Варианты
Параметры
1
2
3
4
5
6
7
8
9
10
L, м
2
2
2
2
2
2
1
1
1
1
T1,˚С
10
10
10
-10
-10
-10
10
10
5
-5
T2,˚С
-10
-5
0
10
5
0
5
-10
10
-10
1. Для K = 1 и разных моментов времени на одном рисунке построить
графики распределения температуры вдоль стержня, включая установившийся
режим.
2. Задать hx = 0,1. Оценить время установления ожидаемого линейного
распределения температуры в стержне в зависимости от параметра K,
принимающего значения K = 0,2; 0,6; 1,0; 1,5; 2,0. Построить график.
3. Построить и сравнить графики, полученные численным методом, и
выполненными программой расчетами по формуле (19), описывающей
аналитическое решение при K = 1 и T0 = 0˚С. Сделать выводы о характере
зависимостей. Объяснить их поведение, анализируя распределения
температуры для численного и аналитического решений.
5.3. Уравнение эллиптического типа
Варианты
Параметры
Таблица 3.
1
2
3
4
5
6
7
8
9
10
а, м
1
1,5
2
1
1
2
1,5
2
1
1
b, м
1
1
1
1,5
2
2
1
1
1,5
2
U, В
10
10
10
10
10
50
50
50
50
50
1. В окне метода Гаусса-Зейделя установить погрешность 10-3 .
2. Построить и сравнить графики, полученные численным методом, и
выполненными программой расчетами по формуле (18), описывающей
аналитическое решение при x1 = 0, x2 = а. Сделать выводы о характере
зависимостей. Объяснить их поведение, анализируя распределения потенциала
для численного и аналитического решений.
3. В случае численного решения при заданных в таблице параметрах для
двух случаев hх = 0,05, x1 = 0,1а, x2 = 0,2а и hх = 0,05, x1 = 0,5а, x2 = 0,8а
занести в “Рисунок” пошаговые графики распределения потенциала.
Объяснить получающуюся асимметрию распределения.
6. Рекомендации по оформлению отчета
Отчет по лабораторной работе должен иметь титульный лист и содержать
следующие данные:
• номера вариантов задания и выбранные параметры;
• результаты решения задач в графическом виде, характеризующие
рассматриваемые процессы;
• сравнение в частных случаях численного решения с решением
аналитическим;
• пример неустойчивого решения уравнения гиперболического типа;
• физическая интерпретация результатов с описанием влияния параметров на
характеристики процессов;
• выводы по существу проделанной работы, а не перечисление выполненных
операций.
Отчет должен быть представлен в виде распечатки и текстового файла.
7. Контрольные вопросы
1. Запишите основные типы уравнений в частных производных второго
порядка и назовите процессы, описываемые этими уравнениями.
2. В чем заключается постановка краевых задач для дифференциальных
уравнений в частных производных второго порядка?
3. Каким образом происходит переход от функций непрерывных аргументов
к функциям дискретных аргументов?
4. Дайте определение явной и неявной разностных схем.
5. Дайте определение устойчивой и неустойчивой разностных схем.
Запишите условие Куранта.
6. Ваши суждения по поводу ожидаемого поведения решений, описывающих
исследуемые процессы: распределение электрического потенциала в
прямоугольной области, установление температурного режима в стержне,
колебания струн при двух возможных видах воздействия.
8. Рекомендуемая литература
1. Вержбицский В.М. Основы численных методов. – М.: ВШ, 2002. – 840 с.
2. Самарский А.А. Введение в численные методы. – СПб.: Лань, 2005. – 288с.
3. Самарский А.А., Гулин А.В. Численные методы. – М.: Наука, 1989. – 432с.
4. Плис А.И., Сливина Н.А. Лабораторный практикум по высшей
математике. – М.: ВШ, 1983. – 208с.
5. Мудров А.Е. Численные методы для ПВМ. – Томск: Раско, 1991. – 272с.
6. Кошляков Н.С., Глинер Э.Б. и др. Уравнения математической физики в
частных производных. – М: ВШ, 1970. – 712с.
Документ
Категория
Без категории
Просмотров
16
Размер файла
182 Кб
Теги
решение, методов, уравнения, конечный, дифференциальной, частных, производной, задачи, краевых
1/--страниц
Пожаловаться на содержимое документа