Численное решение уравнений математической физики в интегрированной среде Mathcad

38 Pages·0.678 MB·Russian

by  Михайлова Н.А. (сост.)

Министерство образования и науки Российской Федерации
Волгоградский государственный архитектурно-строительный
университет

Численное решение уравнений
математической физики
в интегрированной среде Маthcad

Методические указания к лабораторным работам
по дисциплине «Специальные разделы высшей математики»

Составитель Н.А. Михайлова

© Федеральное государственное бюджетное
образовательное учреждение
высшего профессионального образования
«Волгоградский государственный архитектурно-
строительный университет», 2012

Волгоград
ВолгГАСУ
2012

1
УДК 517.5:004.42(076.5)

Численное решение уравнений математической физики в интегриро-
ванной среде Mathcad [Электронный ресурс] : методические указания к лабо-
раторным работам по дисциплине «Специальные разделы высшей математи-
ки» / М-во образования и науки Росс. Федерации, Волгогр. гос. архит.-
строит. ун-т ; сост. Н.А. Михайлова. – Электронные текстовые и графические
данные (0,5 Мбайт). – Волгоград : ВолгГАСУ, 2012. – Учебное электронное
издание комбинированного распространения : 1 СD-диск. – Систем. требова-
ния: PC 486 DX-33; Microsoft Windows XP; 2-скоростной дисковод СD-ROM;
Adobe Reader 6.0. – Официальный сайт Волгоградского государственного ар-
хитектурно-строительного университета. Режим доступа:
http://www.vgasu.ru/publishing/on-line/ – Загл. с титул. экрана.

Содержатся краткие теоретические сведения, необходимые для выполнения лабора-
торных работ по дисциплине «Специальные разделы высшей математики», приведены ва-
рианты индивидуальных заданий для лабораторной работы, сформулированы контроль-
ные вопросы по изучаемой теме.
Для магистров направления «Строительство» очной формы обучения.
УДК 517.5:004.42(076.5)

Нелегальное использование продукта запрещено

2
ОГЛАВЛЕНИЕ
1. ТЕОРЕТИЧЕСКАЯ ЧАСТЬ……..……………………………………… 4
1.1.Введение………………………………….. 4
1.2. Сетки и сеточные функции ………………………………………... 6
1.3. Аппроксимация простейших дифференциальных операто-
ров………………………………………………………………………… 7
1.4. Разностная задача…………………………………………………… 11
1.5. Устойчивость…………………………………………………….. 12
2 ЛАБОРАТОРНЫЕ РАБОТЫ…………………………............................ 15
2.1. Лабораторная работа 1. Решение первой краевой задачи для ги-
перболического уравнения методом конечных разностей в интегри-
рованной среде Mathcad……..................................................................... 15
2.2. Лабораторная работа 2. Решение первой краевой задачи для па-
раболического уравнения методом конечных разностей в интегриро-
ванной среде Mathcad……...............................................................…… 21
2.3. Лабораторная работа 3. Решение задачи Дирихле для уравнения
Лапласа методом конечных разностей в интегрированной среде
Mathcad…….............................………………………………………….... 26
2.4. Содержание отчета…………………………………………………... 34
2.5. Список рекомендуемой литературы…………………………..……. 34
2.6. Приложение…………………………………….................................. 35

3
Цель лабораторных работ - изучение численных методов решения урав-
нений математической физики, выработка навыков применения этих методов
в решении конкретных задач, умения оценить полученные результаты и вы-
брать наилучший метод.

1. ТЕОРЕТИЧЕСКАЯ ЧАСТЬ

1.1. Введение
Многие задачи математической физики приводят к дифференциальным
уравнениям с частными производными. Наиболее часто встречаются диффе-
ренциальные уравнения 2-го порядка.

Уравнением с частными производными 2-го порядка с двумя независи-
мыми переменными х, у называется соотношение между неизвестной функ-
цией и(х, у) и ее частными производными до 2-го порядка включительно:
F(x, y,u,u ,u ,u ,u ,u ) = 0.
x y xx xy yy
Мы пользуемся следующими обозначениями для производных:
∂u ∂u ∂2u ∂2u ∂2u
u = , u = , u = , u = , u = .
x ∂x y ∂y xx ∂x2 xy ∂x∂y yy ∂y2
Аналогично записывается уравнение и для большего числа независимых
переменных.
Уравнение называется линейным относительно старших производных,
если оно имеет вид
a u + 2a u + a u + F (x, y,u,u ,u ) = 0, (1.1)
11 xx 12 xy 22 yy 1 x y
где a , a , a являются функциями х и у.
11 12 22
Если коэффициенты a , a , a зависят не только от х и у, а являются,
11 12 22
подобно F , функциями х, у, и, и , и , то такое уравнение называется квази-
1 х y
линейным.
Уравнение называется линейным, если оно линейно как относительно
старших производных и , u , и , так и относительно функции u и ее первых
хх ху уу
производных и , и :
х y
a u + 2a u + a u +bu +b u +cu + f = 0,
11 xx 12 xy 22 yy 1 x 2 y
где a , a , a , b , b , c, f — функции только х и у. Если коэффициенты –
1 2
11 12 22
уравнения (2) не зависят от х и у, то оно представляет собой линейное урав-
нение с постоянными коэффициентами. Уравнение называется однородным,
если f(х, у) = 0.
В зависимости от знака выражения a2 −a a уравнение (1.1) называет-
12 11 22
ся в точке M области D(x,y) уравнением:
гиперболического типа, если в точке M a2 −a a > 0,
12 11 22
эллиптического типа, если в точке M a2 −a a < 0,
12 11 22
параболического типа, если в точке M a2 −a a = 0 .
12 11 22
4
Уравнения с частными производными 2-го порядка гиперболического ти-
па наиболее часто встречаются в физических задачах, связанных с процесса-
ми колебаний. Простейшее уравнение гиперболического типа
u = a2u или u −u = 0 , (y = at)
tt xx xx yy
обычно называют уравнением колебаний струны. К гиперболическим
уравнениям приводят уравнение электрических колебаний в проводах, попе-
речные колебания мембраны, уравнения гидродинамики и акустики.
Уравнения с частными производными 2-го порядка параболического типа
наиболее часто встречаются при изучении процессов теплопроводности и
диффузии. Простейшее уравнение параболического типа
u = a2u или u −u = 0, (y = a2t),
t xx xx y
обычно называют уравнением теплопроводности.
При исследовании стационарных процессов различной физической при-
роды (колебания, теплопроводность, диффузия, и др. ) обычно приходят к
уравнениям эллиптического типа. Наиболее распространенным уравнением
этого типа является уравнение Лапласа
∆u = 0.
При математическом описании физического процесса прежде всего надо
поставить задачу, т. е. сформулировать условия, достаточные для однознач-
ного определения процесса.
Дифференциальные уравнения с обыкновенными и тем более с частными
производными имеют, вообще говоря, бесчисленное множество решений.
Поэтому в том случае, когда физическая задача приводится к уравнению с
частными производными, для однозначной характеристики процесса необхо-
димо к уравнению присоединить некоторые дополнительные условия.
Рассмотрим простейшую задачу о поперечных колебаниях струны, закре-
пленной на концах. В этой задаче и(х,t) даст отклонение струны от оси х. Ес-
ли концы струны x = 0 и x = l закреплены, то должны выполняться граничные
условия
и(0,t) = 0, и(l,t) = 0. (1.2)
Так как процесс колебаний струны зависит от ее начальной формы и рас-
пределения скоростей, то следует задать начальные условия:
( ) ( )
u x,t =ϕx ,
0
( ) ( )
u x,t =φx .
t 0
Таким образом, дополнительные условия состоят из граничных и началь-
( ) ( )
ных условий, где ϕx и φx - заданные функции точки. Если концы струны
движутся по заданному закону, то граничные условия (1.2) принимают вид
и(0,t) = µ (t), и(l,t) = µ (t),
1 2
где µ (t) и µ (t) - заданные функции времени t.Это так называемые гра-
1 2
ничные условия первого рода.
Всякий физически определенный процесс, развивающийся во времени,
должен характеризоваться функциями, непрерывно зависящими от началь-
5
ных данных. Если решение математической задачи непрерывно зависит от
исходных данных, то говорят, что задача устойчива. В связи с изучением фи-
зически детерминированных явлений вводится понятие корректности. Го-
ворят, что математическая задача поставлена корректно, если решение зада-
чи существует, задача имеет единственное решение, решение задачи непре-
рывно зависит от исходных данных (устойчиво).
Существует множество аналитических методов решения уравнений с ча-
стными производными, таких, как, например, метод разделения переменных
или использование представление решения в виде интегралов, называемых
потенциалами. Однако явное представление решения в виде ряда или инте-
грала не всегда возможно. Представление решения нелинейных уравнений в
аналитической форме возможно лишь в исключительных случаях.
Универсальным методом приближенного решения дифференциальных
уравнений, применимым для широкого класса уравнений математической
физики, является метод конечных разностей (или метод сеток).
Метод конечных разностей состоит в следующем. Область непрерывного
изменения аргументов (например, х и t) заменяется конечным (дискретным)
множеством точек (узлов), называемым сеткой; вместо функций непрерыв-
ного аргумента рассматриваются функции дискретного аргумента, опреде-
ленные в узлах сетки и называемые сеточными функциями. Производные,
входящие в дифференциальное уравнение, заменяются (аппроксимируются)
при помощи соответствующих разностных отношений; дифференциальное
уравнение при этом заменяется системой алгебраических уравнений (разно-
стным уравнением). Начальные и граничные условия тоже заменяются на-
чальными и граничными условиями для сеточной функции.
Естественно требовать, чтобы полученная таким образом разностная
краевая задача была разрешима и ее решение при увеличении числа N узлов
сетки приближалось (сходилось) к решению исходной задачи для дифферен-
циального уравнения. Введем понятия аппроксимации, сходимости, точности
и устойчивости и проиллюстрируем их на простейших примерах.

1.2. Cетки и сеточные функции

Рассмотрим простейшие примеры сеток. Пусть область изменения аргу-
мента х есть отрезок 0≤ x ≤l. Разобьем отрезок 0≤ x ≤lточками x =ih (i=0,
i
1, 2, ..., N; h > 0) на N равных частей длины h=l/N каждая. Множество точек
x =ih, i=0, 1, 2, ..., N называется разностной сеткой на отрезке0≤ x ≤lи обо-
i
{ }
значается ω = x ,i = 0,1,...,N , а число h — расстояние между точками (уз-
h i
лами) сетки — называется шагом сетки.
Отрезок [0, l] можно разбить на N частей, вводя произвольные точки
{ }
x < x <...< x <l. Тогда получим сетку ω = x ,i = 0,1,...,N,x = 0,x =l с
1 2 N−1 h i 0 N
шагом h = x − x , который зависит от номера i узла х. Если h ≠ h хотя бы
i
i i i−1 i i+1
для одного номера I, то сетка ω называется неравномерной. Если h = const =
i
h
6
h =l/N для всех i = 1, 2, ..., N, то мы получаем построенную выше равномер-
ную сетку. На бесконечной прямой−∞ < x < ∞можно рассматривать сетку
{ }
ω = x+ih,i = 0,±1,±2,... с началом в любой точке х, состоящую из беско-
h,x
нечного числа узлов.
Функцию y= у(x) дискретного аргумента x, I = 0, 1, ..., N, называют се-
i i i
точной функцией, определенной на сетке ω .
h
Всякой непрерывной функции f(х) можно поставить в соответствие сеточ-
ную функцию f h, полагая, например, f h = f (x ). Впрочем, в некоторых слу-
i i i
чаях удобнее устанавливать это соответствие другими способами.
Пусть область изменения аргументов (х, t) есть прямоугольник
( )
Д = 0≤ x ≤1,0≤t ≤T, . Построим на отрезке 0≤ x ≤1сетку
{ } { }
ω = x =ih,i = 0,1,...,N с шагом h =1/N и сетку ω = τ = jτ,i = 0,1,...,N с
h i τ j 0
( )
шагом τ =T/N на отрезке 0≤t ≤T . Множество узлов x ,t с координатами
0
i j
x =ih иτ = jτ назовем сеткой в прямоугольнике Д и обозначим
i j
{( ) }
ω = x +ih,τ + jτ ,i = 0,1,...,N, j =1,2,...,N . Эта сетка равномерна по каж-
h,τ i j 0
дому из переменных х и t. Если хотя бы одна из сеток ω или ω неравно-
h τ
мерна, то сетка ω называется неравномерной. Сетка ω очевидно, состоит
h,τ h,τ
из точек пересечения прямых х = х ,i=0, 1, ..., N и прямых t = t, j= 0, 1, ..., N .
i j 0
Пусть y — сеточная функция, заданная на ω . Будем обозначать
h,τ
yj = y(x ,t )значение сеточной функции у в узле (х, t) сетки ω .
i j
i i j h,τ
Непрерывной функции и(х,t), где (х,t)—точка из Д , будем ставить в соот-
( )
ветствие сеточную функцию uj =uj =u x ,t .
i i,h,τ i j

1.3. Аппроксимация простейших дифференциальных операторов

Оператор L преобразующий сеточную функцию у в сеточную функцию
h
У = L у, называют сеточным или разностным оператором. Дифференциаль-
h
ный оператор L, заданный в классе функций непрерывного аргумента, может
быть приближенно заменен (аппроксимирован) разностным оператором L ,
h
заданным на сеточных функциях. Для этого каждая из производных заменя-
ется разностным отношением (отсюда и название «разностный оператор»),
содержащим значения сеточной функции в нескольких узлах сетки. Посмот-
рим, как это делается для первых и вторых производных функции одного пе-
ременного.
{ }
Пусть ω = x =ih — сетка с шагом h на отрезке 0≤ x ≤1. Рассмотрим
h i
первую производную Lv =v′функции v(x). Заменить ее разностным выраже-
нием можно бесчисленным множеством способов. Простейшими являются
замены
v −v
Lv ≈ i i−1 = L−v ,
h h i
7
— левая разностная производная или левое разностное отношение,
v −v
Lv ≈ i+1 i = L+v ,
h h i
— правая разностная производная или правое разностное отношение,
v −v
Lv ≈ i+1 i−1 = L0v
2h h i
- центральная разностная производная.
Здесь v=v(х), знак,≈ означает соответствие или аппроксимацию. При за-
i i
мене Lv =v′разностным выражением L±v допускается погрешность
h i
L±v −(Lv) =ψhназываемая погрешностью аппроксимации оператора L раз-
h i i i
ностным оператором L . Естественно требовать, чтобы при стремлении h к
h
нулю эта погрешность стремилась к нулю. Для оценки ψh надо предполо-
i
жить, что v(х) гладкая функция. Будем говорить, что v(х) принадлежит классу
(пространству) С(m)[0,1](v(x)∈C(m)[0,1]) функций, заданных на отрезке
0≤ x ≤1, если у(х) имеет т непрерывных на отрезке 0≤ x ≤1производных.
При т = 0 получаем класс С(0)[0,1]непрерывных при 0≤ x ≤1 функций.
Пусть v(x)∈C(m)[0,1] где m≥ 2. Разложим v(х) в окрестности точки x = x :
i
v =v ±h +O(h2)
i±1 i i
и вычислим ψh = L−v −v′ =O(h), ψh = L+v −v′ =O(h).
i h i i i h i i
Будем говорить, что разностный оператор L :
h
1. аппроксимирует дифференциальный оператор L на сетке ω , если
h
maхψh = maх L v −( Lv) ,
i h i i
ωh ωh
где v(x)—достаточно гладкая функция, стремится к нулю при h→0;
( )
2. аппроксимирует L с порядком п (n>0), если maхψh = O hn (или
i
ωh
maхψh = Mhn, где М — положительная постоянная, не зависящая от h).
i
ω
h
Обращаясь к формулам для L± , видим, что L−v и L+v аппроксимиру-
h h i h i
ютLv =v′ с первым порядком при v∈C(m)где m≥ 2. Увеличение m не меняет
порядка аппроксимации.
Выражение для L−v содержит значения v в двух узлах х =x и х =x сетки.
i i-1
h i
Говорят, что оператор L−является двухточечным или оператором первого
h
порядка.
Множество узлов, значения сеточной функции в которых входят в выра-
жениеL v , называют шаблоном оператора L в точке x. Очевидно, что шаб-
h i
h i
лон оператора L− состоит из двух узлов x и x . а шаблон L+ - из узлов x и
i i-1 i
h h
x .
i+1
Возьмем, например, трехточечный оператор, определенный на шаблоне
x , x, x :
i-1 i i-1
8
σv +(1−2σ)v −(1−σ)v
Lσv =σL+v +(1−σ)L−v = i+1 i i−1 ,
h i h i h i h
где σ - произвольное число. В частности, при σ = 1/2 получаем централь-
v −v
ную разностную производную L0v = i+1 i−1 которая, как нетрудно пока-
h i 2h
зать, при v(x)∈C(3)[0,1]аппроксимирует v′ со вторым порядком.
Введем обозначения:
v −v v −v
v = i i−1 , v = i+1 i
x,i h x,i h
v −v 1
( )
v = i+1 i−1 = v +v .
x0,i 2h 2 x,i x,i
В тех случаях, когда номер i узла не имеет значения, будем его опускать и
писать v ,v , v .
x,i x,i x0,i
Рассмотрим теперь вторую производную Lv =v′. На двухточечном шаб-
лоне, очевидно, ее аппроксимировать нельзя. Выберем трехточечный шаб-
лон, состоящий из узлов x , x , x , и рассмотрим разностный оператор
i-1 i i+1
1 v −2v +v
( )
L v = v = v −v = i+1 i i−1 .
h i xx,i h x,i x,i h2
Если v∈C(m)[0,1], m ≥ 4, то можно написать
h2 h3 h4
v = v ±hv′+ v′′± v′′+′ v(IV) +o(h4),
i±1 i i 2 i 6 i 24 i
(o(h4)- величина, стремящаяся к нулю при h→0 быстрее, чем hn). Отсюда
следует (индекс i опускаем), что
h2
v −v′′ = v(IV) +o(h2) (1.3)
xx 12 i
т. е. v аппроксимирует v′′ со вторым порядком.
xx
Для аппроксимации четвертой производной Lv = v(IV) выберем пятиточеч-
ный шаблон, состоящий из узловx +kh,(k = 0,±1,±2), положим
i
v −4v +6v −4v +v
L v = v = i−2 i−1 i i+1 i+2 (1.4)
h i xxxx,i h4
В этом случае v −v(IV) =O(h2)для v(x)∈C(6). На пятиточечном шаблоне
xxxx
x +kh,(k = 0,±1,±2) можно добиться аппроксимации О (h4) для v′′ , если
i
v(x)∈C(6). В самом деле, из (1.3) и (1.4) следует, что оператор
h2
L v =v − v =v′′+O(h4)
h xx 12 xxxx
имеет 4-й порядок аппроксимации.
На практике аппроксимация производных на многоточечных шаблонах
используется редко, так как при увеличении шаблона обычно увеличивается
9
объем вычислительной работы и ухудшаются качества получающихся раз-
ностных операторов (в смысле устойчивости).
Порядок аппроксимации разностного оператора L зависит от порядка т
h
дифференцируемости функции v(х). Выше мы фактически рассматривали
максимальный порядок аппроксимации, который не меняется при увеличе-
нии номера т класса C(m), считая, что v(х) - любая функция из C(m). Очевид-
но, что при специальном выборе v(х) порядок аппроксимации может повы-
ситься.
Рассмотрим более сложный оператор
∂u ∂2u
Lu = − ,
∂t ∂x2
гдеu =u(x,t)- функция двух аргументов х и t, меняющихся в области
( )
Д = 0≤ x ≤1,0≤t ≤T .Введем сетку
{( ) }
ω = x =ih,τ = jτ ,i = 0,1,...,N, j =1,2,...,N
h,τ i j 0
с шагами h =1/N и τ=Т/N , построенную выше. Произведем замену
0
⎛∂u⎞j+1 uj+1 −uj ⎛∂2u⎞j uj −2uj +uj
⎜ ⎟ ≈ i i =uj+1,⎜ ⎟ = i−1 i i+1 =uj
⎝ ∂t ⎠ τ t ⎝∂x2 ⎠ h2 xx,i
i i
В результате получим разностный оператор
uj+1 −uj uj −2uj +uj
L uj+1 = i i = i−1 i i+1 , (1.5)
hτ i τ h2
который можно записать в виде L u =u −uˆ , гдеuˆ =uj, u =uj+1. Этот
hτ t xx xx i i
оператор определен на шаблоне, состоящем из четырех точек (x,t ), (x,t),
i j+1 i j
(x ,t ), (x ,t ), (рис. 1,a). Оператор L определен не во всех узлах ω , а
i-1 j+1 i+1 j+1
hτ h,τ
только при 0 < i< N , j >0 т. е.

Рис.1 а) б)
в тех узлах, в которых шаблон состоит только из узлов сетки ω .Узлы (x,t)
i j
h,τ
{ }
назовем внутренними и обозначим ω = (x ,t ),i = 0,1,...,N, j =1,2,...,N
h,τ i j 0
множество всех внутренних узлов. Таким образом, оператор L определен на
hτ
ω т.е. во внутренних узлах. В остальных узлах, называемых граничными,
h,τ
должны быть заданы краевые и начальные условия. Оператор L имеет пер-
hτ
вый порядок аппроксимации по τ и второй по h:
10