МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
Инновационная образовательная программа
«Наукоемкие технологии и экономика инноваций»
Московского физико-технического института
(государственного университета)
на 2006-2007 годы
В.В. Демченко
вычислительный практикум
ПО
ПРИКЛАДНОЙ математике
КОНТРОЛЬНЫЙ листок
СРОКОВ ВОЗВРАТА
КНИГА ДОЛЖНА БЫТЬ
ВОЗВРАЩЕНА НЕ ПОЗЖЕ
УКАЗАННОГО ЗДЕСЬ СРОКА
Колич. пред, выдач.

/2	/3 &>€*
tV.ctf-r-0
Н--
МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
Федеральное агентство по образованию
Государственное образовательное учреждение
высшего профессионального образования
Московский физико-технический институт
(государственный университет)
i,?'
Kit '< i-.'Wff rjj.
B.B. Демченко
: ВЫЧИСЛИТЕЛЬНЫЙ ПРАКТИКУМ
ПО ПРИКЛАДНОЙ МАТЕМАТИКЕ
 *
.Я;
•и;
. ч.И'Х’зЯ
?
Рекомендовано	“!,1; • '
Учебно-методическим объединением
.	высших учебных заведений Российской Федерации
-	По образованию в области прикладных математики и физики
'	в качестве учебного пособия
ГОУ ЬГ;О
"Московский
физнко-техяпческий пнсткпт
иосудаг>стйл!.чы11 vxiiavvcnwi!”
БИБЛИОТЕКА
ИЙ.й
•Э^Мч	-nt •	"
МОСКВА 2007
Приложение иа CD-ROM:
«Пакет программ для
проверки лаб. работ»
удК 519-6(075)
ББК 22.19Я73
дзо
Рецензенты:
Кафедра математических и информационных технологий
Института автоматизации проектирования РАН
Доктор физико-математических наук, профессор А.И. Толстых
Демченко В.В.
ДЗО Вычислительный практикум по прикладной математике:
Учебное пособие. — М.: МФТИ, 2007. — 196 с.
ISBN 5-7417-0186-8
На примере шести актуальных задач нз физики и механики сплош-
ных сред для выбранных математических моделей рассматриваются методы
построения аналитического и численного решений, обеспечивающие задан-
ную точность результатов при расчетах на современных ЭВМ. Изучение
рекомендуемых подходов происходит при выполнении лабораторных работ
из вычислительного практикума и может контролироваться с помощью «Па-
кета программ для проверки лабораторных работ вычислительного практи-
кума», являющегося неотъемлемой частью всего учебного пособия и имею-
щего инновационный характер.
Вычислительный практикум предназначен для студентов, аспиран-
тов, преподавателей и научно-технических работников, использующих в
своей деятельности методы вычислительной и прикладной математики.
УДК 519.6(075)
ББК 22.19я73
ISBN 5-7417-0186-8
© Демченко В.В., 2007
© Московский физико-технический институт
(государственный университет), 2007
СОДЕРЖАНИЕ
Введение.......................
РАЗДЕЛ I. НЕЛИНЕЙНЫЕ УРАВНЕНИЯ И
СИСТЕМЫ УРАВНЕНИЙ	,
...............о
Лабораторная работа № 1 .Распад разрыва в механике
сплошной среды.....................
Список литературы.....................27
Приложение............................28
РАЗДЕЛ II АНАЛИТИЧЕСКОЕ
ПРИБЛИЖЕНИЕ ФУНКЦИЙ............................39
Лабораторная работа № 2. Интерполяция. Сплайны.39
Список литературы..............................45
Приложение...................................  46
РАЗДЕЛ III ОБЫКНОВЕННЫЕ
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ	49
Лабораторная работа № 3. Метод Рунге-Кутты
решения задачи Коши для обыкновенных
Дифференциальных уравнений первого порядка.....49
62
Список литературы.....................
63
Приложение............................
з
РАЗДЕЛ IV. КРАЕВЫЕ ЗАДАЧИ ДЛЯ
ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ
- • оигний ВТОРОГО ПОРЯДКА.
66
Лабораторная работа № 4. Метод прогонки.....66
Список литературы...........................95
Приложение..................................96
РАЗДЕЛ V УРАВНЕНИЯ С ЧАСТНЫМИ
ПРОИЗВОДНЫМИ ГИПЕРБОЛИЧЕСКОГО
ТИПА.........................................Ш
Лабораторная работа № 5. Методы решения
уравнения переноса...................................111
Список литературы....................................132
Приложение...........................................133
РАЗДЕЛ VI. УРАВНЕИЯ С ЧАСТНЫМИ
ПРОИЗВОДНЫМИ ПАРАБОЛИЧЕСКОГО ТИПА. 162
"«ого уравнения теплопроводности......|62
С™с«™гера„ры......
......*....................1 оЭ
Приложение
...........................................186
4
Введение
Стремительное развитие и усовершенствование вычислитель-
ной техники в последнее время, когда её быстродействие порядка X
скольких терафлоп и память величиной в терабайты становятся доступ-
ными для массового пользователя, позволяет по-новому взглянуть на
методы решения математических задач, постановка которых уже из-
вестна и теоретически строго обоснованна, но решение всё ещё не по-
строено. Многие из этих проблем представляют интерес не только с
чисто математической точки зрения, но имеют важное практическое
значение, так как результаты их исследования используются при изу-
чении разнообразных физических, химических, биологических и других
естественно-научных явлений и процессов.
С другой стороны, простое знакомство с методами прикладной
математики, не подкреплённое самостоятельным опытом их реализации
на ЭВМ и получением конкретных решений, не гарантирует в будущем
их эффективного и правильного применения при решении конкретных
практических задач. Возникает необходимость выработки у обучаю-
щихся вычислительной математике целой системы навыков, включаю-
щей в себя: знакомство и умение пользоваться вычислительной техни-
кой и математическим обеспечением к ней, выбор наиболее подходяще-
го численного метода, его алгоритмическую реализацию в виде
программы на языке высокого уровня для ЭВМ, её отладку, проведение
тестовых и серийных расчётов, обработку и анализ полученных резуль-
татов.
Достижению этих целей как раз и призван способствовать вы-
числительный практикум по прикладной математике, который на акту-
альных примерах, имеющих физическое содержание, способствует
формированию и закреплению у обучающихся необходимых практиче-
ских навыков решения современных сложных задач вычислительной
математики и получения конкретных результатов, тем самым подтвер-
ждая обоснованность теоретических положений прикладной математи
ки.
5
раздел I
НЕЛИНЕЙНЫЕ УРАВНЕНИЯ И СИСТЕМЫ
УРАВНЕНИЙ
Лабораторная работа № 1
Распад разрыва в механике сплошной среды
Введение
Многие современные задачи механики сплошных сред связаны с
изучением возникновения и распространения различных разрывных
течений. На разрывах должны выполняться три основных закона сохра-
нения: массы, количества движения и энергии. К ним, как правило, до-
бавляют дополнительные соотношения, которые следуют либо из экс-
периментальных исследований, либо являются обобщением большого
числа натурных наблюдений. В ряде случаев оказывается, что решения
справа и слева от поверхности разрыва однородны, т.е. значения газо-
динамических величин не зависят ни от временной переменной, ни от
пространственных координат. С течением времени меняются только
размеры этих областей, но не значения физических переменных в них.
Это свойство оказывается справедливым в течение некоторого времен-
ного интервала.
мятииГг?™* условиях появляется возможность сформулировать мате-
нений п₽.Ю 3аДаЧу’ [1Редставляющую собой систему нелинейных урав-
оазоынпн КОТОрую можно определить скорость распространения
РассмотпимДОВаТеЛЬН0’ И размеРы соответствующих областей.
Рассмотрим несколько таких случаев.
РазРыва с Уравнением состояния идеального
* <13 «1
Одним из важных
Равного решения являете °ВИИ правильн°го описания поведения раз-
па поверхности разрыва п Вь1ПоЛнеНие основных законов сохранения
щих преобразований. 0ЛУчим полезные соотношения для дальней-
Для простоты останов
^*Те и будем рассматриваться НЭ ОДНомеРном пространственном вари-
°РДинат, которую На3овемаалачУ в неподвижной декартовой системе
Пусть справа от разрыва давление в среде Ро, массовая скорость
t7o> слева — давление Р^, массовая скорость . Предположим, что
сам разрыв перемешается по пространству в лабораторной системе ко-
ординат со скоростью 2)0. В этом случае законы сохранения массы и
импульса на разрыве могут быть записаны в следующем виде:
Р\ ОЛ ~ Л) ) = Ро (Ц) — ),	(1)
Рх +	\ — Do ) = Ро + p0(t/0 — Do) ,	(2)
где р0 и Р] - это значения плотности справа и слева от поверхности
разрыва.
Выражая Ux — Do из уравнения (1), подставляя во второе уравнение и
разрешая относительно (Uo — Do )2 , получаем
(и - d )! =	&
1 “ °’ р^-рЛ
Совершенно аналогично находим, что
Pi(Pi~Z?o)
Вычтем из (4) почленно (3) и преобразуем получившееся выражение:
(Ц2-2ВД-^о+2ад) =
HoHi
Выразим U} - D6 через Uo - Do из (1) и подставим в (5):
(А-А)(А-А)-А+^=-^йм-
Р\
Р\ + Ро
Сокращая на-—, имеем
Р\
Ро
Р\Ро
(3)
(4)
(5)
(6)
7
(7)
Извлечем	Ю,°" "	*
у, -D0«t6>иразр^имотносительно и -и,-
Ппи выводе соотношения (7) были использованы законы сохранения
массы и импульса. Добавим к ним закон сохранения энергии:
.	, Г (г0-о0)!
= (n,-O.V.
(8)
(Ю)
предварительно преобразовав его, где через £0, обозначены удель-
ные внутренние энергии справа и слева от точки разрыва. Так как по
предположению газ идеальный, то удельная внутренняя энергия, плот-
ность, давление и показатель адиабаты связаны соотношением
/ = о,1.	(9)
Вынесем из фигурных скобок в уравнении (8) справа - р0, а слева - Д.
Затем почленно поделим на (1):
+л+({/«-л0)2
Р1 2	0 7^-------------
^eM £0 И £] из (9) и подставим в (10):
fro ~1)р,	V- + --°-	(’
V'o-Opo 2
Выразим отношение пл
отношение давлений ИпТ”00*6” П0 Р^111116 стороны от разрыва через
р р р Г°"адста»“"(3)и(4)в(11):Р
. А(д-л)
2

После умножения нраяой и левой част„
произведение плотностей, приведен™ подобием JT"°™"
выражению	членов приходим к
А _ U + 1>о	(Гг>+1,+	‘') J
А ^-1>о+(го+1)Г(^)^7])#-	<,2)
^0
Найдем теперь отношение (Р, - Р„)/(и, -U„\ тюш Р, - Р, ш
Скорость звука в газе связана с плотностью и давлением выражением
I Р~
Со =	_0- Подставим Со и отношение плотностей из (12) в пра-
V Ро
вую часть (13). В результате имеем 
Допустим, что справа от разрыва значения плотности Pq-, скорости
Давления Ро, а также показатель адиабаты /0 нам известны, а
9
ачения слева от разрыва, имеющие индекс 1 и ско-
нужно определить зн	Соотношения (1) - (14) позволяют
роль Движения фро^^^ известно одно из искомых неизвест-
да, ««м "	„аия.го комбинация « них.
^™пеосновных варианта распада разрыва (1]:
образованием двух ударных волн;
XS».’«у4»»я юл"011" олвоя [’a,p,“c,™:
в-третьих, с двумя волнами разряжения;
в-четвертых, с двумя центрированными волнами разряжения
Остановимся более подробно на первых двух.
1.1.1. Распад разрыва с возникновением двух ударных волн
Пусть в начальный момент все пространство можно разделить на
два однородных полупространства, граница между которыми проходит
по плоскости разрыва. Для определенности будем считать, что газоди-
намические величины справа от этой плоскости имеют индекс 0, а слева
- индекс 3 (см. рис. 1).
Гз’Л’^З’Рэ’^З 1
I
Рис. 1
Со врсм	*
пространствам н^^Г раСПада Р^рыва по правому и левому полу- ‘
и™ индекс 1 газодинамическим1РаНЯТЬСЯ УДарНЫе волны- ПрИ’ '
нои и индекс 2 - за левой (с КИ1^ параметрам за правой ударной вол-
Контактный
10
Запишем три основных закона сохранения массы, количества дви-
жения и энергии на левой и правой ударных волнах и условия непре-
рывности давления и массовой скорости на контактном разрыве:
Р\ №1 —	) = Ро(^о ~ Д) )>

(t/j Do К /?! £t
Рз^Цз ^з)_Рг(^2 А)»
(15)

Рг~Р„иг=и,.
Используя формулу (14), из первых трех уравнений системы (15) полу-
чаем, что
а из четвертого, пятого и шестого уравнений —
(16)
11
Гр
//з  ---скорость звука е
( Рз
Заменяя в (16) Р2 на
уравнений с двумя неизвестными rj
I—
где Су
3 -Зк
в области «3».
Р} и U2 на (7,, приходим к системе из двух
______Jo
2 Го
(17)
3
Р,-Л
и,-и у
Р>
V-P' Г-Р’- а
Введем обозначения i — л ~ r> '	(	—iV
Г Л» vo У
а3 =	• Умножим первое уравнение (17) на Ux — UQ и поде-
v з “ У
лим на правую часть, а затем аналогично второе уравнение (17) умно-
жим на [/] - (/3 и разделим на правую часть:
П -TI ~ П— V_________У - I
1 ’ + уМл	’ «,г I
,, „ - С,т/2
1	-1) ,Л(Л <«,/) )
Возведем в квадрат первое и второе уравнения О В) н вычтем из от Р
первое:
-^з-ОГ-Уз+(7,-1/0) =
= Iy________________2с] (к -	... _
Гз(/з-1№ + а3У) /«(ro-OS+ao^
(18)
(19)
12
Подставляя в (19) вместо U}-U3 и U, -Uo их выражения из(18) и
сокращая на общий множитель, получим
Г 2cf (Г.Х1
1 Г, О') - 1№з - и,)! л|л-.,<чг)
± fzz^rzz._tiL=+1, . т
71 + <V
Обозначим
2С2
^’гзСгз-О^з-У.)2'
е 2С<
и сделаем соответствующую замену в (20):
±У\(1'-1)..±1	(21)
Jx(x+a,Y) +
Из (21) нужно определить неизвестное У. Возведем в квадрат левую и
правую части (21):
е,(У-Х)!	е0(У-1)~ , 2^(Г-Х)(Г-1) т
X(X + a,Y) (1 + а0Г)	,/.«-(% +а,Г)(1+а0У)'
Затем еще раз возведем в квадрат, но уже равенство (22):
13
г/(Т-Т)'(1+«0Т)-’ +x!(x+ajr)2(i +a„yy _
- 2 Vj Х(Х+а, Г)(1 + »0У)(Г - X)2 (Г-1)2 +
+e2X2(X + a3Y)2(Y-l)* -
^XiX + aXXl + a.r^Y-X)2-
~2е,Х2(Х+а!Г)2(1+Й0Г)(У-1)2 =0.
Если в (23) привести подобные члены и сгруппировать их
хл^г.»upnowin ивдимпои. члены и группировать их при одинако-
вых степенях У, то получим алгебраическое уравнение шестой степени:
aJ6+a{Y5+a2Y*+a3Y3+a4Y2+asY + a6=0,	(24)
где коэффициенты a:(i = 0-5-б) выражаются через ранее введенные
константы:
Оо ~ (®о^3 ®3^0) ’
О] = 2| (ttee3 — л3Ае0)[е3(1 —
~айа3Х(аое3+а3Хей) };
а2=е2(ба2/2 -8аД + 1)-
-2е0е3^[а0а3(х2 +4Х +
+ е2Х2(ба2 -Ъа.Х + Х2У
з
аз + «0
ДЛШ Jikt -леи и*» -

X
— ла3л +04 J~o(.o<x3%(аoA’ + (Xз) +
^ ^з £ cZqOCj А" 2Х ( 2аоа 3 + cig^T + (2cXqJT + ос.
+ еоА'[аоаз-2a3(a3+2anA') + 2or Ул-у- v2
а4 = X2 { е2 (а02Х2 - 8а0Х + б)-
- 2е0е3 [аоа3А" ~ 2(Х +	+ ао%)+ X2 +4Х +1]+
+ во(а3 -%а3Х + 6Х2)+(а1 + 4айа3Х + а%Х2)-
- 2е3 (&оХ + 2аоа3 )а" - 2(2а0Х + а3)+1]~
-2е0 а3(2а0Х + а3)-2Х(2а3 +а0А')+Х2] },
as =2Jf3[ е3(а0Х ~2)-е0е3(а0Х-2 + а3-2Х)+
+ ео(а3 -2Х)+(а3 +аох)~
-е3(2а0Х + а3 -2)-е0(2а3 +0^-2^) ],
аб - ЛГ4[(е3 -е0)2 +1 -2(е3 + е0)].
Общие подходы к нахождению корней алгебраических уравнений будут
рассмотрены ниже, а сейчас рассмотрим вариант распада разрыва с об-
разованием ударной волны и волны разряжения.
1*1.2. Распад разрыва с образованием ударной волны
и волны разряжения
Снова, как и в пункте 1.1, будем предполагать, что в начальный мо-
мент все пространство можно разделить на два однородных полупро-
^анства, газодинамическим величинам в которых припишем индекс О
правого и индекс 3 для левого полупространств. Граница раздела
Между ними - плоскость (см. рис. 1 ирис. 2).
15
мтр паспада разрыва по правому полупространству начина-
В реЗУ'Тняться ударная волна, а по левому - волна разряжения. На
ет распространят , ы выполняться три основных закона сохранения
^Г^Хдв-нияиэНергаи:
Pi(U -D0)=p0(U0-D0),
^+Р,(^-/)о)2=/>о+а(^о-А)2.	-
а затем, используя соотношения (25) и (26), получить для волны разре-
жения:
Тогда с учетом равенств (27) и (28) система из семи уравнений может
быть сведена к системе из двух уравнений с двумя неизвестными Р} и
В волне разряжения на характеристике семейства U + С сохраняются
значения инварианта (25) и справедливы соотношения для адиабатиче-
ского процесса (26):
и,+^-=и2+^, 
Zj-1 /,-1
(25)
л-1
' 2
3 А
‘ (26)
ные coZT* гРанице остаются непрерывными давление и нормаль-
ные составляющие скорости:	и Е «
ri=r2>	(27)
Вкон U'=U1'	<28)
мью иеизвестнымГл110;?11^ К°Р"И СИСТемы 113 семи УРавнеНИЙ С
и*" Удобно записать f*’ 1 ’ f1 ’	’ С2 ’ ^2, А • Первые три уравне-
“писать в форме (14):
Система уравнений (29) допускает дальнейшее упрощение, если ис-
ключить из нее Uy:
2С
: иУ-и0=^~
’	•=; Лз-1
17

ГОУ 1)110
“Московский
фнзико* технический институт
(юсу дарственный увнвнреиъ’ <) мг jg
БИБЛИОТЕКА* ”
16
Перенесем из правой части (30) первое слагаемое влево
квадрат и умножим на знаменатель:	’ веДем в
р р	/
Z3 zo	J
2s0
п (31)
2С,
^°-уз_1
ъ-'
-1
л
где £й = Л/ко (/о -0] - начальная удельная внутренняя энергия в
правом полупространстве.
Если ввести обозначения CtQ

1.2. Распад разрыва с уравнением состояния типа
Грюнайзена (Забабахина)
При высокоскоростном ударе конденсированных веществ система
уравнений (15) должна измениться, так как функциональная зависи-
мость между давлением, плотностью и удельной внутренней энергией
оказывается другой, отличающейся от уравнения состояния идеального
газа. Е.И. Забабахиным было предложено уравнение, устанавливающее
связь между давлением и плотностью, как за фронтом ударной волны в
конденсированном веществе, так и при изэнтропическом сжатии или
расширении:

'о
2
у =---
/3-Ц
(33)
2 Ро Рз
1 з ,
Р0’
тогда уравнение (31) можно записать в виде
X2Z2n-2XZn + 1 = |д2-2^v(Z-1)+v2(Z-1)2 Kl + aVZ")
или, раскрывая скобки и приводя подобные члены,

X2Z2n -a0v2AZ"+2 +2a^(fi + v)XZn+x -
-[2 + (^ + v)2a0]xzn-v2Z2+2v(^ + v)Z-
-(д + у)2 +1 = 0.
Г0 Р^Рива удалось систем^* вариантов распада газодинамиче-
мого порядка свести к г>п Ы Нелинейных уравнений восьмого и
паицТеЛЬНЫМИ КоэФФиЦиентами°ьгУ алгебРаическому уравнению с дей-
вов „еСКИх уРавнений вернемся ВопР°сам нахождения корней алгеб-
ХХ““н"“	"гнм с,уч”‘ ₽“пма ри₽и;
,м от Уравнения состояния идеала
где /0- показатель адиабаты конденсированного вещества, /Эо- плот-
ность, Сд - постоянная вещества (условно скорость звука),
Ло = (/о + 0/(/о	плотность за фронтом ударной волны,
р , Р - плотность, давление в волне сжатия или разрежения, если она
формируется в области за ударным фронтом. Если область за фронтом
ударной волны однородна, тогда выражение (33) упрощается и превра-
щается в адиабату ударного сжатия (Гюгонио) для конденсированного
вещества:
А
Р = РйС2д
kg ~
А
(34)
1.2.1. Распад разрыва в конденсированном веществе
с образованием двух ударных волн
Допустим, что в начальный момент происходит соударение двух
однородных полупространств, состоящих из конденсированных ве-
19
is
wn.,„ занимает область положительных значений
Шеег«, >« ” ““"Ответствует индекс «О», и второго - заиимающе-
XXX »»^-» ™еюш'г“индек <,3>>(от'рис- ’
/3;P3;t/3;A I A’^oi^o’/o
Рис. 3
После соударения при достаточно большой относительной скорости по
правому полупространству начинает распространяться ударная волна
£>0 с давлением , плотностью р} и массовой скоростью Uх за ее
фронтом. Аналогично по левому полупространству перемещается удар-
ная волна со скоростью £)3, имеющая давление Р2, плотность р2 и
массовую скорость С/2 за фронтом (см. рис. 4).
Z)3	Do
Хз»Л»^3»А А’^2»Р2’Хз '	А А’^О’^О’/о
контактный
разрыв
Рис. 4
-iff №Т'Я№ '	*
Разделяет эти два полупространства контактный разрыв. Будем считать,
что на каждой из ударных волн выполняются законы сохранения массы,
количества движения, справедливы адиабаты ударного сжатия (Гюго-
НИО), а на контактном разрыве - непрерывны давления и массовые ско-
рости. Система из восьми уравнений с восемью неизвестными D3, р2,
^2 > Р2 > Pi > , р\, Dq принимает следующий вид:
p,(t;i-Do)=po(t/0-JDo),
Р* +	“Do)2 = Ро + p0(U0 - Do)2,
/> = а(а г а)
А^о — Р\
20
Рз + РзФз Д})2 - А + Р2(^2 ~ ^з)2.
р -_а^Р2~Рз)
2 РзК ~ Pl
/2 = /|, t/2 = I7j ,
где а0 = р0С0 (h0 -1),а3 = р3С\ha=(y0+ 1)/(уа -1),
h3 =(/з+1)/(/з-1)-
Первые два уравнения системы (35) согласно формуле (7) можно пере-
писать в виде
U -U =+	(36
UlUo +3	’	(jOj
V Р\Ръ
а четвертое и пятое - в виде
и2 - и3 = + (ЕНВЖНЫ.	о?)
\! РзРз
Вычитая почленно уравнение (36) из (37), учитывая соотношения семь
и восемь из системы уравнений (35), получаем
и ,, _	?(^-^)(а-а) ,эд
Мл (J Т--к /--------------Г /
v лро i
Третье и шестое уравнения из системы (35) разрешим относительно рх
и Р2 и заменим в последнем Р2 на 
pt =	(39)
Ро Р\+ао
PL_Plh3 +а3	(40)
Рз Р+а>
Преобразуем подкоренные выражения в (38), поделив числитель и зна-
менатель первого слагаемого на р0, а второго на р3
21
Подставим (39) и (40) в (41) и введем обозначение
И=(/о-£/,=±
Р<Ж + ао)
(42)
р3(РД + а3)
Возведем левую и правую части (42) в квадрат, затем перенесем в левую
часть слагаемые, не содержащие радикалов:
_	1 Wr Ъ) _	(д - л)=
+	р,(Р,А,+а,)
= ±2/>	<4
Ч АрД/’Д+“о)(ЛАз+«3)
„ v Р\ v	/? — а°___
Введем обозначения 1 --, Л -------—у, р$ -	2 ’
РьРУ РьРУ РьРУ
В =
Ру РоР^2’
в. • (А» - 1)а . <?, = (А, - 1)а • « = р^уг
и возведем в квадрат левые и правые части (43). Получившееся уравн
ние можно записать в виде
£а,У‘-'=0.	<*”
/-0
где До = (<% ^з^о ) >
Д| ~ 1[ h380 (/?3 — й3е)+ h083 fa — A7i0)—h0h3(80h3 + 83h0)—
-8й83(Ий03 + h30o)+8o81hoh3(X +e) ],
(45)
<h -ДО + 52 (h]e2 + P3 - 4рДе) + З2 (p2 + %Хг - 40A) -
—230Л)(2р,/^ + Ajp, —	25э/^,(2рой,+/^Pj —
-26,3, (Xeh^h, + P,p,) + 23,3, (X + е)(РД + 3,p,),
o,=4
A(M’-«₽»A>2)-2W.(P»-«A,)-S,(3,PS-xpa1)-
-23,РЛ(Р, -^)-e3„8,%(Ap, +А,р,) + 3,8,рд(лг + е) ],
Я, = Д!А? + АЧ! + 4ДД,Л,Й, + 51РУ +6^Х2 -
- 2<?, {р,р2 - epjh, - 2eh,p,p,)-
-2S,PMP, -Xh,)-2p,Xh,]-2e6tS,pap,X,
я,	+ Pi^o)+ ДоР^^РйеРг + Xp,i)],
= PlPi 
1.2.2. Распад разрыва в конденсированном веществе
с образованием двух ударных волн н нулевым
начальным давлением
В случае столкновения невозмущенных конденсированных веществ,
когда плотности в них соответствуют исходным значениям, начальное
давление равно нулю. Это упрощает выражение (43), т. к. в нем можно
положить Ро = 0 и Р3 = 0:
23
(46)
2 1Л I... _ Ч/	\ “
¥ -
Чр.л(/,А+а«^1,>+а
Приведем (46) к общему знаменателю и сократим на него, а затем поде-
™ вег уравнение и (д/2)  В результате получим
р(/ J
Р, , а3 '
р> к хД»_ -5-
рТрЛ=,ч> рои*Лрогг
-(V!)P3
•2
(47)
"(Vl)Po
t/2
= +2
x
) . аз
2 n t/2

X
РЛ ! ao
Ipo^2 РоИ2
Введем обозначения Y = !—	PL = _J^o—	о _ ——
Р<Г’	P° А/2’	Рз p0V2
е=Рз_
p и возведем левые и правые части (47) в квадрат. Тогда уравне-
ние относительно У примет вид
=0,
/=0
(48)
гдеао=КЛо-1Хе-(й3-1>0]2,
°' =Л	-ОЧд -(Ч-iyh^ -
-( ,-fo -!)(*, -1)е(^д + V,) ],

24
«2 =«%2 +(л0 -1)V# +(й3 -i)2 д2 -
-2(Л0-1)е2А3(ДЛ3+2А0Д)-
-2(й3 -1)ей0(й0Д + 2й3/?0)-2(й0 -1)(й3 -1>ДД,
«з =2[ е2кок3(0ок3 +й0Д)-(й0 -1)е2Д(/?Л+2й3/?0)-
-(л3-1)еД(Мо+2Лод)],
«4 =е2(р2к2 +к2^2 +4к0км)-2(к0 -1>2ЛА2 -
-2(Л3-1>Д/?02,
«5 = 2е2ДЛ(Л0Л +Мо), «6 =е2Д2Д2.
1.3. Численные методы решения алгебраических
уравнений
В предыдущих пунктах 1.1 и 1.2 с помощью математических преоб-
разований удалось свести системы нелинейных уравнений (15), (25) -
(28), (35) к отдельным алгебраическим уравнениям (24), (32), (44), (48),
из которых нужно определить искомые неизвестные Y или Z. Задачу
нахождения корней алгебраического уравнения можно разделить на две
подзадачи: во-первых, локализации корней, а во-вторых, уточнения
значения каждого корня до нужной степени точности.
Под задачей локализации корней подразумевается указание таких
частей области определения уравнения, в которых существует только
по одному его решению.
Допустим, что имеется алгебраическое уравнение вида
f(z) = aozn +a{zn-' +...+a„_lz + a„ = £a,zn~' =0	(49)
i=0
с действительными коэффициентами , i = 0 -г п (п - натуральное
число больше единицы), решение которого z может быть как действи-
тельным, так и комплексным числом.
Введем в рассмотрение два числа:
25
|о | и 5 = max{|a0'’'ai:’’ 1а"-' где
А . maxlie,,	- •! д	следствием из основной тео-
— СЛеДУЮШ“ УП~
И: Все кории алгебраического ура— (49) Р— -»
комплексной плоскости в кольце
__Ы_ <lz|<l+r-
(50)
Существует еще ряд признаков, позволяющих определить область на-
хождения корней. Приведем два из них.
ными значениями функции /(г) на границах меньше п, то процесс
можно повторить, например, удвоив число частей разбиения
Д/ = 22V. Продолжим этот процесс, пока не будут локализованы все
«-корней уравнения (49). Реализация этого алгоритма на языке про-
граммирования высокого уровня, как правило, не вызывает никаких
дополнительных сложностей.
После того как решена подзадача локализации корней, уточнение
приближенного значения корня до нужной степени точности может
быть выполнено методом половинного деления, методом простой ите-
рации, методом Ньютона или каким-нибудь приближенным методом с
обязательным контролем точности приближенного значения корня и
выполнения всех условий сходимости выбранного метода [3].
Метод Ньютона
Пусть многочлен из (49) имеет действительные коэффициенты
Я (j = 0 ч- и) и Яо)О. Тогда если при z = с выполнены неравенства
f(c)>Q, f'(c)>Q,..., f<n)(c)>0,
то число с служит верхней границей положительных корней уравне-
ния (49).
Список литературы
A',*•****"
Теорема Декарта
Число положительных корней уравнения (49) с учетом их кратности
равно числу перемен знаков в последовательности коэффициентов
°о> Op..., ол_р ап, (причем равные нулю коэффициенты не учиты-
ваются) или меньше этого числа на четное число.
С другими аналитическими подходами к локализации корней можно
ознакомиться в работах [2], [3].
После того как с помощью неравенств (50) или еще каким-либо дру-
гим способом определена область существования корней, подзадачу
л°“ции можно свести к выполнению следующего алгоритма
Вея область деящся на » частей, гас N _
««Г’”™	Ж Для дальнейшего
рассмотрения оставляют те подобласти и«
XW принимает значения разных ?	“™Р“
стей равняется «-степени яп к	“ ЧИСЛ° ТаКИХ подобла-
ча локализации считается решенной^^д0 МН°Г°ЧЛена (49>’ то подзада-
решеннои. Если же число подобластей с раз-
1.	Зельдович Я.Б., Райзер Ю.П. Физика ударных волн и высоко-	|
температурных гидродинамических явлений. М. Гос. Изд.	В
Физ.-мат. лит., 1963.	I
2.	Курош А.К. Курс высшей алгебры. — М.: Наука, 1975.
3.	Демченко В.В. Распад произвольного гидродинамического
разрыва: Учебное пособие. — М.: МФТИ, 1998.
26
27
ПРИЛОЖЕНИЕ
Вариант № 1
Распад плоского разрыва в одномерной газовой динамике может
происходить с образованием двух ударных волн и контактного разрыва.
Припишем газодинамическим параметрам перед первой ударной вол-
ной индекс «0»; между фронтом первой ударной волны и контактным
разрывом - «1», между контактным разрывом и фронтом второй удар-
ной волны - «2», невозмущенному газу перед второй ударной волной -
«3». Тогда справедливы соотношения Рэнкина-Гюгонио на фронте
ударных волн в лабораторной системе координат:
л(Ц _ Do)= P0(U0 - Dq\
р,+л(Ц-по)! = ро + а(ц -dJ,
2
_ (^о Во л А А) +
2
r! + P,(C,-O!)1=n + P!(U!-D,)!,
(и,-Р,)кр,Л’~Р»У]+р,.=
и соотношения на контактном разрыве
р - п
и2=и15
где р - плотность, U- массовая скорость газа, Р — давление,
С - скорость звука в газе, D - скорость фронта ударной волны, у - по-
казатель адиабаты в соответствующих областях. Определить скорости
ударных волн Dq и D3.
Задание № 1
Го = 5/3; р0 = 2,219-IO"4 г/см3; Uo =-1,587-105 см/с;
Ро = 3,7812 • 106 дин/см2;
/3 =5/3; /?3 =2,71-10~3 г/см3; U3 =10 см/с; .
Р3 =5-105 дин/см2.
Задание №2
у0 =5/3; р0 =7,9 г/см3; Uo =0; Ро = 3,04-109 дин/см2;
/3 = 5/3; р3 = 11,37 г/см3; U3 = 5 • 104 см/с;
Р3 = 1,17928-109 дин/см2.
Задание № 3
/о =5/3; р0 = 2,71 • 10'3 г/см3; Ро = 5• 105 дин/см2; Uo = 0;
/з =7/5; р3 = 2,219-Ю"4г/см3; U3 = 1,587-Ю5 см/с;
Р3 = 10б дин/см2.
Задание №4
Го =5/3; р0 =7,9 г/см3; U. =0; Ро = 3,04 • 109 дин/см ;
/3 = 7/5; рз = 11,37 г/см3; U3 = 5 • Ю4 см/с;
Р3 =1,17928-109 дин/см2.
Задание №5
Го =7/5; р0 = 7,9 г/см3; Ро = 3,04-109 дин/см ; Uo = 0;
у3 = 7/5; р3 = 11,37 г/см3; U3 = 5 • Ю4 см/с;
Р3 =1,17928-109 дин/см2.
28
29

Задание М 6
у0 =5/3; р0 = 11,37 г/см3; 1/0 = -2,28 • 104 см/с;
Рй =1,17928-10’ дин/см2;
у, =5/3; р3 = 7,9 г/см3; U} = 2,72-104 см/с;
Ру = 3,04 • 10’ дин/см2.
Задание № 7
Го =5/3;	=7,9 г/см3; Ро = 3,04-\09 дин/см2;
Uo=-2,72-104 см/с,-
/з = 7/5; Ру = 11,37 г/см3; U3 = 2,28 • 104 см/с;
?3=U7928-l09 дин/см2.
Задание Ns 8
у0 =7/5; pQ = 11,37 г/см3; Uo = -2,28-104 см/с;
Po= 1,17928-10’дин/см2;
Гу = 7/5; Ру = 7,9 г/см3; Uy = 2,72 • 104 см/с;
Ру =3,04-10’ дин/см2.
Вариант № 2
Распад плоского разрыва в одномерной газовой динамике может
происходить с образованием ударной волны, контактного разрыва и
волны разряжения. Припишем газодинамическим параметрам перед
ударной волной индекс «0», между фронтом ударной волны и контакт-
ным разрывом - «I», между контактным разрывом и «хвостом» волны
разряжения - «2» и невозмущенному газу перед волной разряжения - .
«3». Тогда справедливы соотношения Рэнкина-Гюгонио на фронте |
ударной волны в лабораторной системе координат:	:
P,(U,-D0)=p0(U0-oj	:
А+рЛц-о.Ьп+Мп.-р.у,	;
30
(t/.-A) A
2
— (Ц) A>)- Po
соотношения на контактном разрыве
и в волне разряжения
и3+-^-=и.
n-i
2С.
+----
где р - плотность, U - массовая скорость газа, Р - давление,
С - скорость звука в газе, Do - скорость фронта ударной волны, /
показатель адиабаты в соответствующих областях. Определить ско-
рость ударной волны Do.
Задание № 1
Уь =5/3; р0 = 1,694 10^ г/см3; Ро = 1,013-106 дин/см2;
Ц>=0;
Г3 =7/5; С3 =3,6537-104 см/с; Р3 = 1,6768-106дин/см2;
U3 =0.
Задание № 2
уй =5/3; рй = 1,694-1О4 г/см3; Рй = 1,013-106 дин/см2;
/з =7/5; С3 =3,6537-104 см/с; Р3 =1,6768-106дин/см2;
Uo =0;U3 =1,229-104 см/с.
31
п=5/3;	U. =0^. = 3,848-10= дан/см =;
/з = 5/3; С3 =1,31478-Ю4 см/с;
щ = 5 • Ю4 см/с Р3 = 1,17928 • 109 дин/см2.
Задание №4
уа = 7/5; А = 1,694 • 10"4 г/см3; Ро = 1,013 • 106 дин/см ;
UQ = Ю'3 см/с,-
/з = 7/5; С3 = 3,6537 • 104 см/с; Р3 = 1,6768 • 106 дин/см2;
U3 = 1,229-104 см/с.
Задание № 5
/о = 7/5; р0 = 10'5 г/см3; Рй = 3,848 • 103 дин/см2
С/о =0;
/3 =5/3; С3 =1,31478-Ю4 см/с; U3 = 5-104см/с;
Р3 =1,17928-Ю9 дин/см2.
Задание Ns 6
/о = 5/3; р0 = 1,694-10^ г/см3; Ро = 1,013 -106 дин/см2 •
Uo = ~Ю 1 см/с,-
г. =S/3; С, =3,6537-10’ см/с; Р, =1,6768.10‘дан/смг;
Задание Ns 7
/о =5/3; рй =1()-5 r/CMs.
/'з =5/3; С3 = 2,53248-IO4
Л =3,04-109 дин/см2.
• -01Л, =3,848-10= дин/см:.
см/с ; из = 0;
32
Задание № 8
/о = 5/3; Ро =10'5 г/см3; Ро = 3,848 • 103 дин/см2;
Ц>=0;
Уз = 7/5; С3 =2,53248-104 см/с; Р3 =3,04-109дин/см2;
U3 = 0.
Вариант № 3
Распад плоского разрыва при высокоскоростном ударе конденсиро-
ванных веществ может происходить с образованием двух ударных волн
и контактного разрыва. Припишем газодинамическим параметрам пе-
ред первой ударной волной индекс «О», между фронтом первой ударной
волны и контактным разрывом - «1», между контактным разрывом и
фронтом второй ударной волны - «2», невозмущенному веществу перед
второй ударной волной - «3». На ударных волнах выполняются законы
сохранения массы, количества движения и справедливы адиабаты удар-
ного сжатия (Гюгонио):
Р\ (^/| — А ) = Рй (Ц, ~ Df) ),
Р,+р|(С/,-£>0)!=Р0+р„(ио
р _ «о(/>|-Ро)
РЛ~Ру
p3(U3 — D3) = рг (и 2 — D3),
Р3 + p3(U3 -D3)2 = Р2 + p2(U2-D3)2,
р ^аЛр2~Р^
Рз^з ~ Рг
а также соотношения на контактном разрыве
Рг=Р^
и2=и.,
где р - плотность, U- массовая скорость вещества, Р - давление, D
скорость фронта ударной волны; у, - показатель адиабаты,
а, = ptC*(ht = (/, +1)/(/; -1), С- постоянная i-
33

вешества, / - <U • еоотве^уюшнх областях. Опряишп. скорое™
ударных волн Dq и D3-
Задание № 1	,	.
Гй =3;	=13.76163г/см3; Uo = -5,322389-10 см/с;
Ро =5,176613 10'2 дин/см2; Со =4,65-105 см/с.
п = 3; р3 = Ю2 г/см3; U3 = 0 см/с; Р3 = 0 дин/см2; ,
С3 =4,65-105 см/с.
Задание № 2
у0 =3; р0 = 102г/см3; Uo = -106см/с; Ро = 0 дин/см2;
Со = 1,972-105 см/с.
Уз =3; р3 =21,80089 г/см3; U3 =-5,322389-105 см/с;
Р3 =5,176613-1012 дин/см2; С3 = 1,972-105 см/с.
Задание № 3
/о =3; р0 = 13,76163г/см3; Uo = 0 см/с;
Ро = 5,176613 • 1012 дин/см2; Со = 4,65 • 105 см/с.
Гз =3; р3 =102г/см3; U3 =5,322389-105 см/с;
Р3 = 0дин/см2;С3 = 4,65-105 см/с.
Задание № 4
Го =3; р0 = Ю2 г/см3; UQ = -4,677611 -Ю5 см/с •	‘
Ро = 0 дин/см2 ;со= 1,972-Ю5 см/с
Хр3=3;Рз=21,80089г/см3;С/3=0см/с-
Рз = 5,176613-ю12 дин/см2- г ,
дин/см , сз =1,972.105 см/с
Задание № 5
Zo = 3, р0 - 10 г/см ; Uo = 0 см/с; Ро = 0 дин/см2;
Со =4,65-105 см/с.
Уз = 3; Рз =13,76163 г/см3; U3 = 5,322389-105 см/с;
Л = 5,176613-1012 дин/см2; С3 = 4,65-105 см/с.
Задание № 6
Уо =3; р0 =21,80089 г/см3; Uo = 5,322389-105 см/с;
Ро = 5,176613-1012 дин/см2; Со = 1,972-105 см/с.
Уз =3; р3 = 102 г/см3; U3 = 106 см/с; Р3 = 0 дин/см2;
С3 = 1,972-105 см/с.
Задание № 7
Уй =3; р0 =102г/см3; Uo = -5,322389-105 см/с;
Ро = 0 дин/см2; Со = 4,65 • 105 см/с.
у3 =3; Рз = 13,76163 г/см3; U3 =0 см/с;
Р3 = 5,176613-1012 дин/см2; С3 =4,65-105 см/с.
Задание № 8
/о =3; Ро =21,80089 г/см3; 1/0 = 0 см/с;
Ро = 5,176613-Ю12 дин/см2; Со =1,972-Ю5 см/с.
/3 = 3; р3 = 102 г/см3; из = 4,677611 • 105 см/с ;
Р3 = 0 дин/см2; С3 = 1,972• 105 см/с.
Вариант № 4
Распад плоского разрыва при высокоскоростном ударе конденсиро-
ванных веществ может происходить с образованием двух ударных волн
и контактного разрыва. Припишем газодинамическим параметрам пе-
ред первой ударной волной индекс «0», между фронтом первой ударной
35
волны и контактным разрывом - «1», между контактным разрывом и
фронтом второй ударной волны-«2», невозмущенному веществу перед
второй ударной волной — «3». На ударных волнах выполняются законы
сохранения массы, количества движения и справедливы адиабаты удар-
ного сжатия (Гюгоиио):
Pi0-А ~ Р^о)~ Ро(^о ~)>
^+а^,-Д)2=Л+а(^-Д)2.
Ро^о Pi
P3(t/3	~^з)’
Рз^Рз^з'^ = + Р^2 ~D^’
р
Рз^з ~ Рг
а также соотношения на контактном разрыве	'
А=Л-
(/2=(/„
где р- плотность, I/ - массовая скорость вещества, Р - давление, D - :
скорость фронта ударной волны; уг показатель адиабаты,
ai ~ РР, (fy _1), =(/,	“0> - постоянная ь вещест-
ва, У - 0,3 в соответствующих, областях. Определить скорости ударных
волн Do и D}.
Задание № 1
=3; р0 =11,346 г/см3; Uo -106 см/с ; Ро = О дин/см2 ;
Со = 1,972-105см/с.
Уз =3; Рз = 7,85 г/см3; 17э = О см/с; Р3 = О дин/см2;
С3=4,65-105 см/с.
Задание № 2
Zo = 3; А = * 1346 г/см3; Uo = -4,677611 • 105 см/с;
Ро = 0 дин/см2; Со = 1,972 • Ю5 см/с.
Z3 =3; р3 = 7,85 г/см3; U3 = 5,322389-105 см/с;
Ру = 0 дин/см2; Су = 4,65 • 105 см/с.
Задание № 3
Zo =3; р0 = 7,85 г/см3; Uo = 0 см/с; Ро = 0 дин/см2;
Со =4,65- 1О5 см/с.
Уу =3; Ру =11,346 г/см3; U3 =106 см/с; Р3 =0 дин/см2;
Су =1,972-105 см/с.
Задание № 4
Zo =3; р0 =7,85 г/см3; Uo =-5,322389-105см/с;
Ро = 0 дин/см2; Со = 4,65 • 105 см/с.
Г3 = 3; р3 = 11,346 г/см3; U3 = 4,677611 • 105 см/с;
Р3 = 0 дин/см2; С3 = 1,972 • 105 см/с.
Задание № 5
Zo =3; р0 = 7,85 г/см3; Uo = -106 см/с; Ро = 0 дин/см2;
Со =4,65-105 см/с.
у3 =3; р3 = 11,346 г/см3; U3 = 0 см/с; Р3 = 0 дин/см2;
С3 =1,972-105 см/с.
37

г
Задание N° 6	2
=3. Ра = 11346 г/см’; U. . О см/с; Л = О дин/см ;
со = 1,972 1 05 см/с.
л = 3; А = 7.85 г/см’; U, = 10* см/с: Р, = 0 дин/см ;
Су =4,65-105 см/с.
Задание № 7
п = 3; р0 = 7,85 г/см3; U6 = 5,322389 • 105 см/с;
Ро = 0 дин/см2; Со = 4,65 • 105 см/с.
Гу = 3; р3 = 11,346 г/см3; U3 = -4,677611 • 105 см/с;
Ру = О дин/см2; С3 = 1,972 • 105 см/с.
Задание № 8
У0 =3; р0 = 11,346 г/см3; UQ = 4,677611 -105 см/с;
?о = 0 Дин/см2; Со = 1,972 • 103 см/с.
Л =3; а =7,85г/см3; U3 = -5,322389-Ю5см/с;
Дин/см2; С3 =4,65-103см/с.
38
РАЗДЕЛ II
АНАЛИТИЧЕСКОЕ ПРИБЛИЖЕНИЕ
ФУНКЦИЙ
Лабораторная работа № 2
Интерполяция. Сплайны
Введение
Часто при проведении научных исследований приходится
иметь дело с функциями, значения которых заданы в конечном числе
точек. Это могут быть результаты каких-нибудь экспериментальных
работ или физических наблюдений; таблицы, полученные при выпол-
нении расчетов на ЭВМ; основные точки будущей конструкции или
восстанавливаемого изображения и т. д. Необходимо как можно точ-
нее по имеющимся значениям представить поведение функции на ин-
тервалах между заданными точками, привлекая всю известную до-
полнительную информацию о свойствах функции, например, о её не-
прерывности и гладкости, возможной величине погрешности
заданных величин, её периодичности, четности или нечетности, мак-
симальных и минимальных значениях.
В данной работе будут рассматриваться способы аналитическо-
го приближения функций, о которых известно, что вся имеющаяся
информация достоверна с нужной степенью точности, а искомая
функция непрерывна и обладает требуемой степенью гладкости, т. е. у
неё существуют во всей рассматриваемой области производные необ-
ходимого порядка. Такой метод восстановления функции принято
называть интерполяцией.
ПЛ. Постановка задачи
Условимся называть совокупность точек, в которых заданы
значения функции, сеткой, а сами точки - узлами сетки. Если решает-
ся задача интерполяции, то узлы сетки называют также узлами интер-
поляции {хя }, п = 0 -г N.
Для простоты остановимся на случае функции от одного неза-
висимого переменного. Решение задачи интерполяции начинается с
39
„ейв0 Кез»»*»* ФУ"™Й	KOTOpM
“X XX»	—чго лю“ о6о6ше“-
л *а" "°-
п=0
лнн а ^0 не должен иметь на отрезке интерпо-
стоянные и хотя один ап * о, н
(г (х (х (Хы ) больше, чем А нулей. Считается, что среди
точек хп нет совпадающих. Существует несколько таких систем, на-
Ii v „2 Л I _ тогда говорят об алгебраической ин-
пример, р, X, X	> •••) м
терполяции, или {l,sinx,COSX,sin2x,COs2x, ...} - тригонометри-
ческая интерполяция.
Задача интерполяции ставится следующим образом: для задан-
ной системы узлов {xt) и соответствующих значений функции
= предложить обобщенный интерполяционный мно-
гочлен (совокупность постоянных ап), который в узлах интерполя-
ции принимал бы значения функции, т. е.
= o-w.	(1)
п=0
Из лекционного курса известно, если выполнены все ранее оговорен-
ные условия, то решение задачи интерполяции существует и единст-
вснно.
п«,ши В дальнейшем остановимся только на алгебраической интерпо-
•ИлЦИИ.	*
П.2. Алгебраическая интерполяция
терполяции сводится к оешеКЦИИ ЯВЛЯК>ГСЯ базисными, то задача ин-
порядка с определителем в!н ЛИНеЙНь1х Уравнений W + 1
коэффициентов ап,п = 0 “деРМ0НД\0™0СИтельн° неизвестных
определитель ВандерМ0нда 3 ЛИНеЙН°Й адгебРЫ известно, чт«
40
!	1	*0	*0	-	х"
I	1	X!	X,2	...	х/
1 х х2 х^
1	ЛН ••• ЛУ
отличен от нуля, если среди точек хк нет повторяющихся,
т. е. X, Ф Xj , если i & j, i,j = Q-^N. Так как определитель систе-
мы линейных уравнений отличен от нуля, то её решение всегда суще-
ствует и единственно.
Известны несколько удобных форм записи этого решения (алгеб-
раического интерполяционного многочлена).
Алгебраический интерполяционный многочлен в форме Ньютона за-
писывается в виде
W л-1
PN (*) = 2 Ь« П (Х " Х‘ ) = Ь<> + (Х " Х0 ) +	„.
л=0	/=0	.	(2)
+ .../>ЛГ(х-х0)...(х-хЛГ_1),
где коэффициенты Ьп определяются последовательно из
уравнений (1):
Положим сначала в (2), что х = х0, и приравняем согласно (1)
Ря (х0) = /о  Имеем b0 = f0. Далее аналогично X = X] и
Ря(х])=/; ... х = х, и P„(xJ = f для i< N. В результате на-
ходим все bn = /(х0,х,,...,х„), n = 0 + N, которые называют
разделенными разностями и-порядка. Полученный многочлен (х)
41
им интерполяционным, так как его
Де8с«№«»	Ь-
степень не выше V н по	закЛ10чение следует из единствен-
удовлетворит УРавнен'^М Э
““ГО,ЛЯ* ’
...+ 
Ш=1УЛ(Х _х.) '°(х0-х1)...(х0-^)
ияО <»0 к «	''	4
/*я
(х - х0)...(х - V1 )(х - XiJj* ~	+ - + (3)
+^' u,-x0)-"(x»"x"-'Xx" -x"+i)-^n 'Хы'
(x-x0)...(x-xw_,)_	„
iV(xN-X1)..-(xN-X?/_1)
Степень многочлена (3) не превосходит и в узлах интерполяции он
принимает заданные значения	=	из (1), следовательно,
является интерполяционным, после раскрытия скобок и приведения
подобных членов совпадает с многочленом, полученным при решении
задачи (1).
П.З. Кусочно-многочленная интерполяция. Сплайны
В ряде случаев алгебраические интерполяционные многочлены
высокого порядка оказываются неэффективными, так как нахождение
их значений требует большого объема вычислений, а точность полу-
ченного результата может оказаться недостаточной. Тогда применяют
кусочно многочленную интерполяцию. Этот подход состоит в сле-
из котопь^сто”"5' интер"0ЛЯции разбиваются на группы, для каждой
торого меньше°сте ЯрСВ°И Интерполяционный многочлен, степень ко-
го по всей совокупностиТлоРвПТЦИ°ННОГО Многочлена’ построенно-
многочленов для ппик/ ’ Исп0ЛьзУется каждый из построенных
исю"ой™й
узлов имеет с другой не бопер Л п°строен- Если любая из двух групп
«я»»™ bS",",6™"«'«ж точки, то «м miX
представляющего собой об-крпи^^^ интеРполяции многочлена,
пионных многочленов - каждого^ ВСех посгР°енных интерполя-
‘ огочлен и есть кусочно-многочленный°еМ 0Трезке- Объединенный
।	пяленный интерполяюг. Но при этом в
42
узлах общих для разных групп может нарушаться свойство гладкости
объединенного многочлена, так как первые производные справа и
слева для этих точек могут иметь разные значения. В ряде задач это
недопустимо. Например, для построения уравнения состояния веще-
ства по экспериментальным таблицам давления как функции плотно-
сти и удельной внутренней энергии (температуры) необходимо, чтобы
в измеренных точках было непрерывно не только давление, но и част-
ные производные от него по плотности и удельной внутренней энер-
гии.
Удовлетворить этому требованию возможно, если ввести в рас-
смотрение для каждой из групп узлов алгебраические многочлены
более высокой степени, чем интерполяционные. В результате решение
задачи о приближении функции становится не единственным, а полу-
ченную свободу в выборе многочлена можно использовать для прида-
ния объединенному многочлену свойства гладкости в узлах задания
исходных данных. Кусочно-многочленный интерполянт, обладающий
не только свойством непрерывности, но и гладкости заданного поряд-
ка в общих для разных групп узлов точках, называется сплайном. До-
пустим, что каждая группа содержит только по два ближайших узла
сетки и нужно предложить алгебраический многочлен степени не вы-
ше третьей, который обеспечивал бы гладкое продолжение на концах
каждого из отрезков. Предположим также, что известен интерполяци-
онный многочлен, построенный по всей совокупности узлов сетки, и
требуется, чтобы значения первых производных сплайна в узлах в
точности равнялись первым производным от интерполяционного мно-
гочлена в этих же точках. Тогда задача построения элемента сплайна
(X) = (а0), + (с?) ), X + (а2), х2 + («3). х3
между х, и х|+1 узлами сетки, где / принимает все значения от 0 до У
- 1, сводится к решению системы линейных уравнений четвертого
порядка:
(а3). х3 + (а2), х,2 + (а,), х, + (а0), =	(х,) = f,
(«3 1	+ (а2 ), 41 + (Я1 ), Х/+< + (а0 ), = PN 41 ) = Z+1 >
З(а3 \ X2 + 2(а2), х, + (а,), = P’N (х, )=P'N,= Р^.,.
З(а3)(x/+l + 2(а? х/+( + (<?] )г — Pf) (xj+j) — Р^_м •
Её решение имеет следующий вид:
43
Выполнение лабораторной работа позволяет овладеть одним из
способов аналитического приближения сеточных функций, основан-
ного на решении задачи интерполяции; приобрести необходимые на-
выки алгоритмической реализации теоретических положений и про-
вести их практическую проверку с использованием ЭВМ; накопить
опыт в использовании алгоритмических языков высокого уровня и
отладке программ.
В процессе работы предполагается решение следующих задач:
1) построение согласно (п. 2) данного пособия одного интер-
2)
поляшюнного многочлена по всей совокупности заданных
узлов интерполяции и соответствующих значений сеточной
функции в них;
нахождение кубического сплайна, используя ранее постро-
енный интерполяционный многочлен и руководствуясь по-
ложениями (п. 3} данного пособия;
44
3) вычисление приближенных значений функции в точках
принадлежащих отрезку интерполяции, с применением
найденного кубического сплайна.
При разработке алгоритма необходимо предусмотреть после оп-
ределения интерполяционного многочлена в одной из известных форм
записи приведение его к алгебраическому степенному виду и выводу
на экран монитора значений коэффициентов при соответствующих
степенях многочлена; у построенного кубического сплайна для каж-
дого элемента S3i (х) между двумя узлами сетки X, и Х(+1 вывести
на • печать	четыре 1 коэффициента
(а0 );, (л()(, (а 2 )(, («з)(,1 = 0 + N -1, а также возможность вво-
да значения X е [х0, XN ], вычисление значения кубического сплайна
в этой точке и вывода результата на экран монитора. При отладке
программы в качестве тестовой задачи можно воспользоваться сеточ-
ной функцией, построенной по тем же узлам, что и в предложенном
задании, взяв в качестве чисел, задающих функцию, значения любого
алгебраического многочлена, степень которого меньше степени иско-
мого интерполяционного многочлена. В этом случае вычисленные с
помощью интерполяционного многочлена значения в любой точке,
принадлежащей отрезку интерполяции, должны совпадать с величи-
нами, полученными с помощью выбранного алгебраического много-
члена. Для проверки правильности нахождения элементов сплайна
53/(х), < = О-ьДГ-1 можно сравнить значения функции и первой
производной от интерполяционного многочлена с соответствующими
величинами, полученными с использованием элементов сплайна, в
узлах сетки.
СПИСОК ЛИТЕРАТУРЫ
1.	Рябенький В. С. Введение в вычислительную математику. —
М.: Наука, 1994. —335 с.
2.	Сборник задач по методам вычислений / Под ред. П.И. Мона-
стырного—М.: Физматлит, 1994. —320 с.
3.	Лабораторный практикум по курсу Основы вычислительной
математики.— М.: МЗ-Пресс, 2001. — 192 с.
ПРИЛОЖЕНИЕ
построить кубический сплайн по заданной таблице
значений функции
							
Л*	0.17453	0,5236	0,87267	1,22173	1,5708	1,91986	2,26893 |
	0,00162	0,00252	0,00498	0,0129	0,03964	0,1207	0,34097 1
Задание № 2							
X	0,17453	0,5236	0,87267	1,22173	1,5708	1,91986	2,26893
	25 10'5	0,00116	0,00361	0,01162	0,03827	0,11933	0,33960
Задание № 3							
X	0,17453	0,5236	0,87267	1,22173	1,5708	1,91986	/ /
3’	38-Ю'6	0,00052	0,00254	0,01013	0,03636	0,11699	
Задание № 4							 1-'
X	0,17453	0,5236	0,87267	1,22173	1,5708	‘ 1	
IL	6-I0"6	0,00023	0,00173	0,00854	0,03373		
46
	ГЗадание №81 МКХ						
I7	0,52360	0,87267	1,22173	1,57080	1,91986		
	0,00000	0,00017	0,00199	0,01282	0,05744		
Задание № 9							
X	0,87267	1,22173	1,57080	1,91986	2,26893	2,61799	2,9670
У	0,00082	0,01039	0,07037	0,32762	1,18669	3,59003	9,4835
Задание № 10
X	0,87267	1,22173	1,57080	1,91986	2,26893	2,61799	2,96706
У	0,00123	0,01343	0,08411	0,37341	1,31146	3,88447	10,10742
Задание № 11
X	0,87267	1,22173	1,57080	1,91986	2,26893	2,61799
У	0,00161	0,01550	0,09139	0,39329	1,35729	3,97819
Задание № 12
X	0,87267	1,22173	1,57080	1,91986	2,26893
У	0,00196	0,01686	0,09511	0,40157	1,37344
Задание № 13
X	0,17453	0,52360	0.87267	1,22173	1,57080	1,91986	2,26893
У	3-10^	0,00018	0,00227	0,01770	0,09688	0,40481	1,37878
Задание № 14
X	0,17453	0,52360	0,87267	1,22173	1,57080	1,91986	2,26893
У	12-10"6	0,00026	0,00250	0,01815	0,09763	0,40593	1,38035
Задание № 15
X	1,22173	1,57080	1,91986	2,26893	2,61799	2,96706
У	0,01834	0,09787	0,40623	1,38070	4,01752	10,3450
47

____Задание № 16
х 11,22173 11.57080 I i f
Jfi01'849'~
flo/wG ,1.91986 \2,26893 [2.6П99
|oTQ98Q210.40638~ 1,38085 4,01768

48
РАЗДЕЛ III
ОБЫКНОВЕННЫЕ ДИФФЕРЕНЦИАЛЬНЫЕ
УРАВНЕНИЯ
Лабораторная работа № 3
Метод Рунге-Кутты решения задачи Коши для
обыкновенных дифференциальных уравнений первого
порядка
Введение
При изучении самых разнообразных явлений окружающего
мира, имеющих отношение как к точным, так и к гуманитарным нау-
кам, исследователи сталкиваются в ряде случаев с тем, что функцио-
нальные зависимости между величинами находятся из уравнений, в
которых присутствуют производные от искомых функций. Наиболее
простыми среди них являются те, что содержат только производные
первого порядка и могут быть записаны в виде
~ ,
ах
где у - искомая функция, х - независимая переменная, f(x,y) - не-
прерывная функция от х и у. Однако получить аналитическое решение
этого уравнения для достаточно произвольной функции f не удается и
только для некоторых частных случаев, с которыми можно ознакомить-
ся в справочной литературе, например, [1], [2], решение задачи пред-
ставляется в виде конкретной функции.
В связи с быстрым развитием электронной вычислительной
техники в последние десятилетия появилась возможность использовать
приближенные математические методы для решения подобного рода
задач. Один из таких подходов называется методом Рунге-Кутты, и он
объединяет целую группу модификаций, связанных способом их полу-
чения.
(1)
in 1 Задача Коши для обыкновенного
J£p—него уравнения первого порядка
ие"р'р“в,,° «*ренцирхемая
функция Ж у) И задано уравнение
^с=^Х,У^'
Согласно теореме существования и единственности для любой точки
(х0,у0)е G найдется решение у = у(х), удовлетворяющее условию
у(хо) = У°> и оно будет единственным. Задача для уравнения (1) с
начальным условием у(х0) = у ° (задача Коши) состоит в нахождении
функции у(х), обращающей и уравнение (I), и начальное условие в
тождество. Допустим, что значения, которые принимает независимое
переменное х, принадлежат интервалу , XN ), и запишем задачу
Коши:
у(хо} = У°,
%o{X{XN.
(2)

Разобьём отрезок [X0,Y 1	„
1 и’ wi иа 7V частей так, что
”+> хп~"„, и = 0-г?У —1 Впя
ности, рассмотрим случай, - ьнеищем> не ограничивая обш-
h» = h = XN-X0= const,
= *0 + nh, n = oTv v,
сеткой к..	узлами сетки, а множество
™>«-E«, кад тачке
«>»тор« ч»с» „ =	поставлено .
1 называть сети "	‘ > То множество значений
ФУваде», опрелелекной па сетке
.-const, “а" Г6"™"е р™°мерное, т.е. все
Назовем точки x
точек {x I _ .
соответствие i
будем
М*.
50
Ш.2. Метод Рунге-Кутты
Приведенный ниже способ построения сеточной функции по диф-
ференциальной задаче (2) называется методом Рунге-Кутты.
Положим, что величина сеточной функции в узле Хо = Хо равня-
ется начальному условию из (2) у0 = у°, а её значения в следующих
узлах сетки П = 1 -г N будем находить последовательно по формуле
yn+i=y(xn+h) = y(xn)+hXbift x^c^.y^h^jfjl,
'= К .	7=1 J
п = 0 N -1,	(3)
где	;	
Л =Ж + c2h,y„ + ha2Jl) = f2(x2,z2),
/з = Ж +c3h,y„ +ha3lfl + ha32f2) = f3(x3,z3),
ft =f(*„+ СЛ У„ + ha^f + haj2 + haj3) = /4 (x4 ,z4),
= xn + csh; z,=yn+ h£ asjfj ;
7=1
а^,Ьп cs - константы, которые подлежат определению.
Предполагается, что у(хп ) уже известно и имеет конкретное значение.
При нахождении ух в качестве у(хп), п = 0 берется
Яхо) = У о =/ ит-д-
Рассмотрим функцию
51
р(Л) = А +h)-y(xn)~h^bif(xi,zi)
(4)
V	-- hs+l; 0 < 0 < 1
~~&Т	(* + D!
и определим порядок точности (сходимости) метода.
Определение. Разностный метод Рунге-Кутты (РК), задаваемый соот-
ношениями (3), имеет q порядок точности (сходимости), если
р(0) = (0) =... =	(0) = о, а(р(^}(0) * 0,q <S.
Покажем, как определяются коэффициенты , Ct на примере ме-
тода РК четвертого порядка точности.
III.3. Вывод формул метода Рунге-Кутты
4-го порядка
Вычислим первые четыре производные no h от функции (p(h)
(4), полагая S = 4 и дифференцируя её как сложную функцию:
....... А / ч 4
из
(5)
-Н±ь^'
IV	4
i=l	Wh L-ii\Ji )hhhh •
В формулы (5) входят производные пт
фвертого порядка включительно. J"*0™* ФУНВДИИ + Д°
У) из правой части дифференциальна ЧСРе3 пРоизв°Днь1е
y)f = f,y'a = f+j-lf Реального уравнения (2):
+2 f" Гх Г" У2	..	/ ,
(6)
52
Ж-fL + з/'/+з/- f+/• f+/; (з/;+5/;/)+
+r„ (37;/+4/;/2)+w,*(/;)’ (/;+/;/) 
Также потребуются первые три производные по h от функций
/,(x„z,), 1 = 1,23,4.
Первые производные:
(/1)л =®’(Л)л ~ kfik с2 +(Л)12 a21fl ’
(/з )л = (/з )х3 сз + (/з )z3 аз1 fi + ^32/2 + ha32 (/2 )л J;
(/Л =(/4)х4 с4 +
'43/3 + ha К (f2 )л + ha43 J.
~(/ )х2х2 С2 + 2 (/2 )х2г2
+ (/4 )г4 a41/l
Вторые производные:
(/i) м
[f)hh =(-А)хЛ С3+2(Л)хл С3 a3\f\+a32f2+‘
+ 032/2 + h°32 {/г )* 1
'А
+ (/к''
(Л)„ =(Л).
с. X
'41У|
г4г(А)л +а4з(Л)* 11
+(Д к (2а42 (л)/ + 2я43 (Л)/ + 4(Л )м + «4, (Л L
Третьи производные:
53
(Н =0-
V1 /ААА	nt
UL =(Аи>Ш« ‘
ЧЛи"'"^’-
и)»' =(/>)„.,., с>+3<ли.с'-
+(А)„.
+з(Л),
+’(/.)вд
'*А
flf
«2 «21/ +
:3 «31/ + a^fz + М (/ )а
2
3
lit Г	('J
|«Э1/1+а32/2 + Аа32(А)А
, С3 2^32 (Л\ +^аз1{/2\
«31/1 + «32/2 + ^«32 (/>)/,
х МЛ)* +Ч2(/2)Ы1 ]+
+(/зК ^2^ +ha32[f2)M ,
(/“)1АА ~(f<)xAxAl. С4 +
X
+
«Р'Л +a^ + a^ + han{f2)h + ha43(f3\
°,,fl +^л+»вЛ+Ч2(Л);
2
'Л
3
54
+3UK-. с<{2а«(Л)* +2a<>W/+*p.!(/.)/+M/X"]j
+
х[ч2(Л)А +MA +h a42(f2)J + a43(f3)J J + (/4)/x
x' 3a42(/2)AA +3fl«(/3)hh +h a42(f2\
\hh +а4з(Л)ш
Для получения разностного метода РК первого порядка согласно
определению из пункта 1П.2 нужно в первое выражение из (5) подставить
h = 0 и приравнять ^(0) = 0, учитывая, что у'х = f из (6). Так как
это равенство должно выполняться для произвольной непрерывно диф-
ференцируемой функции f(x,y) из правой части уравнения (1), необ-
ходимо, чтобы коэффициенты , / = 1-5-4 удовлетворяли соотноше-
нию
Л, +b2 + b3+b4 =1 (/).	(7)
Условие (7) должно быть выполнено для всех четырех стадийных
S = 4 схем РК первого порядка из (3). Задавая три из четырех коэф-
фициентов i = 1 + 4, получим семейство схем РК первого поряд-
ка, но можно использовать имеющуюся свободу выбора коэффициентов
для построения схем более высокого порядка точности. Потребуем,
чтобы (p"hh (0) = О в (5) для любых fx, f' из второго уравнения (6):
2Z>2c2 + 2b3c3 + 2b4c4 = 1	(X)’	(8)
2b2a2i + 2b3(a3l + а32) + 2&4(а41 +а42 + я43)=1 (///).	(9)
Аналогично случаю схем (S = 4) РК первого порядка (3) все схемы
(5 = 4) РК второго порядка точности (3) должны удовлетворять усло-
виям (7), (8) и (9), что дает возможность предложить семейство схем
второго порядка точности, т.к. число определяемых коэффициентов
(йу,bi,ct,i = 1 -г4,J = 1 ч-/ — 1) превышает количество уравнений.
Для коэффициентов схем (S = 4) РК третьего порядка (3), которые
удовлетворяют условию (0) = 0, должны выполняться равенства
55
з+Ь4с42)л1	(f**)'	.
г (а.. + а„) + бЛ («4! + ал +	)] = 1 (ЛЛ*
3[М2.+Ьз(а31+^)
6 Ь3^2аП + ^д(С2а42 + С3Й43
Г/ -л2 ,1
'32
ДЮ)
(И)
(12)
(13)
(14)
Требуя для выполнения (р1^ (0) — 0 того, чтобы коэффициенты при
произвольных
f". r f f* fl. f* fl. f" f. ff'f- f'f'f1-' '
J na> JxxyJ’JxyyJ ’ J yyyJ ’ J xyJ x> J xyJ yJ ' •> ?yJ yJ ’ . (|
fj';,
обращались в ноль, приходим к ещё одиннадцати уравнениям для ко-
эффициентов, которым удовлетворяют четырехстадийные (S = 4) ме-
тоды РК четвертого порядка (3):
Фз(«3! + а32)а32С2 + ^4 («41 + «42 + «4зХ«42«2 + «43^3 )] = 1
	(21)
3 (2^3 («31 + «32 ) «32«21 + ^3«32«21 + 2^4 («41	«42 + «43 ) Х
ХЕ«42«2! + «43 («31 + «32 )] + ^4 [«42«J + «43 («31 + «32 /]} = 1
,	(/;/;/2),	(22)
12^53а32с2 +^(ад2+«43^)] = 1	(/«У/),	(23)
2464«43«32С2 = 1	'	((Л')2Л')	(24)
24Z>4«43«32«2i = 1	(/;)7-	.	(25)
Выполнение равенств (7) - (25) необходимо для того, чтобы четырех-
стадийные методы Рунге-Кутты (РК), задаваемые формулами (3), име-
ли четвертый порядок точности. Девятнадцать уравнений (7) - (25) со-
держат только четырнадцать неизвестных коэффициентов
(«о	1' ~ 1 = 4, j = 1 + i — 1), и, казалось бы, система уравнений
переопределена, и у неё не может существовать решений, определен-
ных в обычном смысле. Но оказывается, что часть уравнений в этой
системе линейно зависима и число линейно независимых уравнений
меньше числа искомых коэффициентов, т.е. можно предложить целое
семейство 4-стадийных методов РК четвертого порядка.
Ш.4. Таблицы Бутчера
Следуя [3], удобно представлять модификации методов РК в виде
таблиц Бутчера следующего вида: !	’
А .  .. Л Л " j. i А • „ . 5 '
>$7

где коэффициенты atj fb^C,, i — 1 -J- S’, j — 1 : i 1 те же, что и в
выражениях (3) для 5-стадийного метода РК. Приведем несколько наи-
более употребительных методов РК 2,3 и 4 - порядков точности, пред-
ставленных в виде таблиц Бутчера.
2-стадийные методы РК второго порядка точности
1) Модифицированный метод Эйлера
=	=/[хп+|>Л +|/Д j/„+l =уп +hf2
о Го
1/2 1/2 О
 0 1
2) Метод Эйлера с пересчетом
Л	+h’y°+ hf^	= Уп + ~(/1 + Л)
о ГГ	2
J_ 1 о
3’"М"Й"“'М"ВД“РКЛ'™«»»Ртеа„,И()ста
И Метод Хойна
( ' 2	2
Л = / xn+~h;yn +~hfi
к	3	3
Л+1 = Уп ^(/i +3Л)
0	0
1/3	1/3 0
2/3	0 2/3 0
	1/4 0 3/4
2)
Уп+t = К + 4 (Z + 2А + /з)
0 2/3 2/3	0 2/3 0 -1/3 1 0
	1/4 2/4 1/4
Л =f{xH,y„);f2 = /^f„ + ^,уп + j/J;
/з = f{xa +h,y„~hfl+ 2hf2);
Уп+i ~ У„+ . (/i + 4Л + /з)
6
0 1/2 1	0 1/2 0 -12	0
	1/6 4/6 1/6
59
4-стадийиый метод РК четвертого порядка точности
,	/ А А
)> fi ~ /1-у ’ У» а ]’
Классический
/,=Ж’^
А =/^.+|л +	-/(*« +А>Л +АЛ)
Л.,=Л+^и+2/> + 2Л+Л)
о	.
' ----1
0 1/2 1/2 1	0 1/2 0 0	1/2 0 0	0	10
	1/6 2/6 2/6 1/6
III.5.	Оценка погрешности численного решения
Важнейшим вопросом численного интегрирования задачи Коши
для обыкновенного дифференциального уравнения является оценка
погрешности полученного решения. Непосредственное вычисление в
узлах сетки нормы разности двух решений, аналитического и численно-
го, не всегда возможно, так как дифференциальная задача в общем слу-
чае не интегрируется до конца. Получить приближенное представление
о допущенной погрешности можно, если известен порядок точности
используемого численного метода. В этом случае сначала решают раз-
ностную задачу на сетке с шагом 2h, а затем повторно — на сетке с ша-
гом h. Обозначим сеточную функцию, полученную на сетке с шагом 2h,
через у , а на сетке h - След решения дифференциальной
задачи на сетке 2h пусть будет у*, а след р<А) на сетке 2h - yw . То-
гда справедливы два равенства:
^=/+0(2^,
/',? = /+СМ,	J (26^
нения (26) второе,Чим«м Выбранного метОДа. Вычитая из первого урав-
или	V V
60
к „У
2к
-y{h}
-1
(27)
Вводя в конечномерном пространстве, соответствующем сетке 2h, одну
из возможных норм, получаем оценку
I <„>	к2** -/! .
И >' ;* ’	~£-	,2S>
где 8 -	точность. Если условие (28) не выполняется, в два
раза увеличиваем число узлов сетки одновременно уменьшая шаг сетки,
и снова проверяем условие (28) на новой сетке. Так поступаем до тех
пор, пока не будет выполнена оценка (28).
III.6.	Лабораторная работа на ЭВМ
При выполнении лабораторной работы происходит освоение одно-
го из методов численного решения задачи Коши для обыкновенного
дифференциального уравнения первого порядка, приобретается практи-
ческий опыт решения математических задач с использованием ЭВМ и
способов оценки погрешности, совершенствуются навыки применения
алгоритмических языков высокого уровня для написания программ и их
отладки, осуществляется проверка теоретических положений численно-
го анализа на конкретных примерах.
В лабораторной работе предлагается решить следующие задачи:
1.	Разработать алгоритм решения задачи Коши для обыкно-
венного дифференциального уравнения одним из методов
Рунге-Кутты и написать программу, реализующую его на
одном из алгоритмических языков высокого уровня;
2.	Отладить программу и провести расчеты до достижения
заданной точности.
При разработке алгоритма необходимо предусмотреть возмож-
ность:
1)	изменения области интегрирования;
2)	задания других начальных условий;
3)	проведения расчетов на последовательно удваиваемых сетках;
4)	вывода в одиннадцати равноудаленных точках, расположенных
равномерно на отрезке интегрирования, значений сеточных
функций, полученных на сетках с шагами 2h и Л, а также их
разности в тех же точках.
61
СПИСОК ЛИТЕРАТУРЫ
Камке Э. Справочник по обыкновенным дифференци-
альным уравнениям. -М.: Наука, 1971. - 576 с.
Зайцев В. Ф.. Полунин А.Д. Справочник по обыкновен-
ным дифференциальным уравнениям. - М.: Физматлит,
2001.-576 с.
Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновен-
ных дифференциальных уравнений / Нежесткие задачи.
-М.: Мир, 1990. - 512 с.
ПРИЛОЖЕНИЕ
Найти указанным методом решение задачи Коши иа заданном интервале
в одиннадцати равноудаленных точках с указанной точностью £.
Задание № 1
Метод Рунге-Кутты 4 порядка точности; X G (1, 2); £ — 10 ~4
2х2уу'х + у2 = 2х3 +х2,
Х1)=1-	"
Задание № 2
Метод Рунге-Кутты 3 порядка; X е (0,1), £ = 10-4 ,
(х2у - 1)ух + ху2 -1 = 0,
Я0) = 0.
Задание № 3
Метод Эйлера с пересчетом X е (1,1,1), £ = 10-4,
(бху + х2 + з)ух + Зу2 + 2ху + 2х = 0,
Х1) = -1-
Задание № 4
Модифицированный метод Эйлера; X е (1, 2), £ = 10-4 ,
(ху-х2)ух +у2 -Зху-2х2 =0,
XD = 1 + V2.	.
Задание № 5
Метод Эйлера с пересчетом х G (1,1,2), £ = 10"4 ,
х(2х!у  1пО) +1)>>, =2у,
XD = 1.
143
Задание № 6
Модифицированный метод Эйлера; х € (О, 1), £ = 104,
у'х ~ху2 ~3ху — О,
Х0)=-з.
Задание №7
Метод Рунге-Кутты 3 порядка точности; х е (1,1,2), S — 10-4 ,
2ххУх ~У2 +5х = 0,1
XD = 1. J
Задание №8	' Л 1 ’
Метод Рунге-Кутты 4 порядка точности;	.
хе(1,2), £ =
хух + хуг =	'
У(1) = 2.	-
Задание№ 9 \ У
Метод Рунге-Кутты 3 порядка точности; X € (1,2), £ = 10 >
х(2х-1)УЛ + / -(4х + 1)у + 4х = 03	; , ч.
У(1) = 2,	J. „•
, ... . . ... ;
Задание № 10
Метод Рунге-Кутты 4 порядка точности; X € (1,2), ’£• = 10-4 ,
2х1у'х-2у2 -Зху + 2х = 0Д
У0) = 0,5.	(
'	 -V , 7 Л
Задание Ns 11	'	1	> г
Модифицированный метод Эйлера; х <= (1,2), е = 10"4 ,
г	' "р;
64
ху'х - ху2 - (2х2 + l)j- х3 - О,
^а)=-з.	.
‘ Задание №12
Метод Эйлера с пересчетом х е (1,2), £ = Ю-4,
х2(ух+ у2)+4ху+ 2 = 0,.
Я1) = -1-	, / -
Задание № 13
Метод Рунге-Кутты порядка точности; х & (1,2), £ = 10-4,
К«,  6'-*2К=*/|
/ х1)=1’5- J
• И .	1
Задание № 14
Метод Рунге-Кутты 3 порядка точности; х е (1, 2), £ — 10"4,
уух + у2 + 4х(х +1) - О,
я1)=12.
Задание №15
Модифицированный метод Эйлера; X G (1,2), £ = 10"4,
 Х3^-X4j2+х2^ + 20 = 0,
Я1) = 4.
Задание № 16
Метод Эйлера с пересчетом х е (2, 3), £’ = 10"4,
X2 (х - 1)л - У2 - х(х - 2)у = 0,1
’	И2) = 4.	J
6$
раздел IV
второго ПОРЯДКА
Лабораторная работа № 4
Метод прогонки
Введение
Большинство физических процессов, протекающих в окружаю-
щем человека материальном мире, имеет нестационарный диссипативный
характер, но в ряде случаев оказывается, что условия, задаваемые на гра-
ницах областей и определяющие решение, не зависят от времени. Тогда
можно ожидать, что через некоторый конечный промежуток времени
процесс станет стационарным и частные производные по времени от ис-
комых функций обратятся в ноль. Например, если рассматриваются про-
цессы теплопроводности или диффузии в конечном объеме вещества,
граничные условия не зависят от времени, а исследователей интересуют
установившиеся стационарные распределения температуры или концен-
трации в изучаемом пространстве, то искомое решение будет являться
функцией от пространственных координат. Допустим к тому же, что за-
дача обладает свойством изотропности по двум ортогональным простран-
ственным направлениям из трех. Тогда математическая модель представ-
ляет собой краевую задачу для обыкновенного дифференциального урав-
нения второго порядка. Остановимся более подробно на трех возможных
случаях постановки задачи [1].
IV.1. Третья краевая задача для стационарного
?пХеНИЯ теплопРов°Дности с непрерывными
коэффициентами
P«nwaWyS^	обыкн°венного диффе-
интервале (0,1) с точностью е ₽	даумя кРаевьми условиями на
66
(I)
-£(1)mx(1) = S2w(l) - e2,
где , §2, £1, £2 - const, k(x), q(x), f(x) - непрерывные ограниченные по-
ложительные на (О, I) функции. Построить аналитическое решение задачи
(1) для произвольных функций к(х), q(x),f(x) не удается, но если заменить
указанные функции константами, тогда можно решить задачу (1) анали-
тически в общем виде и полученное решение использовать в дальнейшем
для отладки численных программ. Поэтому прежде чем приступить к
численному решению задачи (I), сначала рассмотрим вспомогательную
модельную задачу с постоянными коэффициентами k,q,f на (0,1) и с
той же точностью е:
d г Г -I —	Z*
-г ... ах ах
Ц0)4(0) = 51и(0)-Е1, >
-jt(l)wi(l) = 52w(l)-e2,
(2)
где Sx, 82, е2 - const, Л = Л(0,5), <7 = ?(О,5^ f - /(0,5).
IV.1.1. Аналитическое решение модельной задачи
Общее решение дифференциального уравнения из (2) можно
представить в виде и = побщ + ычасгн > где иобш ~ °®щее решение одно-
родного уравнения, а «части _ какое-то частное решение неоднородно-
го дифференциального уравнения второго порядка (2).
Для нахождения общего решения однородного дифференциального урав-
нения второго порядка сделаем в нём подстановку w = е2' и после со-
кращения на придем к характеристическому уравнению:
Ы2-т=0.	(3)
Разрешая (3) относительно X, получим
67
Тогда общее решение однородного уравнения можно записать в виде
Их - Их
+С2£		(4)
где С, С2 - произвольные постоянные. То, что ыобщ из ” общее ре-
шение однородного уравнения, проверяется подстановкой его в это урав-
нение и получением тождества. В качестве частного решения неодно-
родного уравнения можно взять
и Л Л  ' • 
“части — >
q		
и тогда общее решение задачи (2) представляется в виде
Лх	f	.	'
w=C|6 1 +С2е 1 +4=.	(5)
После подстановки в краевые условия выражения (5) приходим к системе
из двух линейных уравнений относительно двух неизвестных постоянных
С, и С2.
(*Л] £])<?! +(^-^^2=5ji-fp
- (^+S2)e\ =s2’	<6)
Решая (6), получаем
С] ’Л^2^1/-у)Л+(кХ2.&])(Ъ2у.£2?) ‘
^г51Яп2+82;Л.^2.61ХГХ1+б2;/17:
^^2-^rn1+S2/i.fn2+S2Xt- Л;
F	 ; . и
и тем самым находим решение задачи (2), подставляя (7) в (5).
IV.1.2. Численное решение модельной задачи
Постановка разностной модельной задачи
Введем в области интегрирования на [0,1] равномерную сетку,
Выбрав в качестве узлов точки х, = lh, I = 0 4- L, где h = —. Если каж-
дому узлу сетки х,, l = 0 + L поставить в соответствие некоторое значе-
ние , I = 0 4- L, то сеточную функцию можно определить как упорядо-
ченное множество значений {/}, заданных на сетке {х,}, I = О-rZ,. Рас-
смотрим конечномерное пространство размерности £+1, элементами
которого будут сеточные функции {ft} = определенные на сетке
{х,} . Если в этом пространстве ввести норму как = max | / |, то
получим конечномерное нормированное пространство размерности £+1.
Для постановки разностной задачи поступим следующим образом: в каж-
дом узле сетки /	заменим вторую производную в дифференци-
альном уравнении (2) конечно-разностным отношением вида
и В результате получим уравнение, приближающее дифференциальное в
узле /:
—	—2,111 +Uj i	_	-------
A-J+l------1----L±-qU =-f-l = UL-l.	(8)
h2	1
Система уравнений (8) состоит из L - 1 линейного уравнения, относи-
тельно L + 1 неизвестного значения сеточной функции и(А) Для коррект-
ной постановки задачи воспользуемся двумя краевыми условиями из (2),
записав приближенные выражения для первых производных:
(9)
к- ~¥о *1’
Объединяя уравнения (8) и (9), получаем систему линейных
уравнений £ + 1 порядка относительно L +1 неизвестного значения
сеточной функции = {«/}, / =0+£-1, которую и будем называть
разностной модельной задачей.
Метод прогонки
Для решения разностной модельной задачи или, что то же самое,
системы линейных уравнений порядка £ + 1 применим метод прогонки,
обратив внимание на тот факт, что каждое уравнение системы содержит
не более трех ближайших значений искомой сеточной функции. Для
удобства последующих преобразований умножим уравнения (8) на Л2, а
каждое из двух уравнений (9) - на h , затем после перегруппировки чле-
нов уравнений введем новые обозначения:
a^k^b^-lk-qh^-c, = k-d, =-J/i2, z = 17£-l;
ao = к; Ьо = -fc - 5th; c0 = 0; d0 = ~£th;
aL=Q-,bL=-k-S2h-cL = k; d, = ~£2h
и перепишем систему (8), (9) в виде
a0ul+b0w0=rf0’
alUl+\	= d/# /	.
bLuL+cLuL-\=dL-
МСГ0Д ПЕИппяГТ0ИТ И3 ДВУХ ЭТапов: п₽ям°й и
(Ю)
обратной прогонки,
сначала первое уравнение (10) разрешают
Уволят обозначения
“O-'y-Ui+i
ьо *0
70
w0 и
a°=-^0 ;>5()=^0 ,U°=a0ul+^0-	(4)
Далее, используя метод индукции, докажем, что систему линейных урав-
нений (10), имеющую трехдиагональную матрицу, можно преобразовать
в систему линейных уравнений с двухдиагональной матрицей. Для дока-
зательства предположим, что I — 1 уравнение системы (10) уже представ-
лено в виде
«Z-l=aZ-lMZ+^Z-l	<12)
с известными коэффициентами и . Покажем, как к аналогично-
му виду приводится и Z - уравнение системы (10). Подставим в него вы-
ражение для из (12) и разрешим относительно ut:
и __ ai и .drcA-\
' bl+clal-l м bl+clal-l
Обозначим
Лт	d, —Cifl] ।
а1 = ~T-----’ Pl = A-----(14)
1 bl+clal-l ' bl+clal-\
и перепишем (13):
иГа1иМ+Рг	(15)
Значения at и Д находятся, так как коэффициенты a^b^c,,
dt задаются при постановке задачи (10), а значения at_} и счита-
ются известными по предположению индукции (12).
У читывая формулы (11) и (15), в которых первое уравнение сис-
темы (10) и Z - уравнение этой же системы приведены к виду (15) с из-
вестными коэффициентами at и /?,, согласно методу индукции делаем
вывод о том, что любое уравнение системы (10) приводится к виду (15).
Следовательно, это справедливо для L —1 уравнения системы (10), и оно
представимо в той же форме с известными коэффициентами и
А-. =
uL-\ =aL-\uL +Pl-V	(1б)
Подставим цЛ_! из (16) в последнее уравнение системы (10) и разрешим
относительно uL :
71
_dL clPl-V	,	(17)
w£ bL +cLaL_x
Нахождением u£ заканчивается прямая прогонка. ____
При обратной прогонке, используя формулы (15) для / - £ -1 - 0 и вы-
численное значение «£ (17), находят все W£, а следовательно, и се-
точную функцию , которая является решением разностной модель-
ной задачи (8), (9).
IV. 1.3. Численное решение задачи с переменными
коэффициентами
При решении задачи (1) с переменными коэффициентами внача-
ле вводят в области интегрирования сетку, например, разбивая отрезок
[0,1] равномерно на L частей так, что расстояние между соседними уз-
лами постоянно и равно h = 1/L. Совокупность точек
хо = 0; = х0 + lh,l = 14- L образует сетку. Подобно тому, как это
было сделано в nJV. 1.2, сформулируем разностную задачу, заменив во
внутренних узловых точках / = 1 ч- £ — 1 в дифференциальном урав-
нении (1) производные конечно-разностными отношениями и использо-
вав следующие обозначения’
^,±Л/2) = Кх,±1/2)^;±1/2;	'
лачи (1) приходим ксис,ек..РаЖеН-И^ ВДИ<^^еРенциальное УРавнение за-
L +1 "	™"“ЫХ УРа’“ШЙ "Ч*™ £ -1 относк-
°П”“,еив’«-очно» функции „,.( = 071,
72
применить на практике известные из курса обыкновенных дифференци-
альных уравнений приёмы нахождения аналитических решений для
уравнений с постоянными коэффициентами, проверить в численном экс-
перименте теоретическое положение о сходимости разностных решений к
следу аналитического решения в случае корректно поставленной конеч-
но-разностной задачи, углубить навыки алгоритмической реализации
сложных вычислительных процессов и провести их практическую про-
верку с использованием ЭВМ, накопить опыт в использовании алгорит-
мических языков высокого уровня и отладке программ. В процессе рабо-
ты предполагается решение следующих проблем:
1)	постановка и аналитическое решение модельной задачи;
2)	численное решение модельной задачи с заданной степе-
нью точности;
3)	численное решение задачи с переменными коэффициен-
тами и проверка сходимости на последовательно удваи-
ваемых сетках.
При разработке алгоритма необходимо предусмотреть
возможность проведения расчётов на последовательно увеличивающихся
сетках, одновременный вывод на печать в одиннадцати равноудалённых
точках с шагом h - 0,1 значений аналитического решения модельной
задачи, численного решения модельной задачи, полученного на сетке с
шагом hx кратным h (h = Nhx, N - целое число), их разности в тех же
узлах, решения задачи с переменными коэффициентами.
СПИСОК ЛИТЕРАТУРЫ
1. Трифонов Н.П., Пасхин Е.Н. Практикум работы на ЭВМ. -
М.: Наука, 1982. -288 с.
2. Самарский А. А., Николаев Е. С. Методы решения сеточных
уравнений. -М.: Наука, 1978. - 592 с.
95
ПРИЛОЖЕНИЕ
Вариант № 1
Найти решение краевой задачи для одномерного стационарного уравне-
ния теплопроводности
^[k(x)^]-q&)u=-f(x}
ах ах
в одиннадцати равноудалённых точках отрезка [0,1] с относительной
точностью 0,0001. Отладку программы произвести на модельной задаче с
постоянными коэффициентами.
Задание №1	'
Краевые условия задачи	' :
Wt(O) = o, 1	'
-	Л(1)Мх(1) = м(1), 
*	(*) = е1; q(x) = е’; /(х) = sin х.
Модельная задача
= 4e-q(x) = 4e;f(x) =sin 0,5.
Задание № 2
Краевые условия задачи
^(0)u/0) = u(0)s
-	fc0X(l) = w(l)-lJ :
’	= ex',q(x) = ex;f(x) =cosx
Модельная задача	J1 cosx-
ВД = Л;,(х)=7;;/(х)=соз015
Задание № j
Краевые условия задачи
к(0)их (0) = 100a(0) ,1
-*(!)<(!) = О, >
2
к(х) = х + l;g(x) = x+2;/(x)==cosx.
Модельная задача
к(х) = 1,25; q(x) = 2,5; f (х)=cos 0,5.
Задание №4
Краевые условия задачи
*(0Х(0) = и(0)-1,
-*(!)< (1) = «(!)-!,	,
2	Y
к(х) = sin х +1; q(x) = sin х; f (х) = е .
Модельная задача
к(х) = sin 2 0,5 +1; q(x) = sin 0,5; f (x) = -fe.
Задание№5	i -
Краевые условия задачи
*(0>’(0) = 0, 1	•
-*(!)<(!) = «(!),/
к(х) = х2 +1; ^(х) = х; /(х) = е~ х.
Модельная задача
к(х) = 1,25; <?(х) = 0,5; /(х) = 1/V^.
Задание № 6
Краевые условия задачи
к(0)их(0) = и(0),
-*(1)и;(1) = и(1)-1,
fc(x) = 1 + cos^ х; q(x) = 1; /(х) = sin х.
Модельная задача
к(х) = 1 + cos2 0,5; 9(х) = 1; /(х) = sin2 0.5.	..
97
Задание №7
Краевые условия задачи
Л(0Х(0) = и(0),
-fc(lX(l) = wOb .
к(х) = е^пх;д(х) = еСО5Х;Лх) = 1.
Модельная задача
^(х) = е“п0’5;^(х) = ес“0,5;/(х) = 1.
Задание №8
Краевые условия задачи
Л(0Х(0)=и(0),Г
-л(1)М;(1)=М(1), ’
£(х) = ecosx;<?(x) = es,nx;/(x) =1.
Модельная задача
&(x) = ecos0,5;g(x) = esin0’5;y(x)=l. . Л
Задание №9
Краевые условия задачи
к(О)ых(О) =0, |
к(х) = х + l;g(x) = ех	.
Модельная задача
к(х) = 1,5; q(x) =	/(х) = е А25.	'5
Задание № Ю
Краевые условия задачи
fc(0)<(0) = u(0)-O
- WxC) = w(l)-l, '
*(x)»l+sin2 x;g(x) =COSX;/(X) = ex ;
98
Модельная задача
к(х) = 1 + sin2 0,5;q(x) = cos 0,5;/(х) =
Задание № 11
Краевые условия задачи
Ц0)их(0) = 0,
-Л(1)< (!) = «(!),/
к (х) = ех; q(x) = 1 + х3; /(х) = е~х.
Модельная задача
£(х) =	= 1,125;/(х) =е-°’5.
Задание № 12
Краевые условия задачи
ад«;(0) = м(0)-1,
-Л(1)<(1) = м(1),
к(х) = ex;q(x) = ех;/(х) = cosx.
Модельная задача
к(х) = л/ё; q(x) = Je;f(x) = cos 0.5.
Задание № 13
Краевые условия задачи
цо>;(о)=и(о),
-А(1>х(1) = «(!),]
к(х) = ex;q(x) = ex;f(x) = е'х.
Модельная задача
к(х) = 4e;q(x) = 4e;f(x) =е~°'5.
Задание № 14	..
Краевые условия задачи	,,, ; .
Модельная задача	n f(x) =e'*'25 
'д’.- J-<-’ 1 '
Задание № 15	, j .
Краевые условия задачи	.	;	..	.	,,
Ц0К(0) = и(0),1	—
k(x) = e*;g(x) = ex;/(x) =smx. ,
Модельная задача	- ; . t -w« •
k(x)= 4e;q(x) = 4e-,f(x) =sm0.5.
Задание № 16
Краевые условия задачи	' ,	- - V
fc(0)ux(0) = 0, 'I	. Pt-в! >сн.^ • :
fc(x) = х +1;qr(x) = х2 + l;f(х)
Модельная задача	.
к(х) = 1,5;$(х) ^1,25;/(X) = l/Ve.
Вариант № 2
Найти решение краевой задачи для одномерного стационарного уравне-
ния теплопроводности
^5*W	~	=~ W
в одиннадцати равноудалённых точках отрезка [0,1] с относительной
точностью 0,0001. Отладку программы произвести на модельной'
lfifr
задаче с постоянными коэффициентами.
Задание № 1
Краевые условия задачи (условия периодичности)
w(0 + 0) = w(l - 0),
А(0 + ОХ (0 + 0) = Л(1 - 0)их (1 - 0),
к(х) = sin 2лх +1,5; q(x) = 1; f(x) = cos 2лх.
Модельная задача
к(х) = 1,5; q(x) = 1; /(х) = cos 2лх.
Задание №2
Краевые условия задачи (условия периодичности)
w(0 + 0) = w(l - 0),
А(0+0Х (0+0) = к(1 - 0)их (1 - 0),
к(х) = cos 2 2лх +1; q(x) = 1; f(x)=sin 2лх.
Модельная задача
к(х) = 2; q(x) = 1;/(х) =sin 2лх.
Задание № 3
Краевые условия задачи (условия периодичности)
w(0 + 0) = w(l-0),
Л(0+0)< (0+0) = Л(1 - OX (1 - 0),J
к(х) = cos 2 2лх +1; q(x) = sin 2лх +1,5; /(х)=sni 2ях.
Модельная задача
к(х) = 2; q(x) = 1,5; f (х) = sin 2лх.
Задание №4	"	'
Краевые условия задачи (условия периодичности)
и(0+0) = и(1-0),	р
*(0 + d)ux (0+0) = Л(1 - 0Х (1 - 0),
*(х) = е“2ж;«(х) = 1;/(х) = е”2“: 1
Модельная задача
к(х) = 1; q(x) = 1; /(х) = 1 + sin 2лх.
W
Задание №5	.
Краевые условия задачи (условия периодичности)
и(0+0) = и0~°)>
к(0 + 0)< (0+0) = А(1 - 0)< (1 - °)>.
k(x) = V,q(x) = e JW-e
Модельная задача
A(x) = 1;9(x) = 1;/(x) = 1 + cos2®c. .
Задание № 6
Краевые условия задачи (условия периодичности)
u(0+0) = u(l-0),
А(0+0)и j (0+0) = Л(1 - 0)< (1 - 0),
к(х) = е"“”2 2хх; q(x) = 1; f(x) = cos 2лх.
Модельная задача
А(х) = 1; q(x) -1; f(x) = cos 2лх.
Задание Nq 7
Краевые условия задачи (условия периодичности)
w(0+0) = «(1-0),
А(0+OX (0+0) = А(1 - O)uj (1 - 0),
fc(x) = cos3 2лх + 2; q(x) = esin 2ях; /(x) = 1.
Модельная задача
Л(х) = 3; g(x) = 1; /(х) = cos 2лх.
Задание № 8
Краевые условия задачи (условия периодичности)
w(0 + 0) = u(l-0),
к (0 + 0)ил (0 + 0) = k(\ - 0)u j (1-0),
,,	= ln(sin 2т+2); q(x) = 1- f_ i
Модельная задача	’/' ) Ь
*W = ln2;9(x) = l;/(x) = sin2w.
102
Задание № 9
Краевые условия задачи (условия периодичности)
w(0 + 0) = w(l-0),
Л(0 + 0)w j (0 + 0) = Л(1 - 0)w j (1 - 0),
к(х) = cos2 2лх +1; q(x) = ln(2 + sin 2лх);/(x) = cos2ax.
Модельная задача
k(x) = 2; q(x) = In 2; /(x) = cos 2лх.
Задание № 10
Краевые условия задачи (условия периодичности)
w(0 + 0) = w(l-O),
Л(0 + 0)w/0 + 0) = Л(1 - 0)w/l - 0),
A(x) = sin4 2лх +1; q(x) - cos3 2лх + 2;/(х) = e0082®.
Модельная задача
к{х) = 1; q(x) = \f(х) = cos 2лх +1.
Задание №11
Краевые условия задачи (условия периодичности)
w(0 + 0) = w(l - 0),
Л(0 + 0)w (0 + 0) = *(1 - 0)< (1 - 0),
А(х) = esin 2”; q(x) = 1; /(х) = ecos2m
Модельная задача
к(х) = 1; q(x) = 1; /(х) = cos 2лх.
Задание № 12
Краевые условия задачи (условия периодичности)
w(0 + 0) = и(1 - 0),
Л(0 + 0)М; (0 + 0) = Л(1 - 0)« х (1 - 0)J
к(х) = е-5Ш 2лх; q(x) = I; f (х) = cos 2лх.
Модельная задача
к(х) = 1; q(x) = 1; /(х) = cos 2ях +1.
103
Задание №13	?
Краевые условия задачи (условия периодичности)
и(0+0) = «(1-0)>	;,s
ЦО+0)«х (0+ 0) = *0 ~	~
к(х) = cos2 2т+1; q(x) = е*2тJ /W = ес“2ж,.,
Модельная задача	 <-ы
к(х) = 2; q(x) = 1; f(x) = cos2^
Задание№14	• . ’	..и»  >
Краевые условия задачи (условия периодичности)	.г, < .*> яир!
w(0+0) = w(l-0),	. . v f.b,
Ц0+0)ч^(0+0)=Ц1-0)ч^(1-0)>]	,	. >
к(х) = cos3 2лх+2; q(x) = sin 2т+ 1,5; /(х) = sin 2т,
Модельная задача
к (х) = 3;д(х) = 1,5; /(х) = sin 2т + 1.
Задание№ 15	, • .
Краевые условия задачи (условия периодичности) f;
w(0+0) = M(l-0),	, ]
fc(0+0)«; (0+0) = k(\ - 0)ux (1 - 0),
Л(х) = 1; q(x) = ln(2+sin 2лх); f(x) = cos 2m.
Модельная задача
^(*)=q(x) = In 2; /(у) = cos 2m.
Задание № 16
Краевые условия задачи (условия периодичности)
к(0+0) = и(1-0),
*(0+0)<(0 + Q) = *(!_())„ (1_0) ’
к{х) - sip 2дх +1;,^(х) cos3 2т+ 2' f(x) -1
Модельная задача	M-W-1-
er .«hi
>v
‘W-UsM-WW ..........
104
Вариант №3 с;	. s ,,
Найти решение краевой задачи для одномерного стационарного уравне-
ния теплопроводности
О) -^] ~ q(x)u = - Дх)
ах ах
в одиннадцати равноудалённых точках отрезка [0,1] с относительной
точностью 0,0001. Отладку программы произвести на модельной задаче с
постоянными коэффициентами.
Задание № 1
Краевые условия задачи и(0) - 0; м(1) = 1.	v.-j/ X
Дополнительные условия в точке разрыва :	>  j .-»
w(x0 ~0) = u(x0 +0),	-	’ h.-.-.хюГ.
Л(х0 - OX (х0 - 0) = k(xQ + 0)их (x0 + 0),J
x < x0 = 0,525; k(x) = x2 + 1;^(х) = е-х; Дх) = 1;
x > x0 = 0,525; k(x) = x; q(x) = e-x; Дх) = x3.
Модельная задача
=0,525;A(x)=A(^); <Хх)=<з(^);/(х)=/(дь).	“	/
   , i >.
Задание № 2
Краевые условия задачи и(0) = 1; w(l) = 0.	< v . . > л
Дополнительные условия в точке разрыва
w(x0 - 0) = м(х0 + 0),	|	' .хГ.
Л(х0 - 0Х (хо “ °) = к(хо +	(хо +°)>/
х<х0 = 1//л/2;Л(х) = х2 +0,5; <?(х) = 1; Дх) = 1;
х > х0 = 1/V2 ; к(х) = х2 + 0,5; q(x) =	;/(х) = cos х.
f '
Модельная задача	. z yj .zN
xi =]/>/2;А(х)=ад<Хх)=4(хь);/(х)=Л^Х
Задание № 3
Краевые условия задачи и(0) = 1; и(1) = 2.

Допоянюельиие условия	в точке разрыва	'	у	'
ф»-о>=«<1'«+®	[	I,
^x0-ox(x0-o)^(x0+o)B;(xfl+o),j	\
х С х0 = 1/75 ;*(х)= sin х2 +1; q(x) = х; f(x) = 1;
х>х0 =l/V5;jt(x) = sinx2+l;^(x) = x3;/(x) = x -1.
Модельная задача	(
^=^/5;#х)=Ш<Хх)=	Г ,„г
Задание №4	"	- i ..t. . . -f».<
Краевые условия задачи и(0) = 1; »(1) = 0. и >	>t. глгп . .	.
Дополнительные условия в точке разрыва >
и(хо~О) = и(хй+в), ,	1 ;-'
к(хй - 0)их (х0 - 0) = к(х0 + 0)их (х0 + 0),J
х<х0 = 1/л/2;Л(х) = е’х;(7(х) = х2;/(х) = 1; t
x>x0=l/V2;<:(x) = l;i7(x) = e'?;/(x) = cps.Xj<	j.:
Модельная задача	. Л
^=V>/2;^(x)=^);g(x)= ^);/(х)=Я^),
..	. лчЛ
3adanueNs5	j ;-;  <<
Краевые условия задачи и(0) = 1; м(1) = 0>. 	>. ; нч: ’
Дополнительные условия в точке разрыва ’ '
«(хо-0) = м(хо+0),	()
Чх0 -О)ах(х0 -О) = А(хо +О)их(хо +0), ,f >
< х0 = 1/з/2; Л(х) = 1; ^х) = I;/(х) = ;
^>х0=1/Л;к(х) = е“"-;17(х) = 2;/(х) = е\
Модельная задача
ч ^»V^;^x)=^Vr<?(x)=^);/(x)=/(xi)),
'"'fl'. Hi’nt.j
106;
Задание №6	'	<
Краевые условия задачи w(0) = 0; «(!) = !.	'	; < .
Дополнительные условия в точке разрыва	• <
н(х0-О) = и(хо+0),	1
к(х0 - 0)и j (х0 - 0) = к(х0 '+ 0)их (х0 + 0),/
х < х0 = 0,525; Л(х) = е~^; q(x) = х2; f(x) = sin х;
х > х0 = 0,525; Л(х) = х; </(х) = х 2; /(х) = sin х.
Модельная задача
=0,525;А(х)=^); q(x)= <Х^);/(х)=/(^).
Задание № 7	‘	"
Краевые условия задачи м(0) = 0; м(1) = 1. к " ;	! '
Дополнительные условия в точке разрыва
и(х0-0) = w(x0+0),
А(х0 - 0)w j (х0 - 0) = Цх0 + 0)wx (хв + 0),
X < х0 =1/Л;Л(х) = г‘пх;д(х) = 2;/(х) = е\ ,
х>х0 = 1/V2;Цх) = 1;?(х) = 1;/(х) = е.*	-  .у х «
Модельная задача	*	<
= 1/Т2;Л(х)=А(^); <j(x)= <Х^);/(х)=/(^).
Задание№8	!
Краевые условия задачи и(0) = 1; и(1) = 0.		
Дополнительные условия в точке разрыва	S1	а
«(*о-°) = и(хо+°)»	|	,4
Л(х0 -0)wx(x0 -0) = Л(х0 +0)их(х0 +OXJ
х < х0 = 1/-Ло ; к(х) = cosх + 2; <?(х) = 1;/(х) = 2х;
х>х0 = l/VF0;£(x) = sinx + 2;«7(x)= 1;/(х)=?0.	;• s > с .
Модельная задача
3b=VViO;A(x)=^); q(x)= ^);/(x)=/(Ao>
107

„ ,п ,	*Mww
Краевые условия задачи и(0) = 0; и(1) - • .
дополнительные условия в точке разрыва	, ?1,
и(Х0-0) = «(*0 +°)>	: I
к(х0 -ОК(хо -0) =	+0КОо +P)J
х < х0 = 1/7з ;*(х) = 1;?(х) = х2;Л*) = х2 ."J
X > ХО = 1/^; ад =<х;-?(х)=х;/(х)=1.
Модельная задача
^=1/7з;Л(х)=А(хь);4(х)= <Х^);Ях)=/&)*.•« ввн.ии{."оМ
Задание Ns 10
Краевые условия задачи и(0) = 0; u(l) = 1.	“ <?/
Дополнительные условия в точке разрыва,. . '»№;4
и(хй-0>) = и(хй+0),	,-4	;»ii-OiirR
k(x0 -0)мх(хо -0) = к(х0 + 0)w/xe +0),J	' ,
х<х0 = 1/72;к(х) = х2 +О,5;<7(х) = е-Х ;/(x) = cosx;
х > х0 = 1/72; к(х) = X2 + 0,5; q(x) = I; /(х) = 1, J
Модельная задача	.;	т
^=№*(*>=*<^Wx)= falfW-fW' '
Задание № 11
Краевые условия задачи «(0) -1; и(1) = 0.	>> g./
Дополнительные условия в точке разрыва  ,	-<rs' ;н<  >v i
“(xo~0) = w(x0+0),	"I
i(x0 -O)uI(xo -0) = £(x0 +O)uj(xo +0),j lx
x < x0 = 0,125; A(x) = x+1; 9(x) = x2; /(x) = cos
'	x>xo=0>125^(x) = x2;q(x) = x2;/(x) = cosx.
Модельная задала
>i=ai25;A(x)=^); <Kj)= ^);/(x)=/(jtbK

Задание №12
Краевые условия задачи и(0) = 1,' u(Y) = I.	'	'
Дополнительные условия в точке разрыва "
w(x0 -°) = у0о +°)>	1
*(х0 - 0)wx (х0 - 0) = к(х0 + 0)wx (х0 + 0),J
!х<х0 = 1/V3 ;i(x) = l;^(x) = x2;/(x) = sinx; 5
х > х0 = 1/л/з; fc(x) = ecosx; q(x) = х2; /(х) = sin х.
Модельная задача
^=0,125;*(x)=^(V, Ф)=
Задание №13	’
Краевые условия задачи м(0) = 0; и(1) = 1.’	’’ •
Дополнительные условия в точке разрыва '
и(х0 -О) = а(хо +0),	I
iЛ(х0-О)мх(хо-О) = Л(хо + О)мх(хо + 0),
' х<х0 =0,525;fc(x) = x + I;<7(x) = e_JC;/(x) = 1;
х > х0 = 0,525; к(х) - х; q(x) = е"х ;f(x) = х3.
Модельная задача
=0,525;Л(х)=^); <j(x)=
Задание № 14
Краевые условия задачи w(0) = 0; w(l) = I.
Дополнительные условия в точке разрыва
и(*о-°) = ы(хо+0)’
Л(х0 - 0)нх (х0 - 0) = k(xQ + 0)wx (х0 + 0),
х < х0 = 0,525; к(х) = х; q(x) = е~х; /(х) = х3;
х>х0 = 0,525; £(х) = х2 +1;^(х) = е-х;/(х)=1.
Модельная задача
х^ =0,S25;k(x) =к(х0); <j(x) = <Х^);/(х) = Д^).

Задание №15
Краевые условия задачи и(0) = 0; “(1) = ® 
Дополнительные условия в точке разрыва	,,	.
w(x0 -О) = и(хо +0),
к(х0 - 0)<(х0 - 0) = к(х0 + 0)их (х0 + О),
х < х0 = 0,125; к(х) = х +1; q(x) = е ~х; /(х) = cos х;
х > х0 = 0,125;А(х) = х2;^(х) = е-1;/(х) = 1.
Модельная задача
=О,125;Ш)=*С0; (&)= ^b);/(x)=/W-
Задание № 16
Краевые условия задачи п(0) = 2; «(1) = 1.
Дополнительные условия в точке разрыва
u(xo-O) = w(x0+0),
^(хо ~®)их(хй-0) = к(хй +0)их(х0 +0),
х<х0 =1/7з ,k(x) = e'x;q(x) = x3,f(x) = x2 -1;
х>х0 = 1/73;А(х) = е“х;<7(х) = х;/(х) = 1.
Модельная задача
*=VV3;*(x)=^);^= ^);/(х)=Д^).
ПО
РАЗДЕЛ V
УРАВНЕНИЯ С ЧАСТНЫМИ
ПРОИЗВОДНЫМИ ГИПЕРБОЛИЧЕСКОГО
ТИПА
Лабораторная работа № 5
Методы решения уравнения переноса
Введение
В ряде разделов физики для описания функциональных зависи-
мостей между величинами используются уравнения с частными про-
изводными первого порядка, например, система уравнений Эйлера - в
механике сплошных сред, уравнения переноса излучения - в атомной
физике и т.д. Наиболее актуальные и сложные проблемы, как правило,
имеют нелинейный характер и при их решении возникают определен-
ные сложности, преодолеть которые в некоторых конкретных случаях
удается с помошью методов вычислительной математики.
Для задач с коэффициентами при частных производных, яв-
ляющимися функциями только от независимых переменных, широко
применяются аналитические методы решения, которые обычно изу-
чаются в курсах по обыкновенным дифференциальным уравнениям.
Аналитические подходы оказываются чрезвычайно полезными также
при выборе модельных задач, так как на их основе происходит тести-
рование численных алгоритмов и отладка программ, написанных на
одном из языков программирования высокого уровня.
V.I. Математическая постановка задач
В пособии представлены два варианта постановки смешанной
задачи для уравнения переноса с неотрицательным и неположитель-
ным коэффициентами. Это разделение обусловлено особенностями
задания граничных условий и способами численного решения соот-
ветствующих разностных задач [1]. Для удобства изложения материа-
ла в качестве области интегрирования выбран единичный квадрат, что
не ограничивает общности формулировки начальных и граничных
условий для произвольного прямоугольника.
Первый вариант постановки задач имеет следующий вид:
111
—+a(x,l)— = b(x,t); c?(x,/)>0,
dt - dx
м(х,О) = ^(х); 0<х<1,	*
	i^O,/) = !/(/); 0£t<l,
где a(x,t), b(x,t), 0>(x), {/(/) - непрерывные функции своих ар-
гументов в области интегрирования, обладающие требуемой степенью
гладкости.
Другой случай записывается в виде
А«	Л.
— + о(х,/)— = b(x,t); a(x,t)<Q,
dt	dx
w(x,0) = p(x); 0<х<1,
u(l,t) = p(t); 0£t<l.
(2)
V.2. Аналитическое решение задач
Для построения точного решения задач (1), (2) необходимо рас-
смотреть вспомогательную систему уравнений (см. [1]):
Решение первого уравнения из (3), которое содержит диффе-
ренциалы независимых переменных, задает в плоскости (x,t) линию,
принятую называть характеристической. Второе уравнение с диффе-
ренциалом искомой функции и определяет изменение ее вдоль харак-
теристики. Решение первого уравнения	(характеристика) и
второго - v2(x,t,u) представляют собой два независимых первых
интеграла системы (3).
Тогда общее решение уравнений (1), (2) есть произвольная
функция от этих первых интегралов:
F(y\(x,t),v2{xsiy) - 0 или v2(x,t,u) = F,	(x,/)) . ,
112

Если удается разрешить последнее равенство относительно не-
известной функции li(x,t), то получаем общее решение в общем ви-
fle:a(V) = ^(x^vi(x’O)-
Для выделения из общего решения единственного частного не-
обходимо дополнительно учесть начальные и граничные условия.
Сначала удовлетворим начальному условию, полагая в первых
интегралах t = 0:
vI(x,0) = vI°,
v2(x, 0,u) = v2.
Разрешая (4) относительно X, и , приходим к выражениям для
х и искомой функции и через значения первых интегралов:
* = <(«), '
u = W°2(v?,v°2).
Далее подставим равенства (5) в начальные условия (1), (2), за-
менив предварительно значения первых интегралов с индексом 0 на
соответствующие первые интегралы системы (3):
(4)
(5)
[V, (х, t), v2 (x, t, u)] = <p { w°[v, (x, t), v2 (x.t, w)]} .
(6)
Полученное выражение (6) и есть искомое решение задачи в
области, определяемой характеристиками, выпущенными из концов
отрезка х е [0,1] при / = 0 и продолженными до пересечения с гра-
ницами области интегрирования.
Для нахождения решения в оставшейся части области интегри-
рования, где оно зависит от граничных условий, нужно положить в
первых интегралах х0 = 0 для задачи (1) или х0 = 1 для задачи (2):
v,(x0,r) = vI0,
v2(x0,/,w) = v20.
Выражая t и и из (7), приходим к следующим зависимостям:
U = ^2o(VlO’Vx)-
‘ Меняя в (8) v10 на v,(x,/) и v20 на v2(x,t,u) и подставляя в
граничные условия задач (1), (2), получаем соотношение
113
(7)
(8)
ди , .ди ,. .	. .. _
—	+ a(x,t)— = b(x,t); a(x,t)>0,
dt - дх
и(х, 0) = р(.г); 0 < х < 1,
i/(0,7) = $/(/); 0^1 <1,
где а(х,/), b(x,t), <р(х), lp(f) - непрерывные функции своих ар-
гументов в области интегрирования, обладающие требуемой степенью
гладкости.
Другой случай записывается в виде
ди , , ди	, . п
—	+ а(х, t) — = b(x, Г); а(х, t) £ 0,
dt дх
и(х,0) = р(х); 0^х<1,
w(l,/) = ^(/); 0£/£1.
(2)
V	.2. Аналитическое решение задач
Для построения точного решения задач (1), (2) необходимо рас-
смотреть вспомогательную систему уравнений (см. [1]);
Решение первого уравнения из (3), которое содержит диффе-
ренциалы независимых переменных, задает в плоскости (х,/) линию,
принятую называть характеристической. Второе уравнение с диффе-
ренциалом искомой функции « определяет изменение ее вдоль харак-
теристики. Решение первого уравнения vz(x,/) (характеристика) и
второго - v2(x,/,m) представляют собой два независимых первых
интеграла системы (3).
Тогда общее решение уравнений (1), (2) есть произвольная
функция от этих первых интегралов:
^(vi(^,O>v2(x,f)) = 0или v2(x,/,«) = ГДуДх,/)) .
m
I
Если удается разрешить последнее равенство относительно не-
известной функции u(x,t), то получаем общее решение в общем ви-
де: w(x,/) = F2(x,/,v1(x,/)).
Для выделения из общего решения единственного частного не-
обходимо дополнительно учесть начальные и граничные условия.
Сначала удовлетворим начальному условию, полагая в первых
интегралах t — 0:
vi(x,0) = v1°,
v2(x,0,w) = v2. J	(4)
Разрешая (4) относительно x, и , приходим к выражениям для
х и искомой функции и через значения первых интегралов:
1
W = W°(V,°X)- J
Далее подставим равенства (5) в начальные условия (1), (2), за-
менив предварительно значения первых интегралов с индексом 0 на
соответствующие первые интегралы системы (3):
w'2[v1(x7),v2(x,/,!/)] = <p{wl°[v1(x,/))v2(x7,w)]}.	(6)
Полученное выражение (6) и есть искомое решение задачи в
области, определяемой характеристиками, выпущенными из концов
отрезка х е [0,1] при t = 0 и продолженными до пересечения с гра-
ницами области интегрирования.
Для нахождения решения в оставшейся части области интегри-
рования, где оно зависит от граничных условий, нужно положить в
первых интегралах х0 = 0 для задачи (1) или Хо = 1 для задачи (2):
V	1 (Xq>0 ~ V|0,
V	2(x0,r,w) = v20. J
Выражая t и и из (7), приходим к следующим зависимостям:
^ = WlC>(VlO'V2o)’ 1	(8)
w = w2o(vlo,v20). J
Меняя в (8) v(0 на vf(x,/) и v20 на v2(x,t,и) и подставляя в
граничные условия задач (1), (2), получаем соотношение
113
которое и определяет решение, зависящее от граничных условий, в
оставшейся части области интегрирования.
Таким образом, выражения (6) и (9) полностью определяют ис-
комое решение во всей области интегрирования.
Рассмотрим конкретный пример.
Пример 1
х2
и(х,0) = х+~; 0<х<1,
(Ю)

Запишем уравнения, из которых определяются первые интегра-
лы:
OD
х + е
dx = — .	(12)
х
Уравнение (И) является обыкновенным дифференциальным
уравнением первого порядка, линейным относительно x(f) , решение
которого можно получить в явном виде:
v((x,Z) = С° = хе"' -1 = const .
Следовательно, в области
D° = { (x,Z): 0 < t < 1; 0 х < 1; te' < искомая функция
u(x,f) зависит только от начальных условий, т.к. характеристика,
проходящая через начало координат, есть кривая х = 1е' . Напротив,
для области PQ = {(xj):0^r	x^te'} определяю-
щими являются граничные условия. Интегрируя (12), находим, второй
х2
v2(x,u) = u-— = Co - const.
Теперь для нахождения решения в D° вычисляем первые инте-
гралы при t = 0:
vf = X,
о X2
V, -и-----,
2
о о
выражаем X, и через V] , v2:
- ° ш 
и V* +1 2 J
(13)
и подставляем (13) в начальное условие задачи (10), заменив на
v,(x,0 и v° на v2(x,u) :
и - — + (хе"' - Г)2 - = хе~' -t + ~(хе~' -I)2.	(14)
2	2	2
Приводя подобные члены в (14), получаем решение в D°:
х2
и =— + хе -t.
2
Для определения решения в Do подставляем в интегралы
(7) х = 0:
У»=«. .
находим t и U :
' = -vio>
M=V2O’ .
(15>
Затем в (15) меняем vI0 на v}(x,t) и v20 на v2(x,u) и под-
ставляем в граничное условие задачи (10):
115
«•Gii:
х2 t	(16)
И"Т=
Ра,^» <1« «и-™"»	”МУ,““ "С'°“УЮ ФУ”"Ю
«(х,/) В области Da :
х2	i
и^ — + хе -t.
Следовательно, аналитическим решением задачи (10) в области
po=|(x^):0<Z<l;0<x^l} является функция ,
2
V3. Разностные схемы для численного решения
смешанной задачи уравнения переноса
V3.1. Построение разностных схем
Прежде чем приступать к разработке разностной схемы, необ-
ходимо в области интегрирования задачи ввести сетку. Это означает
разбиение всей области на некоторое число подобластей, которые все
пронумеровываются. Например, в задачах (1), (2) отрезок X G [0,1]
разбивается на L частей так, что Xz = lh', Lh=A, I = 0,Z, а отрезок
t e [0,1] - иа N частей: t" = tn, zN = 1, n = 0, W. Пусть X изме-
няется по оси абсцисс, a t - по оси ординат. Тогда пересечение I -
вертикальной и и - горизонтальной прямых определяет точку с ко-
ординатами (x/,t ), которую будем называть узлом сетки или узло-
вой точкой. Множество узловых точек, принадлежащих области ин-
тегрирования, образует сетку {х„Г }^л=0. Если каждому узлу сет-
ки (х,,Г ) поставить в соответствие некоторое число и”, то множест-
f 1L N
во всех значений [uz jz=0 п=0 будет определять сеточную функцию.
«	“ На ДаНН°Й Сетке сегочных Фмий образу-
даумя элементами 'Т°^гранство’ 601111 °но линейно и между любыми
ДВУМЯ элементами его определено понятие расстояния или метрики.
116
I
Рассмотрим один из способов построения разностной схемы иа
заданной сетке, выбрав в качестве примера задачу (1). Заменим част-
ные производные в узлах сетки конечно-разностными отношениями
следующего вида:
du(Xl,t") (ди\я (.v и^-иЧ
---  =   = ( и, ), « —-------------
dt	v 'Л т
/ = 2,Дл»0,ЛГ-1;
8“(Jt-,~) l'g»Y , ,, з»,--4<+<г
fa (fa J,	2h
l = ZL,n = 0,N-l.
Тогда в каждом из указанных узлов сетки вместо дифференци-
ального уравнения можно записать соответствующее разностное
уравнение:
. .з<-4<,+< i=^l
т ' 2h	' « = 0,jV-1
где a” =a(x/(Z") , Ь” =b^xl,t"'j .
(18)
Если к ним добавить выражения, приближающие начальные и
граничные условия задачи (1), то получится система линейных урав-
нений (18), (19):
u,° =<p, =<р(х,), Z = O,Z,
, \	---- (19)
“о ~ V = V(*" )> п = 1- N-
Число уравнений в этой системе равняется
(Z — 1)7У + Z +1 + N = Л7, + Z +1 = Z(N +1) +1 , а число неизвест-
ных значений сеточной функции - (Z + 1)(Л^ + 1) . Для корректной
постановки разностной задачи требуется добавить еще N уравнений.
Что это за уравнения и не будут ли они противоречить уже записан-
ным?
Прежде чем ответить на эти вопросы, исследуем уравнения
(18), (19) на локальную аппроксимацию, подставив след аналитиче-
/г	J nji.A'
ского решения Цы], j о вместо искомого решения {и, j/=0 я=0 и
вычитая из левой части уравнений правую:
П7
Множество значений получившейся сеточной функции называют ло-
кальнои невязкой q
Так как след определен не только в узлах сетки, но и в окрест-
ностях этих точек, то, предполагая требуемую степень гладкости ана-
литического решения, разложим значения следа в ряд Тейлора отно-
сительно тех узловых точек, в которых разностные уравнения при-
ближают дифференциальное:
2
[«L =[«]"-[<]"А+[<]" у+2(^3).	-	(2D
[м1',-2=[мГ-кГ2/'+[с]/в2й2+о(йз):
Затем подставим выражения (21) в (20), приведем подобные члены и
сократим на общие множители:
sf" -[“<]" +Mj+O(r2) + a” [M;J + О(/г2)-Ь/п,
i=2,L;n = 0,N-l;Sf/°=0, / = о7;^/о" =0, n = ij/.
(22)
В первой строке (22) сумма членов
!оч^сё:китаОРЯет' * ^альн
В результате получаем, что
М +а”[игх$ -Ь” = 0, так как
ому уравнению в каждой узловой
Sf” =Qfr + h2), l = 2,L, n = O,N-\,
5/o’=0, w =
8f,°=0, l = O^L.
Вводя в конечномерном подпространстве правых частей (23) норму,
приходим к тому, что из (23) следует ||<^(А)|| < Стт + Chh2, где Сг,
Ch - постоянные, не зависящие от г и h.
Разница в порядках аппроксимации по т и по Л сохранится и
при исследовании сходимости, что приведет к дополнительным объе-
мам вычислений на ЭВМ для обеспечения заданной точности. Если
повысить порядок аппроксимации до второго по г , то этого удастся
избежать.
Рассмотрим один из возможных способов повышения порядка
сходимости разностного метода. Для этого вначале выразим вторую
частную производную по времени через вторую частную производ-
ную по пространству, продифференцировав основное уравнение из (1)
один раз по времени, а другой - по пространству:
д2и , .д2и da(x,t)du db(x,t)
dt2 ' ' dtdx dt dx dt
d2u , .d2u da(x,t)du db(x,t)
dxdt ' dx2 dx dx dx
Умножим второе уравнение из (24) на a{x,t), вычтем из первого
уравнения и разрешим относительно второй производной по времени:
б2 и	2/	\д2и
-т~г = а (x,t)—-
dt2	v	'dx2
(24)
dx
db(x,t)
dt - \ / qx
Если в (25) заменить вторую производную по х, используя формулу
численного дифференцирования
f d2»Y	rf ~2ц/-| + Ц/-2
(S’:'	^СЬС2
а первую производную -
dt ' dx
db(x,t) z
J____1____L _ Zl( V t
(25)
h2
119
duX

2h
Г I'm In
И вычесть из левой части уравнений (18) выражение ~1ы« ]/ » то по-
лучим	и
“Г1-"” /„?“*
-2-------(а,)
)' ;
(26)
'=2’£;в=0’ЛГ-1-
„При исследовании на локальную аппроксимацию уравнений
(26), (19) находим, что норма невязки удовлетворяет оценке
|5/w|^Crr2 + Chh2, если Г порядка h или, наоборот, h поряд-
ка т .
Возвращаясь к вопросу о задании дополнительных уравнений
для корректной постановки разностной задачи, будем исходить из
того, что они не должны ухудшить порядок сходимости, который по-
лучен для уравнений (26), (19). Поэтому наиболее естественным спо-
собом является использование разложения значений следа
[wf, n=),N в ряд Тейлора до второго порядка по h включительно
относительно граничных точек (х0 ,t"):
1иГ=[«1+[<]ой + Щу+0(Л3)	(27)
с дальнейшей заменой частных производных по пространственной
переменной на частные производные по времени в (27). Значения
изаестны из граничных условий (1). Выражения для полу-
чаются из дифференциального уравнения (1):
« = Ш	<28>
во	'
120
гдеМ=(^)Л0> n = l,N.
Несколько сложнее находится выражение для jj. Восполь-
зуемся равенством (25), но только разрешим его относительно част-
ной производной по пространству второго порядка:
1
д2и
a2(x,t)\dt2
да(х,/)
dt
dZ>(x,r)l д2и
~ дх2’
da(x,t) du db(x,f)
дх
дх dt
(29)
' ' дх
заменим первую производную по X на соотношение, аналогичное
(28) из дифференциального уравнения (1), и подставим в (29). Окон-
чательно получаем
ГлЛ» 1	/ .. \п Г/ • и/ I Гтп / •	С’ " '
'о
(30)
Подставляя выражения (28), (30) в (27), получим недостающие N
уравнений для корректной постановки разностной задачи.
Выпишем разностную схему, которая по построению должна
аппроксимировать дифференциальную задачу на решении U со вто-
рым порядком и по h, и по т :
иГ] ~и1
т
+\а1
f['
2 и1? -2иТ ,+и"
| J---Ц—Ь2г4.
2h2
.	~4-\+v”~2 _
x"Jj 2h
Н31)
- а"
~Ь1
/ = 2,Д» = 0Л-1,
uj=q>lt l = 0j; ufi=pn, n = lJV,
121
(31’)
Аналогично разностной схеме (31), (31’) второго порядка ап-
проксимации на решении и по h и по Г строятся методы и более
высокого порядка точности: третьего, четвертого и т.д.
V.3.2. Исследование устойчивости разностных схем
спектральным признаком
Для применения спектрального признака при изучении устой-
чивости разностной схемы (31) нужно в основной подсистеме уравне-
ний при I = 2, L, п = О, N — 1 осуществить подстановку
«Г =(«),"+д;,	(32)
где (и)" - точное решение системы уравнений (31), а Д” означает
погрешность вычислений в соответствующих узлах сетки. В результа-
те все слагаемые, содержащие точное решение разностной задачи
(31), сократятся с правой частью уравнений, т.к. это решение обраща-
ет каждое из разностных уравнений (31) в тождество. В итоге полу-
чим подсистему уравнений относительно Д", I = П = О JV — 1 •
/ = 2,£, n = \, N.
га
При фиксированном времени /" функцию Д”, I = О, L из (33)
с добавленными граничными значениями Д"о, Д" можно рассматри-
вать как зависящую только от переменной х и представить в виде
конечного ряда Фурье:
a;=Z^'.	«и)
£=-/»
где X* - амплитуда к -гармоники в момент времени t”, а макси-
мальный номер гармоники т определяется из выражения
2т + 1 = L, L - натуральное нечетное число, соответствующее
наибольшему номеру узла сетки по h . Разложение (34) справедливо,
если сетка введена на отрезке х' ё [0,2л£/(£ + 1)].При помощи ли-
нейного преобразования х' = 2л£х/(£ +1) отрезок X е [0,1] пере-
водится в х' е [0,2я£/(£ + 1)].
Тогда (34) переписывается в виде
д; = £х;<?Z+I-	<35>
Л=-/я
Аналогичное представление сеточной функции Д"Я, l~Q,L
можно выполнить и для момента времени / , но амплитуды у гар-
моник будут иметь другие значения Х^+1:
л1*1
дг' = £	i+1-	(36)
Подставим выражения (35), (36) в уравнения (33) н вынесем
общий множитель за скобки:
1 <
>23
Jiflsi упрощения дальнейших выкладок воспользуемся принци-
пом замороженных коэффициентов и будем считать, что коэффициент
переноса а" является некоторой постоянной величиной и не зависит
ни от X, ни от /, т.е. а" = а, (о, =0,	= 0. Обозначим через
= Д*, ак = Ink/iL + 1). Умножим (37) на какую-то гармо-
нику е,р1, р~ -т,т, просуммируем по всем I -Q,L и поменяем
в двойной сумме порядок суммирования:
Л»
Хч
1-1	21 - 2е“'а* + е’,2“‘	3 - 4e"to‘ + е-,2в* )
г	2/i2	2й J
^>=0.	.	.
Iх-
(3&)
1=0
ОД * -р,
(L + X),k = -p.
(39)
у /k+p)S
Рассмотрим подробно сумму
1=0
fk+p^ -Л_________—
Принимая во внимание (39), получаем, что в сумме по к в (38)
остается только одно слагаемое к = —р , т.е.
4,-1	+еаа" 3-4е^-'„ ....
—-------а1------------------г + а----------------— - v, (40)
. t	2/1	2/i
для Vp, р = — т,тно в общем случае т ~	, L —> оо, следо-
вательно, т —> оо, В дальнейшем поэтому будем полагать, что (40)
выполняется для \/а = а_р е(-аэ, + со). Величина показывает,
во сколько раз меняется амплитуда к -гармоники при переходе с П
временного слоя на и +1, и явно не зависит от номера слоя. Поэтому,
согласно спектральному признаку устойчивости, чтобы погрешность
Д всегда оставалась ограниченной при стремлении t —> со или
Й -> 0, необходимо потребовать выполнения для V а е (— оо, + оо)
Ща)|<1.	(41)
124
Разрешим (40) относительно X = Х_р , обозначая а = а_р :
X =! - -^(З - 4е~'° + е-'^ ) + £^(j _ 2е~'а + е-па ).
2hx	' 2й2 v	>
Из (41) следует, что верно неравенство
|х|2 = XX < 1,	(42)
где X - комплексно-сопряженное к X выражение:
Х = 1-—(3~4е‘а+ei2a) + ^~-(l-2e'a+е‘2а) =
2/Л	1 2h2 k	>
= (1 -30 + 202) + 4р(1 - 0)е'°-0(1 - 20)е/2“,
где /3 =--. С учетом введенных обозначений неравенство (42) при-
2h
нимает следующий вид:
[(1 - 0)(1 - 20) + 40(1 -р)е-*-0(1 - 20)e~/2ajx
х[(1 - 0)(1 - 20) + 40(1 -0)е'“ - 0(1 - 20) е'2а] < 1.
Перемножая, приводя подобные члены, находим, что
l-60(l-0)(l-20)2 + 80(l-0)(l-20)2cosa-
-20(1-0)(1-20)2cos2a<l.
Сокращая единицы в левой и правой частях неравенства (43),
деля его на 20(1-20)2 .получаем
-(1 - 0)(3 - 4cosa + 2cos2a) < 0
или, преобразуя,
2(1-0)(1-cosa)2 >0 .	(44)
Сокращая в (44) на 2(1 - cosa)2 , имеем
0=“^1-
Неравенство (45) означает, что спектральное условие устойчи-
вости (41) будет выполнено при любом а , если
, 2h
т<— .	(46)
f'T' -  • '	а	 . :
125
Соблюдение неравенства (46) является необходимым условием
устойчивости разностной схемы (31), н оно должно в точности вы-
полняться при построении сетки, так как заранее неизвестно, по каким
именно гармоникам будет раскладываться случайная ошибка вычис-
лений Д" на каждом из временных слоев.
Когда коэффициент переноса a(x,t) является функцией неза-
висимых переменных X, t, в качестве а выбирают максимальное по
I значение а” на данном временном слое п: а = ШЭХ а” и меняют
- ,
Т от слоя к слою: г. <---—, ,	,,
• Л	и * .	।	Il } I
таха, •	4-	!
л
< В этом состоит основная идея принципа -замороженных коэф-
фициентов.	'	;	'
V33. Точность полученных численных решении. ,,,
Из общей теории численных методов [2-3] известно, что для
сходимости решений, получаемых на последовательно увеличиваю-
щихся сетках, разностные схемы должны: во-первых, аппроксимиро-
вать исходную дифференциальную задачу на решении и; во-вторых,
быть устойчивыми. Исследование аппроксимации в большинстве слу-
чаев не вызывает никаких дополнительных трудностей, а определение
условий устойчивости даже для задач с переменными коэффициента-
ми требует дополнительные предположения, например, принцип за-
мороженных коэффициентов, как это было сделано ранее в разделе 2,
не говоря уже о нелинейных уравнениях или разрывных решениях.
Применение же при расчетах разностных схем, не удовлетворяющих
условиям устойчивости, может приводить к непредсказуемым послед-
ствиям, когда при увеличении сетки сначала наблюдается улучшение
точности, а затем, несмотря на убывание невязкн, погрешность нарас-
тает. В таких ситуациях бывает практически невозможно оценить
близость полученного результата к следу аналитического решения,
отображенного на конечномерное пространство, если это решение не
найдено. Использование таких расчетов для анализа и объяснения
физических явлений сопряжено с риском возникновения больших
просчетов и ошибок. С другой стороны, сходимость численных реше-'
нии по сетке до шага определенной величины h, а также порядок
«мелениях ^°ЭВМИ СХеМЫ М0ЖН° проверить непосредственно при
М, последовательно увеличивая сетку в опреде-
<
ленное число раз и оценивая норму разности полученных решений в
выбранном конечномерном пространстве.
Пусть вначале численный расчет выполнен на сетке с макси-
мальным расстоянием между узлами, равным hn, где п - некоторое
натуральное число больше единицы. Затем вычисления были повто-
рены, но шаг сетки уменьшен в П раз. Допустим, что из исследований
аппроксимации известно, что данная разностная схема имеет к-
порядок аппроксимации по h. Отобразим второе решение на конеч-
номерное пространство, в котором решалась первая задача, взяв в ка-
честве отображения значения сеточной функции в соответствующих
(л)
узлах и обозначив его как '. Первое численное решение запишем в
виде y{nh), а след аналитического решения задачи в том же простран-
стве- через [у1,А. Тогда -	-
i;	1 L Jnh v	(47)

где предполагается, что постоянная С практически не меняется при
изменении шага сетки в И раз. Тогда, вычитая почленно в выражени-
ях (47) из первого уравнения второе, находим, что величина погреш-
ности приблизительно равна
	и -i
’ ’ Следовательно, уменьшение шага сетки в п раз должно приво-
дить к улучшению точности в пк раз для схем с к -порядком аппрок-
симации по Л и устойчивых. Если же схема неустойчива или про-
грамма не полностью отлажена, то, начиная с вполне конкретного h,
сходимость перестает наблюдаться и нарушается соответствие между
Порядками аппроксимации и сходимости, в чем можно убедиться,
проводя вычисления на последовательно увеличивающихся сетках.
Этот прием позволяет практическими расчетами на ЭВМ кон-
тролировать правильность теоретических исследований устойчивости
используемых разностных схем или качество отладки применяемых
программ. Наиболее эффективным данный подход оказывается в тех
случаях, когда модельная задача имеет аналитическое решение и для
любого конечномерного пространства может быть построен след.
л
V.4. Лабораторная работа нй ЭВМ	, ,	«	>
нЯЧЯтъ выполнение лабораторной работы рекомендуется с на-
Иачать выполнпредварительно ознакомившись с
хождения аналита гдавЫ	’пособия и, если необходимо, с до-
по^^ной литературой. В качестве конечномерного пространства,
на котором будет происходить сравнение получаемых на последова-
^вьТудашшаемых сетках численных решений, предлагается вы,
брать	____ ,
Dh = {(x,t = l'):xl~hl,hL-l,l-0,LiL = 10] . т.е. равномер-
ную по х сетку с одиннадцатью узлами, расположенными на времен-
ном слое t = Т = 1. Следующим шагом является построение следа,
т.е. вычисление значений аналитического решения в одиннадцати,
равноудаленных точках при t = 1. Удобно представлять результаты
работы в виде таблицы, в первом столбце которой приведены значе-
ния координаты х узлов сетки, во втором столбце - аналитическое
решение (след), в третьем столбце - отображение численного реше-
ния, полученного на кратной сетке на выбранное конечномерное про-
странство, в четвертом - модуль разности между следом и отобра-
женным решением. Максимальное по модулю число в четвертом
столбце характеризует норму, введенную как максимум сеточной
функции по всем узловым точкам даного конечномерного простран-
ства.
В приложении приведены три варианта разностных схем, отли-
чающихся друг от друга порядком сходимости и содержащие по ше-
стнадцать заданий. Наиболее простыми являются задачи из первого
варианта. Для их решения в каждом задании приведена разностная
схема, которая после аппроксимации начальных и граничных усло-
вии, заданных при постановке дифференциальной задачи, позволяет
^и«°вать алгоритм послойного перехода по временным слоям, на-
°Г0’ Н “ ?*ебует привлечени* никаких дополнитель-
ных разностных уравнении.
Дания^рХма СЛ°ЖНЫМИ 0Казываю™ задачи из второго за-
рой порядок по х в то 7 Т0Мэ ЧТ° разностнь1е Уравнения имеют вто-
ко первый. Поэтому требуетсяТ ДИ(^еренцнальное Уравнение толь-
коррекгной постановки разностной*™* дополнительноГо условия для
в п.У.3.1. (см. вывод формул (31))	например> как эт° сделано
Еще более сложными являются примеры из третьего варинта,
так как для корректной постановки разностных задач требуется зада-
ние двух дополнительных граничных условий. Для этого можно будет
воспользоваться тем приемом, который был использован в n.V.3.1.
при выводе разностной схемы второго порядка аппроксимации и по
Г , и по h , но при этом учесть большее число членов разложения в
ряд Тейлора значений следа в приграничной области.
Запишем выражение, аналогичное (27), оставив в рассмотрении
члены до третьего порядка по h включительно:
ь2	>
Для нахождения [м' , можно воспользоваться форму-
лами (28), (30). Третью частную производную по X будем выражать
через третью частную производную по времени, последовательно вы-
полняя следующие преобразования: во-первых, вычислим частную
производную по времени от левой и правой частей (25):
, д3и	2 33и	_ дад2и	Г да	да~\д2и	’
—г = а-------2а—-у- ——а—	—
dt3	dtdx2 dt дх2	dt	дх dtdx
г	,	,	(49)
d2a	д2а	да да ди	д2Ь _ д2Ь дадЬ
dt2	а dtdx	dt дх дх	dt2	dtdx dt	дх
—	—*	.
(3&)
д26
Затем из второго уравнения системы (24) определим
,j д2и . ( \д2и да ди db
(i lt , dxdt Х’ дх2 дх дх дх
и продифференцируем (50) один раз по X;
д3и , лд3и -дад2и д2а ди
dx2dt a Xf дх3 дх дх2 дх2 дх + дх2 ’
Подставляя выражения (50) и (51), умноженное на a2(x,f), в (49),
получаем '	‘	~	
(Г:. Го/И
•0HWI.I.!	-:»Z 03,tsq .'Rl	s 7. «;•?» .-ЫН'.ОлГ.^' ьЫНТЬаИ
-.ч-г- C-UT-ial	•I?) KWHjUb-•:’.»№ >	» ДОГ • л	.
-аг.'Яйаоьок..) - (Rier, . $ -t ш,- -,<и a*.w Oiytustjf! , оту ,иьт оч
dx. dx2
d2a
dt2
д3и	з/ л°'и^
-—= -a (x,r)—r +
dt3 V &
, d2a n
~a^^ + a <*
(XUX
d2b , л d2b
dt ' dtdx
Разрешая (52) относительно третьей частной производной по
ря ее значение при х = 0 и t - tn, имеем
+ a2
dt
dx
da 8Z>
dx dx
(52)
и бе-
&b ^da
6?
2
x

где значения найдены в (28), a в (30). Теперь [ы]| из
(48) с учетом (53) полностью известны до четвертого порядка по Л и
» Г 1"
можно определить дополнительные условия как U( = |WJ|, не при-
нимая во внимание слагаемые в (48) более высокого порядка по h,
чем третий, для всех П = 1, N . Для того чтобы найти недостающие
N дополнительных уравнений, повторим аналогичные преобразова-
ния, но для й:
и1=и4<г?^и?|-=+ы:.;у+о(к*)л5<)
Частные производные по х, входящие в (54), ранее уже были вычис-
лены (см. (28), (30), (53)). Соотношения (53) и (54) различаются толь-
ко тем, что в правую часть (53) входит h, а в (54) - 2h. Следователь-
130
но полагая и" ~ с точностью до членов четвертого и более вы-
соких порядков малости по h, получаем недостающие N уравнений,
которые завершают постановку разностной задачи для этого варианта
заданий.
В процессе работы над заданием предполагается выполнение
следующих этапов:
1)	аналитическое решение дифференциальной задачи и по-
строение следа в выбранном конечномерном пространстве;
2)	вывод дополнительных разностных уравнений, где это необ-
ходимо;
3)	исследование разностной схемы на аппроксимацию и спек-
тральную устойчивость, используя принцип замороженных коэффи-
циентов;
4)	разработка алгоритма решения разностной задачи и его реа-
лизация в виде программы, написанной на алгоритмическом языке
высокого уровня;
5)	отладка программы с учетом вычисленных значений следа
аналитического решения;
6)	проведение практических исследований вычислительных
свойств используемой разностной схемы на последовательно удваи-
ваемых сетках:
а)	соответствует ли порядок сходимости численного ре-
шения порядку аппроксимации на решении, установленному в
пункте V.3.;
б)	насколько сильно зависят численные результаты от
выполнения условий устойчивости и что происходит, когда они
нарушаются;
в)	как изменится порядок сходимости, если при задании
дополнительных разностных уравнений оставить меньшее чис-
ло слагаемых в разложении в ряд Тейлора значений следа;
г)	предложить и опытным путем проверить свои вариан-
ты задания дополнительных граничных уравнений;
д)	выяснить, насколько сильно значения решения при
t = 1 зависят от начальных или граничных условий, и объяс-
нить полученный результат.
13Г
СПИСОК ЛИТЕРАТУРЫ
1.	Демченко В.В. Уравнения и системы уравнений с
частными проихводными первого порядка: Учебное пособие. -
2 изд. - М.: МФТИ, 2004. - 116 с.
2.	Магомедов К.М., Холодов А.С. Сеточно-
характеристические численные методы. - М.; Наука, 1988. -
290 с.
3.	Толстых А.И. Компактные разностные схемы и их
применение в задачах аэрогидродинамики. - М.: Наука, 1990 -
230 с.
'Ь’/.Зг (
 
132
ПРИЛОЖЕНИЕ
Вариант № 1
Найти аналитическое и численное решения смешанной задачи
для уравнения переноса в квадрате 0 < X, t <1 согласно номеру за-
дания и сравнить их значения в одиннадцати равноудаленных точках
в момент времени t = 1.
Задание № 1
Дифференциальная задача	•
—+—= cosx; 0<f£l, 0<х£1,
dt дх
м(х,0) = ln(l + x2) + sinx;w((V) = ln(l + J2).
Разностная схема
Dk ={(х,,tn): х, = hl,hL = l,l~О,L; Г =т; тУ = 1,п = 0,у|,
< = м” -«/"_|)+rcosx/,/ = l,Z,« = 0,y-l,
wz° =ln(l + x2)+sinx,,/ = 0,Z; и„ =1пГ1 + (/л)21,л = 1,М
Задание Ns 2
Дифференциальная задача
= sinx; 0</£1, 0^х<1,
dt дх	.
м(х,0) = х2 +cosx; w(l,/) = (l+/)2+cosl.
Разностная схема
D* = {(х, Z) . X/ = hl,hL = 1,1 = (П; /" = пх, xN = 1,п = OjV},
м,”+1 =< +i(<,-uz”)+Tsinx,,/ = 0>Z-l,n = 0^-l,
л
w,° =х/ +cosx,J = 0,Z; «2 =[l+/"]2+cosl,n = l,M
183
Задание №3
Дифференциальная задача
^f_2—=0; 0<Г<1, 0<х<1,
dt дх
u(x,0) = cosx; w(l,r) = cos(l + 2/).
Разностная схема
Dh ={(*,,/").• х, = hl,hL = \,l = 0,L; t" = nx; xN = 1,и = 0,у| ,
«Г = «Г+2~(u,n+l -«;),z=o,£-i,n=ay-i,
h— ___________________________________ •
ui=COSX!,1 = 0,L- w/', = cos(l + 2r”),n = ljy.
Задание №4
Дифференциальная задача
—-(2f+3)— = 0; 0<Z<l, 0<х<1,
St v 'dx
u(0,x) = ln(l + x2y, w(/,l) = ln 1 + (1 + ЗГ + Г2)2
Разностная схема
Рк=((х„г").х, = й/,й1 = и=О;
/° - л ,»+i _ v V1'	________.1
~0,/ -Lv2X=u=o,y-il.
*=о л=о
«Г =ц’ +

ч°

Задание № $
Дифференциальнаязадача
2\21
,и = 1,У.
134
^+(/ + 6)— = 0; 0<f <Z 1, 0<х<1,
dt дх
t2
«(x,0) = x; w(0,/) = —-—6t.
Разностная схема
А = {(х/ ’ *”) ' xi = hl.hL = 1,1 = Q.L;
п	N-1	_____)
t° =0, Г+|	=1,и = 0Л-1к
*=О n=O	J
<+1 = < -^(/в+б)(< -	), i=й,п=оГлПС
и° =x„l = oj; и0" = -o.5(r")2 -6tn,n = ш
Задание № 6
Дифференциальная задача
—+ (х + 4)—= 0; 0 < t < 1, 0<х<1,
dt v ’дх
м(х,0) = (х +4)2; u(0,/) = 16е-2/.
Разностная схема
Л = {(М" )-х/= Ы’М =
Г°=0, /в+| =£т*;£т, =1,и = 0,ЛГ-1
jfc=O п=О
и”+1 = и” ~^(х1+4) («/" - <>) и=й,«=о,у-1,
и,0 =(x/+4)2>Z = oTZ; w0" =16е'2'",л = ш
Задание № 7
Дифференциальная задача
135
<^+4^=х; 0</<1, 0<х<1,
dt дх
д-2
w(x,O) = sinx +—; w(0,/) = -sin4A
8
Разностная схема
= пх; тМ = 1,п
- 4	- <,)+гт,, I = 1, L, п = О, N -1,
О	X,2	--- ------------------
и, =^l+-L,l = Q,L-u”=-s^tn,n = \,N
О	’
Задание №8
Дифференциальная задача
ди ди ,
"^,O) = x + l;M(o,/) = e'_z
Разностная схема
й* =,<V/') :Х'=ЫМ = и = «' =»х; Л = U = о^>
n = 6j\M,
"/°=х/+15/=^7 в f
’ о е -г
иГ' =и‘
Задание Mb 9
Л д“*Ф'Р«ад1льнм,ад,та
и
-t.
дх
1 2
2 ’

= hl,hL = 1,1 = 0,L;
t° = 0,r"+1 = Z	= 1,и=0Л-Д
A=0 kn=0	J
n+l = мл _ MX/ + e‘" VM« _ un^ j +
n y 7	'
+ ^nxi[xl+e‘"'j, l = 1,L,n = G,N-l, ►
M/° = x/+0.5x/2,Z = 0?Z; < = -t",n = IJV.
Задание №10
Дифференциальная задача
ди . ди Л , Л ,
— + е— = t\ 0<t<l, 0<х<1,
dt дх	
м(х,0) = (х-1)2;м(0,/) = е2' +0.5/2.
Разностная схема
А ={(х/,/").х/ = W,A£ = l,Z = 0j;
r°=O,r"+1= Е тр^т„ = 1л = 0,ЛГ-1к
к=О * n=Q
U"' = "" ~	“ U'-') + V"’ 1 = П =
и? =(х,-1)2,1 = 0Л; и”ъ =е21' +0.5(г”)2л = 1Л-
Задание № 11
Дифференциальная задача
~г~+е~х — = х; 0 < t < 1, 0<х<1,
dt дх
u(x>0) = e2x +(x-l)ex;u(0,!) = (t-l)2-1.
Разностная схема
137
А = M,hL ~ ~
? = О, Г'	гя = 1л = OJV-11
4=0	»=0	J
и;+'=и;-^е-д'(«Г-<,)+гя.
а
»>e2x>+(x/-l)eJ»,/ = AZ; г<
Задание № 12
Дифференциальная задача
^._е-^ = х; о<1<1, 0<х<1,
dt дх	Г
«(х,0) = е2х -(х-1)ех; ы(1,1) = (е + 1)2.
Разностная схема
А = {( V") • х, = hl.hL=1,1 = ОЛ- f.=tn. xN = 1, n = ОЛ} ,
«Г‘ = и” +уе-х' («; -М") + TXl,l = 0,1-1,и = ОЛ-1,
м,° =е2* -(ху —1)ех', I = Oj; unL = (е + 1я)2 ,п = W
Задание № 13
Дифференциальная задача
—= е'; 0<f<l, 0<х<1,1	' >
dtdx	У
м(х,0)=х+1; w(l,t) = e’ +1-1. J
Разностная схема
А ={(хР1"): х, = hl,hL = 1,1 = ОХ 1’ = т; xN = 1,« = О,#}
“Г =«; + v(«M -«;) + те1" ,1 = 0,£ - 1л = ол -I	>
u?=x/+l,l=0,L;«2 = e'’+1”-1л=ГХ	
138
Задание № 14
Дифференциальная задача
— -(х + £Г')-^- = х(х4-е"'); 0</£1, 0<х<1,
dt х дх 4	’
х2
н(х,О) = х——; w(l,r) = e' +Г-0.5.
Разностная схема
Dh = {(x,Z): х, = hl,hL = l,7 = oj;
t° =0,Г+| =£ь;£т„ =±1,л = 0,АГ-1 ,
4=0	n=0
»Г,=«Г + т,(х,+е-'')(<-«;)+
ч-ТоХ^х, 4-е-'’), l = f),L-\,n = Q,N•
u* =xl-0.5x2,l = Q,L; u"=e1' +tn -0.5,n = l,N.
Задание №15
Дифференциальная задача
~-(x + 4)—= 0; 0<Z<l, 0^х<1,
dt V ’dx
w(x,0) = (x + 4)2; n(l,z) = 25e2'.
Разностная схема	____
Dk = {(x,j"): x,=hl,hL =	= »x;t:N »1,л = ОЛ) ,
<*• = и” +-(Xl + 4)(<,-и”),I = Q,L~\,n = 0Л-1,
М,0 = (х, + 4)2,1 = о7; UHL = 25е21", п = 1, N.
Ш
Заданием 16
Дифференциальная задача
= 0<^1, 0<Х<1>
dt дх	•
/	\2 t
«(x,O)*(x+lf;«(I-')=(e’+1) +7j	'	'
Разностная схема	___
Д	=М,Л£ = М = 0Л‘
t° = о,г1 * 2л;1л = 1’и=Q’N ~1 [ ’
t=0	»=о
<=«;+^/(<-«;)+^"./=оХЛл=°л-1,
«®=(х/+1)2,/=о1; <=(/ +l)2+0.5(f )2,и = 1Л.	।
Вариант № 2	]
Найти аналитическое и численное решения смешанной задачи т
доя уравнения переноса в квадрате 0 < х, t < 1 согласно номеру зада-
ния и сравнить их значения в одиннадцати равноудаленных точках в J
момент времени t = 1.
Заданием 1	1	I
Дифференциальная задача	'	т
ди ди	Y - * '	|
—- = smx; 0</<1, 0<х<1,	1
dt дх	\	I
м(х,0) = х2+cosx; M(l,/)=(l+r2)+cosl.	'	'	I
Разностная схема	'	1
D*={U><',).x/ = W,/iL = l,/ = O>Z,r =nr;iJV = l,n = OJv} ,	1
Д:
140
uf =u” -^(<2 -4C +х)+^т(<2 -X +«;
rsinx, +0.5Г2,/ = 0,L-2,n = 0,N-l ;uf =xf + cosx,,
1 = 0,L; unL =0 + f" j2+cosl, n = l,N; u"L_x =?,n = l,N.
Задание № 2
Дифференциальная задача
—-2—= 0; 0<Г<1, 0<х<1,
dt дх	•
«(x,0) = cosx; w(l,f) = cos(l + 2f).
Разностная схема
Dh = f(x.,f]:x,=hl,hL = 1,1 = 0.L; t" = m;xN = l,n = 0,a4
2
-К., +зи;)+2р-(с1 -x,+«,”),
/ = 0,£-2,л = 0,У-1;w° = cosx„ I = 0,L,
unL =cos(l + 2f"), n = l,N', unL_x =?, n = \,N.
Задание № 3
Дифференциальная задача
^ + 4 —= 0; 0<Г < 1, 0<х<1,
dt дх
w(x,0) = sinx + 0.125x2;w(0j) = -sin4r
Разностная схема
Da = {(jc,,f ). x,= hl,hL = 1,1 = O. t" = m;xN = l,n = Ojv} ,
Ml
и?' =<+2;(-<г+К,	+«,")+
+txl-2t\l = 2,L)n = 0,N-1 ‘,uf = sinx; +0.125Х,2,
/=о7;^“ =-яп4/\ л=1Л;м!п«???, «=i?M
Задание № 4
Дифференциальная задача  > •
—+(/+6)—= 0; 0<f£l, 0<х<1,
dt к ’ дх
и(х$) = х;и(0,/) = -0.5Г2 -6г.
Разностная схема
" . f =0,Г' -t/blX =1>» = 0Л-1 ,
*=0 п^О
. -
< =«; 4 т^п+6fe(-<2+4«м -х)+
2	) 2л	., i
(f+^\ т2	____ _______
-Л -2U". + и*), / = 2,L,п = О,N-1, 
«°=хр/=^Х=-0.5(Г)г-6г% л = Ш .
«,"=?,л=й7.	-  ; -	-0
Задание № 5	‘	'Л’ *•'  !
Дифференциальная задача
Su / .ч ди	)	>.
—+(х+4)—= 0; 0</<1, 0^х<1,	; 1 '
01 их	1
м (х, 0)=(х + 4)2; и (0, t) = 1 бе"2'.
Разностная схема
Dh = {(x,Z).x, =A/>AZ = 1,/ = o7№ =лт,т# = 1,л = 0Л}
<’ = «Г +(х/ +4)^(1-°-5г)(-ы/-2+К, -зы;)+
+(х, +4)г ^-(<2	+«,”), / = 2Дл = ^П,
И/° = (xt + 4)2,1 = оГЕ;«о = 16e~2f", п = \N,
М]» = ?,и=йл
Задание№6
Дифференциальная задача
—-(2г + 3)—= 0; 0<г<1, 0<х<1,
dt v } дх
м(х,0) = 1п(1 + х2); w(l,f) = In l + (l + 3/+z2) .
Разностная схема
Dh = {(х,, tn) .• х, = Ы, hL = l,l = 0Л;
/o = 0,Z’+)=^т*,^тя=Lл = O>^-lk
*=0	|>=0	j
<	~^п +з+г„)(<2 -4«;+1 +зм;)+^-х
X(2Г + з)2 (<2 -2и"м +u;),l = 0,1-2,П = 0Л-1,
4° =1п(1 + х2), / = ОГЕ; и[ =1п|1+ 1+3/" +('")2]
я=Г^;м;_1 =?, n = UV.	. ;	;
Задание № 7
Дифференциальная задача
- tffe
~+~ = cosx;	0<х<1,
dt дх
w(x,0) = ln(l + x2)+sinx;w(0/) = ln(l + Z2].
Разностная схема
А ={(V")• *, = M,hL = \,l = 67; г" = пт; xN=\,n = Ojv} ,
2
=“’+^(-“«+4“«-х)+^-(<г -2<,+«;•)+
in	in	"
г2	---- -------
+ rcosx, +—smxt, I = 2,L,n = 0,N-1,
u°t = ln(l+x2)+sinxz, I = 0,L;uq =ln
л=1Л;«Г=?>»=Г^ i? 	' <r;
Задание №8	.<	• / '	1 ;;;-1
Дифференциальная задача
—+e'x — = x; 0<f<l, 0<x<l,	. :
dt дх	У
и(х,0) = е2*+(х-1)е*;м(0,г) = (г-1)2 -1.
Разностная схема	s
Dt = |(x, ,t"); x, = hl, hL = 1,1 = 67;'" = nv,’ xW = 1, n = 67^
=u‘+SJe’'(1+7c''V-<I+Kl -3<j+
f2	* /	x ________ i
+ 2^2e ‘ {u"-2	+w")+r^x/ ~2e~X' J’ l~ "
n = 0 л -1 ;U/° = e2x'+ (X/- l)eI(, / = oT",
n=l,N;«|"=?, n=TJv. Л
144
Задание № 9
Дифференциальная задача	,,
^-+— = е1\ 0<Г ±£1, 0<х<1,
dt дх	{
м(х,0) = х + 1; u(O,t) = e‘-t.
Разностная схема
Dh ={(хрГ); х, = hl,hL =1,7 = 0,£; t" =т; тУ = 1,л = 0,^ ,
Т /	\ Т1 I	\
°"”	+4“’' -3“')+^(“« -2“«+"')+
W (1 +£j, / = 2Лп = О,ЛГ-1 ;ut =x,+l,Z=o7;
и” = е*" -tn,n = VN',u” = ?, и = ijv.
Задание №10
Дифференциальная задача -
—+(х + е')— = х(х + е'); 0<Г<1, 0<х<1,
dt V 'дх v f	1
w(x,0) = x + 0.5x2; u(0,t) = -t.	,<
Разностная схема
= {(x, ,t”)x, = hl, hL = 1,7 = Oj;
t° =0, t"+i =^т*,^тп = 1,л = 0,N-1 !•
	»=0	J
145
и/*1 = и; +Л.(/-^х, + х, ](-<2 + 4<г-3«/")+
2h\	2 у
г’(х+/Г	:
+----(и/Я-2~^и1-) +U/ )+	!

2х,(х, + е'
Г 2L '
и=0,Л^—1 ',Uf = X/ +0.5xj\ I — О,Цио
и’=?,п=1Л
+-е2'"
\,l = 2,L, .
Г, « = !,#:
Й. :
•- п ; -
Задание № II
Дифференциальная задача
^ + e'^.s/; 0<Г<], 0<Х<1,,
dt дх	,
г\	US  • •	\	 ‘ •
и(х,0) = (х-1)\м(0,г) = е2'+у.	Ъ
Разностная схема	’Л,.1	- р).)а
А ={(х(,Гв).х/ = /г/,/г£ = 11/ = О£;
Y ' ' '-Л A(V
г° =о,Г' =£vXtb =1,я = 0Лй1,
А»О w=0	-	J .
х(<2-Xi +«n+rBG”	/ = 27Е,я = ОЛ-1,
\ “* /
ч° =U-i)2, /=о7;«0п=г’ +о,5 (г)2,
«Г=?,«=П7.
146
Задание № 12
Дифференциальная задача
=	0<t <1, 0<Х<1,
dt дх ,	(
2 t1
u(x,O) = (x + l) ;w(l,f) = (l + e'j +—.
Разностная схема
D„ = {(х,, t” ): х, = hl, hL = \,l = О;
. , t° =0Г+1	=1,п=0Л-1 >,
Л»0	w«0
-	,	( т \	V Г2
< = и” —1 + -^ ef (i£2 - 4и^ + 3< )+^е21" х
1	1 2h\ 2 J v	’ 2h2
uf = (x, +1)2,1 = 0, Z; unL = (1 + / J* + 0.5(Г )2, n = 1 ,N\
^=?,n = l,N.
Задание №13	.	.	.
Дифференциальная задача
du	Qu
--—е — = х;0<Г<1, 0<х<1,
dt	Qx
w(x,0) = е2х — (x-l)ex;w(l,f) = (e + t)2.
Разностная схема
= {(*/>'")•’х< = hl’hL = *-Z = Г	= OJV} ,
147
<=“ (’4е" )(•“ ~4"'"'+Зк')+
z/t \ z у
-.2	x f t 'X 1 
+^eJ' (<2 -x., +“,")+ П *, *7«ч . I = o.i-2,
Z/>	\	**	/
n=O,N-\; u* = e2x’ -(x,-])/', / = 0Д;
«”=(е+^)\л = й^;«"ч=?, w = ijVC. ч . ~	. /
Задание №14	<	'	'
Дифференциальная задача
= O«Sl, 0<х<1,1	BI - ,
dt дх	У '/
u(x,O)=x+l;w(l,f) = el+/+1. -
: - . - <- ' j ; r : >  ,’1C - ' i ’
Разностная схема
Dh =^x(,r’).-x/ = M,hL~l,l=.0,L;t" = пт; гУ = 1,и = 0,^ у, I i
_	2	1
<=«^^(uu2-K.+3«r)+^r(c2-2<1+<)+
+r^ +^У’, l = o’i-2, n = 0Л-1; W,0 = x, + l,z = oj;>
“!=е,’+^+1,«=1Л;«;_!=?,л=ijv. ' " Л'
'	- '   'Л . ,;! -• ;• •	i i
Задание Ke 15	t и., .v;
Дифференциальная задача , \\,,	„ у p
"^^X+e'^ = X(X+e''); 0<^l, 0<x<l,
x2
и(х,0) = ——+x; u(lj) = e'+f_Q 5
Разностная схема
1W
Dh ={(х/Ди).х/ = hl,hL = l,l = Q,L;
it N-l	___
Z=OZ+'=£т*;£тя=1,л = 0Л-1к
4=0	л=0
=	+ e'' +S2X')("“	+3“>)+
^-(x/ + e~l") (<2 -2u"+, + uttt) +
2h '	'
x, +e ] +—12х,+2х,е +e ) ,
*	/ 2 ' *	* у 3
l = 0,L-2,n = 0,У-1;«/° = -0.5x2 + x„ I = oj;
u” = e‘" +tn -0.5, n = 1Л; u”, = ?, n = ij7.
Lt	7	г г L,—l	'	'
Задание №16
Дифференциальная задача
—--(х + 4)—= 0; Ocf^l, 0^х<1,
dt v ’дх	Ч
w(x,0) = (x + 4)2; w(l,^) = 25e2'.
Разностная схема
Dh' = {(X/, t” ): X/ = hl, hL = 1,1 = Oj; t" = m; zN = 1, n = 0Л
“"+l ~u”	+4)(“/"г -4м/"| +3u”) + —2-x
2h\ 2 J v	2h
x(x, + 4)2 (u;+2-2u?+l +u”),l = 0^2,n = 0,ЛГ-1,
““ = (xt + 4)2,1 = OTT; uttL = 25e2'', n = ijV; <, = ?,
n — \,N.	 •	.
149
Вариант №3	»	/
Найти аналитическое и численное решения смешанной задачи
для уравнения переноса в квадрате 0 < х, / < 1 согласно номеру зада-
ния и сравнить их значения в одиннадцати равноудаленных точках в
момент времени t - 1.	•
Задание Ns 1
Дифференциальная задача	 j j ‘	. i ......
—-2— = 0; 0<t<l, 0<х<1,
St дХ	'	‘ '  	. j ;	.
a(x,0) = cosx; w(l5r) = cos(l + 2r).
Разностная схема	••	\	\
^={(x/^")-^=^A£ = l,Z = 0,Z;f=«T, Ty = l,w = 0j^] ,
+^к>+18<	+
2т2 /
+YVmw+4<2-5<> +2«г)+
4т2/ „	_____ ’________
+ 3й2 ^М?+3 “^/+2 +	" й"	= 0, Z - 3, п = 0, jV-1, 
и‘ =c°sx/5 I = QfL;u[ = со50 + 2г"), я±=1Д7; * \ -
С|=?5п=1Л;м2_2=?, л =	) ч
= .j  Ч1 и; ", ; - i Л :
V ь f	* Л	j
Заданием?	'	'* f ‘ и" ;5
Дифференциальная задача
4 :и <- « М (J- ./)•<
dt Qx ~х; 0<*Sl, 0<х£1,
"(M)=sinx+0.125x’;a(0,,)_sin4,.| '
150
Разностная схема
Dh = |(xz, tn)  x, = hl, hL = 1,1 = 0, L; f = пт; xN=1, n=0, N
“Г = “'+7г(2“’->-9“'-=+18“'-'-Н“”)+ ‘
3hv	7
• +^.(-<,+4И-г-5С,+2И;)- .
32т2/	\
~ +Зи« "3<-+B’ +
3h '	'
' 1	+ТХ/-2r2,1 = 3,L,n = Q,N
=sinxz +0.125xz2, / = 0,L;Ug =-sin(4/”), n = X,N;
u” =?,n = bN;u” =?>n = ijV.
Задание № 3
Дифференциальная задача
— + (г + 6)—= 0;	0<х<1,
dt v дх
u(x,0) = x; u(0,t) = -6t-0.5t2.
Разностная схема
Dh = {(x/Z ) .• x, = hl.hL = ],Z = 0Д;
fo=0,^,=2vgTjr=l,« = 0JV^L
4=0	л=0	J
151
-Ц«Г)+
* «М 2 J
+ jLff +б)(гв +t* + б)(~«М + 4и/-2	+ 2W/ )"
2Й ________________________________ . ?•
+ бУ(-<з+И%_3мГ-1+м")’^ = 3’^’	>
6h3 V 7 V ________	2	____
п=ОЛЧ;«,° =	=°’£’ “о =-6/" ^5(7") ’п =1 л’
«;=?, «=ш«;=?,«=1Л- f- i
Задание 4
Дифференциальная задача
—+ (х + 4)—= 0; 0<Г£1, 0<х<1,	. . _>q - J; Г
dt v 7 дх	I
w(x,0) = (x+4)2; м(0д) = 16е"?'.
Разностная схема	.,
Dh =(h-f")• = w.w.=1/=OZ;f = nx;xN = l,n=О), ,
7„
x(2<3-9<2+i8u/n_I-iiu;)+
+W^X,+ 4)211" rX-<3+ K-2 - 5«м + 2u”)	X .
XU+4)3 (-««+3C2 -3C. + M;), I ~3j,n = Qjri;
“'° =U+4)\/ = 0Л; «0й =16^, n = ijv,
<=?,n=i,iV;^=?,n=ijv.
133
Задание № 5
Дифференциальная задача *	‘
^•-(2f + 3)—= 0; 0<f<l, 0<Jx<l,
dt дх	> I
ы(х,0) = 1п(1 + х2); м(1,г) = 1п 1-ф+Зг+/2)2 .
Разностная схема
ДА =^(х,х; = hl,hL = l,l = 0,L; Г = пх; xN = 1,и = 0,л|
=и”+^(2z‘+3+г)(2в“ +18< -1 ’“-)+
+^(2'" +3)(2Г+3+2г)(~<> +4“« -*й, +2«;)+
All
Л и;,-з<2+к,/=о,£-3,
0/7	'
п; и° = 1п(1+х/), 1=оХ ’
unL = 1п{1 + 1 + f +
,»=1,М
U”L_, = ?, п = 1 Л; unL_2 =%n=l,N.
Задание № б
Дифференциальная задача
ди ди	.
"г~ + —- = cosx; 0</<1, 0<х^1,
dt дх
w(x,0) = ln(l+x2)+sinx; и(0,/) = 1п(1+/2)?
Разностная схема	___
А «{(х/У^х, =hl,hL = l,l = O^L,in=m,-xN = ln = 0,Ni


м->=1<+Х(2<3-К2+18<1-'И<)+ .	: .
6Й
+7Г7г(-Ч-з + ^w/-2 ~ $ui-i +	) ~
2h v
з	-	" '. 7' '; : • ;)	7
-77т("М« + ЗМ’2 "3м'-> + U> ) +
6й v	?/ .
,	Т3	---••	----~
+г cos X) + 0.5г sin х(---cos х(, I - 3, L, п = 0, N -1,
6
wf = ln(l + x2)+sinx/, 1 = 0,L;Uq
< =?, n^lN-u” =?, n-UL
1
и = 1,К;
ЗаданиеM7	. ;	.... 1 .7'-'
Дифференциальная задача	” • ‘
—= sinx; 0<t<1, 0<x<l,	; :	;
dt dx
u(x,0)=x2+cosx; «(l,t)=(l+f)2+cosl, ”.	' ; , f !
Разностная схема
Dh = ((*//)  x, = hl,hL =\,l = ^L;Ojv} , '
<,=4”^(x.1-k1+ik.-h<)+^x . 7.
х('мм+4«”+2~5ы"+1+2«/"^ + rsinx(+0.5r:!cosx?-+‘
+^5-(<3-3Cj + ЗмГ+1-М;)-~sinxf,7=0^3, - 7
n = 0,N -1, u° = x2 + cos x( ,1 ~ О, L; ueL = (1 +1 ” )2 +.
+cosi,«=1JV; =?, w=ijV; «2.2=?, «=мГ.

15»
Задание № 8
Дифференциальная задача
^+2—= 0; 0<t <1, 0<х<1,
dt ох
w(x,0) = cosx; w(0,f) = cos2z.
Разностная схема
Dh	Xj = hl,hL = \,l = 0,L; t" =т; zN = 1,п = 0,N
< =<+^(2<. -xi+ix -11»;)+
2r2 /	x
+~^2~\и/-з + ^u,~2 ~5ui-t + ^ui ) —
 Iff	+3“’1  3"«+W>” - 57=1
u° = cosx/} I = 0,£; Ug = cos2f, n = l,N,
w,” =?,H=ijV;w2"=?,H = V^
Задание №9
Дифференциальная задача
^•-4—= x; 0<f<l, 0<х<1,
dt dx
«(x, 0) = sin x - 0.125x2; и (1, t) = sin (1 + 4t) - 0.125.
Разностная схема	___
={(x,Z).xz =W,AI = lJ = aZ;/’=«T;'^ = I’n = 0’jV}
И5
<=<+$(KJ-’“M + I8“«-1I<)+' ' >
' 3лv
or2	 .	• 4	i
+^-(-4%+‘t«M-5"«+2“i")+
Л v
+ 'Т7Г‘(м/+з ~^um + ^um ~ui )+rx!	»
z=o,z-3,»=o,y-i,
uf = sinx(-0.125x2, l~G,L; u”L = sin^l + 4r"J-0,125,
»Ww=?,n = U;«; 2 = ?,и = 1Л?	'tV
Задание № 10
Дифференциальная задача
ди ди 1 п п 1 1 ,И 4 ‘	' 1
— + — = е; 0<t^l, 0<х<1,
St дх	У	г<1, •
w(x,O) = x+l; w(O,t) = e'„ ,
Разностная схема
А ={(*/.<")••*/ =W,hL=ll = ^L;f =nx;xN = 1,п = 0^N.
,n J
< = <+£ (2<, - 9u;_,+1 &c, -1К
on
+w ^u”~3+4m"-2 ~ 5m”"’ + 2m" } ~
+X-г -3u;_, + u”)+re1
1+-+ —
2 6 J
l^3,L,n^Q,N-l
Задание Ns 11	.
Дифференциальная задача
^4x+e/)zT = X(X + e'); 0<^L 0<х^1,
dt v St v 1	, 
ф,0) = х+0.5х2; w(0, =	?_
Разностная схема
A=hl,hL = \,l = 0,L;
= £v£T„=U = 0,tf-l >’
,	i«0 Ы)
<‘=«r+ .	....
(2u“ ~9u"2+ ISu" -Hu"
- Tn*i )(-w/-3 +4«/"2 -5<! +2«;)
c +/V
~75	(_w/-3 +3«/_2	+u/)'
6n '
')_£L^+2x,e'-+2x,i)+^x
12^	4 6
x(4x, + 3/ )p = 3,L,n = 0,N~ 1, uf =x,+ 0.5x2,
l = Oj-, u”=-t\ n^N', u;=?,n = UN; u" = ?,n = l,N.
Задание M 12
Дифференциальная задача
^—(x + 4)~ = 0; 0<t<l, 0<x<l,
dt ' !dx
"(x,0)={x+4)2;«(l,/) = 25e2‘.
Разностная схема	___ ____________________
Dh = {(х,Г) •• х, = М, hL = U = 0,Z; f1 =пт; №*1,п = 0,
м;,+*=м;+
г ( т т1
6h{ 2 6
+^т(1 + ^)(х/+4)2(~им + 4им "5м"+1 +2м/")
2п
+
(х, +4)(2<3 -9<  + 18«”, -!!«;) +
им Зим + 3w,+1
Z = 0,Z-3,n = 0,W-L
и’ =(х,+4)2, / = 0Д; и” =25е2'\ и = Ш . ? < ~
и;_1=?,п=1Л;<2 = ?,» = 1Л.
Задание Ns 13
Дифференциальная задача
ди , ди
--e~ = t- 0<(<l,0Sx<l,
и
2

Разностная схема
f - О /Я+1 _ V г=}	; . .	т--_
°'	_2^£т„=1,и=о,лг.-1
*“в п=О
•4мг
i

158
s»—«*r
+77,1+V+7’К (2<= ~9"« +l8“« 4!“”)+
Ort z о у	' I
+^(1 + ^)e2'" (~UM +4m/% ~5uM +2<)+
zrt	'
+ -Ц-е3'’ (w/.j -3«,"2 +3«," -и,а)+т(t" +—1 >
6Л3	\ м м м !)	2 J’
/ = 0,L-3,n = 0,N-l,
M/° =(x,+1)\ Z = 6^L; m” =(l+/)2+0.5(f)2, и = 1Л;
unL_} = ?, n=ijv; =Ъ n=Гж
Задание №14
Дифференциальная задача
~ + 8~ = Г; 0<t<. 1, 0<х^1,
St дх	•
w(x,0) = ех; «(0,r) = е-*+0.5/2.
Разностная схема
Dh = {(х,Z):х, =hl,hL = 1,1 =	tH ^m;rN^l,n = О,ЛГ|,
159
=<+^(2«Д]-9«,:г+18<-И<) +
Зй	,	'
+V’f'”'-j+4“'-i'5“':'+2“"^ *	£
-^-(-<1 +3»’, -3»,:, + <)+«• + 0.5r2,'	/
'	-• > /«ЗД,» = 0,#-1,
U,° = e\ I = OJ>, u” = e^1" + 0.5('")2 > « = W
<=?,« = 1Л,и2л=?,и = Ш	.r ,
•	• «•  t -v ® I , ( _ -t- у !
Задание №15	' - v Л'- f ‘	' n
Дифференциальная задача
=	0<Г<1, 0<х<1,
dt дх	>
ы(х,0) = х+1; uCl,t] = e' +/+1. ,	\ ' Ц '< п
Разностная схема	v. ц--,
х, =hl,hL=\,l	t" =n<т^Г = 1’n = О?Л} l'
< - <	+18< -11<) +'. ’
t2 t ‘	; - v- r3	.-л'1 •
+ 2/? \~Um + 4"i+2 ~ ^um + 2«") +	x
,	v (	f2
x(«;+3 -x+I+з<, -«;)+rii+o.5r+—j
l ~QyL-\n = O,N-1;«? = XZ +1,1 =~^L‘,
unL=/ +гп+1,и=1Л;«2_1=?,п=Ш
<_2=?,n=ijv.

160
Заданием 16	< I / ч
Дифференциальная задача
^“-6—= 0, 0<7<1, 0<х<1, И f , /
dt дх .	. - . . -	> •_< ,
м(х,0) = 1п(1 + х); м(М) = 1п(2 + 6г).
Разностная схема
Da =^(х/,/").х/ = Л/,Л£«=М = 0,1,7" = rtx;xN = 1,« = 0,2vJ ,
<’ =»;+^(2»;., -9<; +18<, -11»;)+
+£5-(-<,+4»;.г-к1+2<)+
п
~3”“ +3”'*'
- и° =ln(l + xz), Z = O,Z; = 1п(2+6Г), »=1,У,
Wi"_1=?,/J=ijv>M;_2=?j7J=ijv;	J
; •	2	. 	'М ''	' ’	’-'Г
'  '
• ’ * '	(- г	,у -	‘	‘
~i'. V - : > -	м ,1 '!. 	- 	4	’		:	J_
.	^.ц,! I-.	{,,: . 	. 	'	
’(«г	.	»	: :1. 	.  •	•<	'	‘	•
_ .	.ч		’	' '
лч .	. . ,	.	L	*	'
- хчч> !Л. ..	•••  :	'
'Л>-	'ч '.ч'.г-., нц;,	' 	 "	! -;
UWluiK,,,-I ..	-	' I 1	’
161
РАЗДЕЛ VI
УРАВНЕНИЯ С ЧАСТНЫМИ
ПРОИЗВОДНЫМИ ПАРАБОЛИЧЕСКОГО
ТИПА
Лабораторная работа № 6	'	.<
Методы решения квазилинейного уравнения
теплопроводности
Введение
Процессы диффузии вещества и диссипации энергии постоянно
происходят в окружающем человека материальном мире и
свойственны всем известным средам: газовой, жидкой, твердому телу,
плазме. При определённых условиях их вклад оказывается основным
для описания физического состояния того или иного материального
объекта, например, в случае диффузии нейтронов при ядерных
реакциях или распространения волны электронной теплопроводности
при воздействии мощного лазерного излучения на вещество. Для
точного количественного анализа необходимо использовать строго
обоснованные математические методы, которые в данном случае
связаны с решением уравнения с частными производными второго
порядка параболического типа, называемого уравнением
теплопроаодностн или диффузии в зависимости от физической
природы процесса. Сложности, возникающие при математическом
решении задачи, обусловлены рядом факторов:
во-первых, нестационарностью и пространственной
многомерностью происходящих явлений;
во-вторых, актуальные для современной
проблемы порождаются нелинейными явлениями
науки и техники
—......и процессами, что
требует новых математических подходов для их изучения;
в-третьих, ускорение научно-технического
ис™е™«яТТСЯ С0КРащением сроков выполнения
исследовательских и опытно-конструкторских работ.
методов0 ппи^ ВеДё^ К б°Лее ШИР°КОМУ применению
методов прикладной математики сен------
современных быстродействующих ЭВМ.
прогресса
научно-
численных
основанных на использовании
162
VI	.1. Квазилинейная пространственно одномерная
задача
Из классической теории плазмы известно, что коэффициент
электронной теплопроводности пропорционален электронной
температуре в некоторой положительной степени. Задачи
распространения тепла или диффузии вещества, в которых
коэффициент теплопроводности или диффузии зависит от значений
искомой функции, будем называть квазилинейными. Вначале
рассмотрим пространственно одномерный случай, когда по двум из
трёх ортогональных пространственных направлений искомая функция
оказывается однородной. Поставим для него в единичной области
(О « t « 1, 0 « г « 1) смешанную задачу с краевыми условиями
первого рода в одной из трёх систем координат: декартовой,
цилиндрической или сферической:
ди I д . v „ди.
— =-------trvu^—)
dt rv dr	dr
u(0,r) = <p(r),	0<r<l,
u(t,O) = y/Q(t),	0<t£l,
u(t, I) = </)(/>,
где v = 0 для декартовой, v = 1 для цилиндрической, v - 2 для
сферической систем координат, р.>0 -const., <p(r), i/r0(t), ^,(/)
непрерывные функции своих аргументов, обладающие требуемой
степенью гладкости в области интегрирования.
Построить точное решение задачи (1) для произвола
функций </>(r), ys0(t), ^(0 не удаётся, но для некоторых частных
случаев можно найти решение в аналитическом виде, что
чрезвычайно полезным для выбора модельной зада пр
вычислительных программ.
Р)
i А '
VL1	.1. Частные аналитические решения
Будем искать решение уравнения из задачи (1)
« = Т(ОЛ(г), следуя методу разделения переменных,
подстановки этого выражения в уравнение получим
dt rv ar *
ГбЗ
(2)
(3)
Деля правую и левую части этого выражения на
и = T^'^Rfr), имеем
1	dr-Q = —!——(/•''Я'(г)—) =- к
T^'(t) dt R(r)rv dr dr
где к - некоторая константа.
Рассмотрим вначале первое из уравнений (2)
—+ кГ*' =0
dt
и найдём его решение, так как переменные разделяются:
Т --------------------!_____.
(СО+ДА//
Второе уравнение из (2) запишется в виде -'
d . rv dR?"	:
—(----r~~)+^ = 0.
dr p + 1 dr
: (4)
(5)
Будем искать частное решение (5) в виде R — r и сделаем эту
замену
+ + + (6)
Потребуем, чтобы показатели степени у двух слагаемых в (6)
были равны
m(p + l) + v-2~v + m, .
и разрешим получившееся равенство относительно т:	(	.,.
т = 2/р.	,	.	,	(7)
Учитывая (7), найдём нз (6) выражение для fa .	'
-	•	(8)
д*
В результате частное решение уравнения (8) Из задачи (1)
можно записать в виде
__________г7*
// « В (9)
р	-	•	'.-.Л
где Со ~ произвольная постоянная. !	ч.-, - '
Заметим, что в декартовой системе координат при v = 0, если
пространственную координату обозначить х, частное решение (9)
представимо в бп™₽ «*•"-
представимо в более общем виде:
(С.+хЛ
/	V '	'••• ••••’
1деС( - произвольная постоянная.
Зная частное решение (9), сформулируем окончательно задачу
(1), решением которой оно является
">	м(0,г) = Л?‘А 05r5i,
и(/,0) = 0, Oct <1,
f	/]'А
И
, или, повторяя те же рассуждения для (10), имеем
«(0,х) = (С, + х^С'/^ 0 5x51,
и(/,0) = C%[CQ - 2^tjl/f А 0 <t < 1,
0 = (С, +1 А(С0 -
*
Если постоянные Со, входящие в решения (9),
положительны, то за конечное время
.. * /
(")
(12)
(Ш
• С°Р
2(^/ + Х^ + 2)
Функция и обращается в бесконечность, 2^ иодирования,
постепенным ростом градиента функции в	„тгугствует и
Пр» С, ..	»„« нер» ««»»»
теплом, волна не выходит за эту грань, хота промсх
163
градиента а следовательно, и теплового потока внутри области
интегрирования вплоть до бесконечности. Полученные частные
аналитические решения описывают эффект так называемых
“тепловых кристаллов”, когда за конечное время происходит
накопление энергии в фиксированной области пространства без
распространения тепловой волны за её границу.
Данное свойство решения имеет важное значение для
тестирования вычислительных алгоритмов, так как позволяет
проверить сходимость численных решений к следу аналитического
решения для всё более увеличивающихся тепловых потоков и
градиентов температуры.
VI.1.2. Численное решение задачи
Для численного решения задач, аналогичных (1), обычно
используют либо явную, либо неявную четырёхточечные разностные
схемы, которые, имея одинаковый порядок аппроксимации по каждой
из независимых переменных, существенно различаются по
устойчивости. Для того чтобы понять, насколько это существенно для
рассматриваемых задач, начнём изучение разностных схем с
линейной задачи, записанной в декартовой системе координат,
полагая при этом коэффициент теплопроводности постоянным и не
зависящим от искомой функции и:
ди гдги
— ~а ------,
dt дх2 .
м(0,х) = (г>(х), 0<х<1,.
“('>0) = ^о(О, 0<t£ I,
«(М) = ^(0.
Введём в области интегрирования сетку
: х> =*,1 = 0,L,hL = l;f = иг>и = (
запишем сначала явную разностную схему;
= a2	« = 0JV-L
Г	L2	>---------.
г о ___ А ( = 1,1-1,
а затем-неявную;	1	’
- : I.
(13)
•И
(14)



166

(15)
<-u” я=ол~1
T	h2 ’ l = \,L-\,
ч°=^,/ = 0,1;МоЯ=Го%«2^^Я’и = 1’^
Если исследовать разностные схемы (14),	(15) на
аппроксимацию, то окажется, что в обоих случаях для нормы невязки
справедлива оценка
||5/а>||<стт+са/Л
где CT,Ch - некоторые постоянные, не зависящие ни от т, ни от Л.
f Для исследования устойчивости основных уравнений
разностной задачи (14) спектральным признаком (см. [1]) необходимо
в них сделать подстановку ut = Х. е и потребовать, чтобы
|Хл+1/Х"|<1 для любого ае(-со,+со). В результате после
сокращения на e,aJ и обозначения X = Xn+I / К приходим к
выражению
. • 4та2 . 2 а •
Х = 1——sin2—.	(16)
h 2
Спектральный признак устойчивости будет выполнен, если X из (16)
4та2 . ,а ,
удовлетворяет двум неравенствам —1<1----------—sin —<1 для
h 2
Vaef-oo.+oo^-
Правое неравенство выполнено при всех a, а из левого неравенства
получаем, что
т<^у.	(17)
2а
Это означает, что при построении разностной сетки расстояния между
соседними узлами по независимым переменным г и х должны
удовлетворять неравенству (17), если в этом конечномерном
пространстве решается задача (14). В противном случае при
нарушении неравенства (17) сходимость численных решений к следу
аналитического решения будет отсутствовать.
Проверим спектральную устойчивость для разностной схемы
(15). После аналогичных преобразований, выполненных ранее для
явной схемы (14), имеем
16?
1
Требуя выполнения епеирдакого признак, уиойинностя дня (18).
получаем да яераьевстм -1S	7а^' КОТОРМ'
,+~7Гял 2
выполняются при Vaef-oo.-K»?-
Эго означает, что устойчивость будет существовать при любом
соотношении между т и h, или, другими словами, разностная схема
(15) безусловно устойчива.
Возникает вопрос: насколько сильно ограничение (17) влияет
иа проведение вычислений по схеме (14) по сравнению с (15)? Для
ответа иа него определим число узлов разностной сетки, задаваясь
некоторым h и считая, что для схемы (15) x~h, а в случае (14) -
x~h2. Тогда для явной схемы общее число узловых точек 5 ~Z3, а
для неявной - S -L2. Другим важным фактором, характеризующим
свойства однородного численного алгоритма решения задачи,
является среднее число арифметических действий, затрачиваемых на
каждый узел сетки при нахождении неизвестного значения сеточной
функции иа новом временном слое. Для обоих рассматриваемых
случаев оио приблизительно равно десяти, если неявную разностную
задачу (15) решать экономичным методом прогонки.
Обозначим через К}А среднее число арифметических
действий с плавающей запятой, которое необходимо выполнить для
нахождения неизвестного значения сеточной функции на новом
временном слое в одном узле сетки, если использовать разностную
схему (14). Аналогично, К.5 - среднее число арифметических
действий для схемы (15). Разберём два случая: в первом общее число
узлов по простраистаениой переменной
L = 100, а во втором - L = 106, Тогда Ql4(l00) * К{4 * 10й - общее
число арифметических операций с плавающей запятой, выполненных
при решении задачи с L = 100 по схеме (14), a Q5(100) » tf15 * 104 -
по алгоритму (15). Если предположить, что при вычислениях
используется современная ЭВМ с быстродействием 1012 операций в
Ш
секунду, т0 разница в реальном времени будет
Ситуация заметно меняется при
незначительной.
jr = Ю6. Для разностной схемы (14) ^(Ю6)»^ *10'», а
варианта (15) - Ql5Q0 )«Kl5*]021 т.е. в0 вариан1е
арифметический процессор будет загружен несколько десятков
секунд, а в первом - несколько десятков суток. Из приведённого
анализа видно, что при использовании алгоритма (15) существенно
сокращаются объёмы вычислительных работ при расчётах иа мелких
сетках. Учитывая результаты анализа устойчивости разностных схем
(14) и (15), в дальнейшем остановимся более подробно на аналоге
разностной схемы (15) при решении задач (11) и (12). Сложность по
сравнению с линейной задачей (13) заключается в нелинейном
коэффициенте теплопроводности. Если его значения вычислять по
значениям сеточной функции иа новом временном слое п + 1, то тогда
алгоритм прогонки уже не применим для решения разностной задачи
в силу нелинейности системы получившихся уравнений. Выбор
других алгоритмов ведёт к значительному увеличению времени
расчётов. Выход из этого положения состоит в использовании
итераций по нелинейному коэффициенту теплопроводности и неявной
аппроксимации градиентов искомой сеточной функции. Выпишем
разностную схему в конечномерном пространстве:/)^ )•
Xj =/Л,? = О,£,ЛЛ = 1;/л =w,M = 0,y,ry = [C0///2(^v+/z + 2)] Д,
Д > 0 - const} для задачи (II):
-U1 =	_(в;.,
т	2г/ h2

2г* h
n~Q,N-1, J~ min£:max
U, ~UI
4+1
И/
<£,£-COnSt >X-0,
~V ..o - l =
+	'	°	Ц9)
Для решения системы уравнений (19)	прогонки, алгоритм
слоя иа п + 1 можно воспользоваться
ДО
применения которого был ранее подробно разобран в [2], иди
ознакомиться с соответствующим материалом по дополнительной
литературе [3-5]. Известно, что алгоритм прогонки принадлежит
группе экономичных численных методов решения систем линейных
уравнений, причём среднее число арифметических действии,
затрачиваемых на нахождение неизвестного значения сеточной
функции в узле сетки, приблизительно равно 10. Алгоритм прогонки
вычислительно устойчив, если матрица системы линейных уравнений
разностной задачи удовлетворяет свойству диагонального
преобладания. Следовательно, основное отличие в применении
метода прогонки для решения задач (15) и (19) состоит в том, что при
рассмотрении квазилинейной задачи для перехода с и временного
слоя на п + 1 требуется многократная прогонка, т.е. использование
этого метода J раз, пока не будет выполнено условие
|u,t+1 —и,*	------
-j~- <£, /=1,£-1.
расчётов, итерации быстро сходятся и двух-трёх
вычислении, как правило, оказывается достаточно для
заданной точности е.
шах
Как показывает опыт
численных
повторных
получения
VI.2. Квазилинейная пространственно многомерная •'*
задача
Хотя пространственно одномерные задачи теплопроводности
и представляют определённый интерес для исследователей как
наиболее простые, на которых выясняются основные свойства
решаемых проблем, но всё же более часто при проведении
практических работ приходится иметь дело с пространственно
многомерными случаями в декартовой, цилиндрической и
сферической системах координат. Остановимся подробнее на каждом
из перечисленных вариантов и рассмотрим несколько
fnWrnaUATBait-iA *— ---
Vi «мерными случаями в декартовой, цилиндри
сферической системах координат. Остановимся подробнее
из перечисленных вариантов и рассмотрим
пространственно двумерных задач.
VI.2.1. Декартовая система координат
Запишем смешанную задачу для квазилинейного
теплопроводности с краевыми условиями первого рода
D={0,х,у): 0< i,x,yS 1} декартовой системы координат;
уравнения
в области
ед
ди	д	р ди	д	ди
<.	—=—(« —) + —(ir—),
dt дх дх	ду	ду
и(0, х, у) = tp(x, у), 0< х, у < 1, 0 < t < 1,	.	(20)
u(t, 0, у) = (/, у), u(t, 1, у) = (/, у), 0 < у < 1,
u(t,x,O) = %0(t,x), u(t,x,]) = X,(t,x), 0 < x < 1.
Получить общее аналитическое решение задачи (20) с
произвольными начальными и граничными условиями пока не
удаётся, но и частные решения играют важную роль, т.к. на них
можно производить тестирование алгоритмов и отладку программ, не
говоря уже о том интересе, который они могут представлять для
физического описания процесса. Разберём один из возможных
способов построения частных аналитических решений для
многомерных задач.
Построение частных аналитических решений
г--- 
Следуя методу разделения переменных, подробно
разобранному для пространственно одномерной задачи, представим
решение в виде
и = (Сг + х + у^[С, -,	(21)
где Сг, С, — произвольные постоянные. Непосредственной
подстановкой выражения (21) в уравнение с частными производными
из задачи (20) убеждаемся в том, что оно обращается в тождество.
Тогда задача (20), для которой функция и из (21) является решением,
записывается следующим образом:
dt дх дх ду ду
0£х,у<1, 0<Т<1, u(t,0,y) = (C2 +уУ/р[С, -
О £ у <> 1, ы(/,1,у) = (С. +1 + J)^[CZ
OS y< i,u(t,x,0) = (C2+x)^[Cl-^^]^,0<x<l,
(22)
H%rr 4^ + 2)fV
w(r,x,l) = (C, +X+1) (G_ u J
r*
. (22’)
ff?'.
Алгоритм численного решения
Определим разностную задачу на конечномерном пространстве:
Dk ={(х„уя,(йУ.х, = hxl,l = 0Л hx £=1; ут = hym, т = 0,М,
ЪуМ = l;f = пт, п ~ fyN, tN = С,/л	+ 2)- Д, Д - const} ,’
введённом в области интегрирования. При её построении учтём
результаты анализа устойчивости, выполненного в пункте VI. 1.2.
В результате получим неявную разностную схему, аналогичную (19).
Хотя эта система уравнений и будет линейной относительно искомых
значений сеточной функции на к + 1 итерации с учётом того, что
коэффициент теплопроводности аппроксимируется по величинам с
предыдущей итерации, непосредственно использовать метод
прогонки для её решения невозможно, так как матрица системы
уравнений имеет не трёхдиагональный вид.
.к<.гжгхе...
(Ч.ш
Т	2hX
2^	2Л„2	-
?2^
---- ___
С 4
= шшАутах
<е, г-const к и!
^Г~

’ uQ
’ “/.I
<-«-1, т^М-1, „ =ojvTii =(Ci+Xj + У1т/.с%
(23)
1*2
$ =(СЛ^)^[С-4(^ + 2)Г'/^0,/ = 1,1-1, 1	;
/	, -у _________ (23’)
<Сг+ х' +1)/“[С/ -4(д + 2)Г'/А]^,Я = 0Л-1.]
Для преодоления этого затруднения при решении (23), (23’)
используют метод расщепления по пространственным направлениям,
когда в пределах одного временного шага т сначала дают теплу
распространяться только в направлении х, вводя как бы по координате
у теплоизолированные границы, а затем - наоборот, тепловая волна
движется в направлении у, а границы по х - теплоизолированы. В
результате на каждом из этих полушагов для решения
соответствующей подсистемы разностных уравнений можно
применять метод прогонки. На первом полушаге при фиксированных
значениях п и т
(т = 1 -ь М - I) решают методом прогонки подсистему линейных
уравнений относительно значений й,т, где I = 1,Z-1, и учитывают
граничные условия при I = 0 и I = L:
М. - и1т [(“/+> т У + (и/т )* ](«/+1 „ ~ Ц )
г,лг tjti <•'	4 /,/и' jv /+1,ж
г~~-	2Л;
\т *(С, +Уж)^[С;-4(д + 2)Г-*7^, т = 1,Л/-1,
= (С, +1 + /'[С, -4(// + 2)Г ’ ///]'^> п = 0,N-\.
Эта подсистема уравнений является частью более полной задачи (25),
(25’):
<-.• : Г •	2ЛЛ2	' .
(25)
№
к<г+(<,гкС-С
2A;

< E, s- const >
(25’)
= ы,
(26)
k=n,J, J = jminimax
Hn*’ = ц7*',/ = 1,1-1, /и = 1,Л/-1, n = 0,N~\,
и1я^Сг+х,+уя/,‘^,1=И,пг = 0^,	_____
u£ =(C, ★y^lC. -4(P + 2)f*' 1^, m = 0,M,
u% <Cf 4.1+y^lC, -4<^ + 2)Г’ lfx]^,n = Q,N-l,
=(Cx+x,)^[C,-4(д + 2)Г7/д/\ /=1,Z-1,
<« =(C,+x/ + l/"[C(-4(д + 2)/я+7д]^, n = 0,tf-l.
Если в (24) ввести обозначения ______________
с:
Ь, =1-ц -с„«0>ж 4 (Сг + У^[С, -4(// + 2У"7>^]
«д_. =(С, +1+у./*[С, -4(я +2)Г’ /^]
то соответствующие уравнения перепишутся в виде
^О.ж ~ Ы0,ж ’	' '
а,“/+|.и + ^«(.м +	= df,l = 1,L -1, ;
Далее подсистема (26) решается методом прогонки для всех
т ~	1	-s-	М-\. На втором полушаге значения
u/.m' / ~ Ы -1> т = \,М -1 считаются известными при тех значениях
и и Л, что и на пераом полушаге, и последовательно для значений
/ 1 Z- - 1 решают вторую подсистему уравнений относительно
►
174'
.[«,г4<гкс-О
г	2h;
2й.2

(27)
С =(С, +»,/'[€, -4(А + 2)Г'
Ц.М =	+ Х1
Обозначим
а„ = 4(«*„+1>
д,и = 0Л-1.
=4(«,\ Г+(«*»-. П
Ът = 1" ат ~с„,и10 = (С, + х^[С, -4(jj + 2)Г'/р]%,
И
(28)
«/X = (Сг +х< + J/4C, -4(// + 2)Г' ///]/", dm =uljK
перепишем подсистему уравнений (27):
Ч.О — Ч.О’
+ЬЛш =dm, m =	.
Ul^f =И1,М>
где I = 1 -г L - 1. После решения (28) методом прогонки определяем
все M*m4 = l,Z-l, т = 1,М-\ и сравниваем со значениями w*m с
предыдущей итерации:
Ч.ж -Ц
W,**1
/.да
где а - требуемая точность. Если неравенство (29) невыполнено хотя
бы для одного I и т, то к увеличивают на единицу и выполняют новую
итерацию, состоящую из двух полушагов, и так до тех пор,
будет удовлетворено неравенство (29). Последнее к,
и**' -ик
и! т	т	.	__________
й**7
••/.да
определения значений искомой сеточной функции на
временном слое и + 1, полагая и"*1 =«/*',/ = 1,^-1, от = 1>^_1-
Этим завершается переход с и временного слоя на и + 1 для всех
J - mini: max
<£,£- const , используют
(29)
пока не
равное
ДЛЯ
новом
. max
175
„ = о — jV - l. Заметим, что, как правило, требуется не больше двух-
трёх итеряпий по к для. перехода с одного временного слоя на
следующий.	':
'  Ч ,	4,( И
VI.2.2. Нилинярнческая система координат	---
Частные точные решения задач
Приведём постановки смешанных задач для квазилинейного
уравнения теплопроводности с краевыми условиями первого рода в
области	цилиндрической	системы,
координат/) = {(/>г,<р):Q<t< Cl/t/2(^ + 2),G<
у =	) + -Я-(и' ~-),u(Q, rtq>) = (rcosp/', 0< г <],
( dt гдг___дг г ар д<р <	।
О £	* /2. ffl) = cos^[C, -	, u(l, 0, р) = 0,
ТЧ.’
0< кС,ц12(ц +2),0< <р< я12„С, -const, и(г,г,^/2) = О,
•iftr.OJ.Z’IC,	0<Г<^/2(р +2),0<'^1
:	J  Г "К; Rf<i;
,э ' ди д	id
dt гдг^™ дг^+ггд^и д<р^
<•XV,<z>) =	1,0£ 0><?г/2,
«0,0,«>) = о,«(/,], р)=sin^[Cz	+
:ЛЧ''
Ос t<C,fil2^ + 2),	у>< л 12, С, -const, •
»0,Л0)=о,U(t,r,x /2) = Лс -	+
 ц 	?
°<г<С,д/2(А+2),0<г<1.	; s:
Рейеннем^адачи (30) будет фуявд^ю.ч^и 1:Г1
-Л, .< <3°>
 J~' "‘и •>:)?.
Н'.-цн/уадэдп
4-.|4jr
-. • Tti. XtHin ... v
-I -'	. ,.
Г!:};,-;	91>фК>
/ "Л U, ;
.. •	j • 5<jKj
М* J {«. И .1 Д,,. . •;rpi|	J/f; it’
(32)
азадачи (31)-	J
A ,
где С, - произвольная постоянная.
Численный метод решения
Принимая во внимание анализ решения многомерной задачи
теплопроводности в декартовой системе координат
(см. 2.1.2), выберем для решения, например, задачи (30), разностную
схему с итерациями по коэффициенту теплопроводности н с
расщеплением по пространственным направлениям. Введём в области
интегрирования конечномерное пространство (сетку):
А = {(г/><Р»S)-ri = Ы l = 0,L,hrL = l;p„= hjn, m = 0,M,
h9M-n/2;tn -nr, n = 0,N, tN = Ctft/2(jU + 2)-&,Д-const^
и выпишем систему разностных уравнений:
%WO'U)''	К^,.„ -»/,„)
1 г	2rt h2	.
_^-#.КиМ,)Д +	~ui-},rly {=т=1.Л/-1, > (33)
2г,/?
“о т ~ 0» й, „ = cos^^M [Ct - 2{ji + 2)t"*1 /ft] , т=О, М.
мг1П	г L,Jt!	I	»
Эта подсистема уравнений является частью полной задачи (34), (34 ).
Й1.т ~	+ (и*шУ^и^1.- ~Ы4».)	.
г	2гД2	.
Г1-\/2((и*тУ + (Ц|-1,т)/'К“|,м ~Ч-1,м)>/ = 1>£-1,т = 1,Л/-1,
2гД2
Ям.
(34)
г	2ri
[(Ч*, )" + (“/*»-! )дКил»	i ь ijfTj, m = 1.Л/-1,1
/ p
-< =Лс -2(д+ 2)/л+,///1A,	=0, / = 1,1-1,
«,‘+1 -u?
ur
"l.m
minimax
< e, e- const •
Hl'BliG»: «
’> л;:
(34’)
' С'="л«’/ = 1>£"1’w = 1,W'“1, ” = 0,Af'4i	‘
u°m=r^‘sin^‘</>„,l ~0,L, m = 0,M-	"'	’•'• r ji.s-t
Если в (33) ввести обозначения
'V' " *Q
9 =-*/-^К«‘.я)д +«1.пГ1г/2г/Л£ и0.м =°» di = <«> ;
^.и=со8^[С;-2(^ + 2)Г7А]^,т = 1,М-1, ' А* ""
то соответствующие уравнения перепишутся в виде (26) и для их
решения может быть применён метод прогонки.
Вторая подсистема уравнений (34), решаемая на дтором полушаге,
представлена ниже:	-	—
Л
С-д>.и j«,r-«*:') ..	: ,й
г	2г/й2
К^Г+С^-ГК^'-С-,} , — —•
--------->/ = и-1,И = 1,Л/-1,	Ч( (35)
• 9	* ‘	•	»*.•
^•к=,}^С' -2(д+ 2)f+l	= o,/«l,£-i.
" '	1	- 1 ’ I ;	,	 ' 1' " ! u ?0.i
Обозначим
»-=*,,)' f«,)'lr/2r,=йи.
-2(А+2)'“/*,м =О,/=йП.
Тогда задача принимает вид (28) и решается методом прогоики
Итерации по к заканчиваются, когда для всех / = LZ-], m =
выполняется неравенство (29). Затем полагают			
•*/,« /.« т	min£; max /,ж		£ e,E-const •,
		с	

и переходят к расчёту нового временного слоя, последовательно
меняя п = О 1.
VI.2.3. Сферическая система координат
Постановка и частные аналитические решения задач
.. Рассмотрим несколько постановок смешанных задач
для квазилинейного уравнения теплопроводности с краевыми
условиями	первого рода в области
D = {(t,r,0) :0 < t< С,ц!2{ц + 2), 0< r< 1, 0< 0< яг/2}
сферической системы координат:
Su 1 d . 2 и ди. 1 д . . Л р ди
— = -г—(г и* —) + —;---------(sinflu —),
dt г2 dr dr r2sin# d0 d0
u(0,r,(9) = C;'/Ar^cos^(9,0< r< l,0£ 0< я/2,
u(z,0,(9) = 0,u(/,l,(9) = cos^[C/ -2(a + 2)-]A
0< Z<C;///2(/z + 2),0£ <p<> nil, C, -const,
u(t, r,0) = r^[C- 2^ + 22£]~^, и(/, г, п 12) = 0,
0</<C>/2(^ + 2),0<r< 1	J
и в области
° = {(Лг>0) :0 £ z< Qu/4(/z + 1).Q^ г 1.0^	/r/2f *
(36)

Г7>*. »’ л«* dff
С, -sonstXr,г,0)=0,u<f,г,я!»='%&
.Я>К«С^/4(Я+Д0<Г,*:.к.-;% <я<лои г.5»^п v
ч ' -V ->-0 “л ; ч(37)
Решение задачи (36) есть функция	. _	,
Р8)
...	. .	Z*-..с.ег ч K^uoisfiiwls
а доя задачи (37) оно представимо в виде
- н=АЛ[С,	•**мЛ (39)
Ц • • - ,;тн.'.ни»до w-л
•пиы,-	'	• •Ои'ЛК	HMKHK31W
где С.-произвольная постоянная.	,
-О1.. Ч >	: (V.	Д),	'-Л
Разностная схема и метод решения задач-’'Ж-п.-рфэ
Алгоритм, который будет .предложен и разобран ниже,
применим как для решения разностной схемы, построенной по задаче
(36), так и по (37). Выберем в качестве примера постановку (37) и
введём в области интегрирования конечномерное пространство
(сетку):	.	' j- .	...;v.
fy. ={(г/:rz = ЛЛI = 0,1, hr L = 1; 0т = hem, т = 0>М,
heM = nl2',t" —fit, n~d,N,tN = Clpl^[p. + Y)-^, A —const) ,
где A>0, C, >0-const. '	: ' (V -'>«
Воспользуемся разностной схемой, аналогичнойг уже ранее
рассмотренным для задач в декартовой и цилиндрической системах
координат. В ней предусмотрены итерации по квазилинейному
коэффициенту теплопроводности н расщепйение алгоритма по
пространственным направлениям с тем, чтобы добиться
2^
f.w	)
экономичности при вычислениях (см. 2.1.2, 222). Вначале приведём
разностную схему для первого полушага задачи (37):
~Ы1.т _ГМ/2((иМ,тУ	~й/т)	 ...	
г	' 2rtf .	’ ~
,, //-1/2 [(Ц/.и У* +	А К«/,д ~	)	—-—-	------
--z=l,z-l, « = 1,М-1,
«о,„ = 0Л.« =sAm[C/-4(A + l)r’//z]A^ = MZ

Эта подсистема уравнений является частью более полной задачи .д'*
«/,т ~«/,т _ Г1+1/гКиМ,тУ + (**/,»)'“~ **/,„)
г ".	>.
- ймл.),
= | 2гД2 ’* ’ ’
й0> =0,й£>Л-sin^.fG -4(^ + 1Х^1Ам=М^
С1 sin6>m+V2[«.,r+«)^](^,
Г -т Л 2r,2sin(9A2
: sin6>m_V2[«r-Cix t_—
2J72sin0„^	, ч. ’	s •-   м >«
к - п, J, J = mini: max
7.я>
М*»
< £,£- const
С' =«/*',/=1.1-1. m = и = 0Л~1. . . . ..
’ v',f. 
«/ж	=	m=0,Af,	. г ;
Если в подсистеме уравнений (40) ввести обозначения 
nil')':.
.<	(41)
Н
С si -is«М
0
—
'а^-г^К^У	h>=i'a>	* •••’ ’ №
ULт = siifidJC, -4(д + 1К+17а]« = ’»М-Ь
то соответствующие уравнения перепишутся в виде (26) и для их
решения может быть применён метод прогонки.
Вторая подсистема уравнений (41), решаемая на втором полушаге,
представлена ниже:	Л	. '
' ЬТС
Г	2r,2sin^
Sin0 ,д[(<Г + (< .)*](«/;' -«‘*L>)	—-
т~ y*l' l,m' v l,fHl ' J' l,m I,Hi—’is	j ।
(42)
it.O-
4- .И
2rfsiD0„hg
m = l,JK-l,
uifi	-4(^+l)Z’+7//f^,7 = 1,2,-1.
Обозначим	f	., л u
=-sm^,+l/2[«»+ir +(<y]r/2r(2sin^hj,im =1-аи-ся, '
cm = -sin^^Ki^y- + («*„_,/}T/2r^0mh2e, d_
uim = rfr[C,-Wji + 1Хя+|///]^, u,0 = 0,1 =
Тогда задача (42) принимает вид (28) и для её решения можно
использовать метод прогонки.
Итерации по к заканчиваются, когда для всех I = \, t -1, т = 1,Л/ —1
выполняется неравенство (29). Затем полагают
„«I _ цУ+1
Ut,m
и,4+1 -W*
<.и и!м
*
= mini: max
йе,е -const >,
«Л1
менмТ=Т^1УлЬЧё1У. Н0В0Г° Временного слоя> последовательйо
fc r'.ij'i.j .• f/	Г > с .
УЬЗ.Оценка погрешностичислеиных решений
Для оценки погрешности получаемых решений разностных
задач строят их отображение на выбранное конечномерное
пространство. Удобнее всего это делать, если шаги сетки по
независимым переменным оказываются кратными. В этом случае
отображением будет множество значений сеточной функции, взятых в
соответствующих узлах сетки. Аналогично можно построить
отображение на то же самое конечномерное пространство частного
аналитического решения задачи, которое будем называть следом.
Сравнивая след и отображение численного решения в одной из
эквивалентных норм выбранного конечномерного пространства, по
величине их разности можно судить о погрешности решения
разностной задачи. Если норма разности следа и отображённых
решений разностной схемы стремится к нулю при стремлении шагов
сеткн к нулю, то это означает сходимость численных решений.
Предельное решение можно рассматривать как решение исходной
дифференциальной задачи.	J	,
VL4. Лабораторная работа на ЭВМ	/
ь.		:
•	: Приступая к выполнению лабораторной работы из приложения
(см. ниже) согласно номеру варианта и задания, следует вначале
ознакомиться с содержанием раздела VI. 1. настоящего пособия, если
задание относится к первому варианту, и раздела VI.2., если она
содержится во втором варианте, а затем иайти частное аналитическое
решение смешанной квазилинейной задачи для уравнения
теплопроводности. Для пространственно одномерной задачи типа (19)
в качестве конечномерного пространства, на котором будет
происходить сравнение следа и отображённых численных решении,
рекомендуется выбрать
£>А={(х„7'):х/=/А,/ = 0Л0,ЮЛ = и7’ = Сод/2(й// + д + 2)-А,
Д > 0- const}, а для пространственно двухмерных задач, нагфимер,
подобных (25), -	___
Dh={(Xl,y„,T)-.Xl=lhx,l = 03,5hx = l,ya = lh,,1*0,5,	-
5АЛ =!;?’ = С,/г/4(/z +2)-Д, A>О-const}.	.
Удобно представлять результаты работы в виде таблицы, в
первом столбце которой приведены значения координаты х узлов
сетки во втором столбце - аналитическое решение (след), в третьем
столбце-отображение численного решения, полученного на кратной
сетке на выбранное конечномерное пространство, в четвертом
модуль разности между следом и отображенным решением.
Максимальное по модулю число в четвертом столбце характеризует
норму, введенную как максимум сеточной функции по всем узловым
точкам данного конечномерного пространства. Для пространственно
двухмерных задач рекомендуется оформлять результаты расчётов в
виде трёх таблиц: отдельно для следа, отображённого численного
решения и их разности. По строкам и столбцам каждой таблицы
приводятся шесть координат независимых переменных, принад-
лежащих области интегрирования задачи, а на пересечении строк и
столбцов - значения соответствующей сеточной функции.
В ходе выполнения практических расчётов на ЭВМ при
решении смешанной квазилинейной задачи уравнения
теплопроводности приобретается индивидуальный опыт написания и
отладки вычислительных программ, реализующих усложнённые
варианты алгоритма прогонки, проходит опытную апробацию
принцип расщепления по пространственным направлениям и
итерационный метод учёта квазилинейных зависимостей в исходных
дифференциальных уравнениях. При этом оценка точности
полученных численных решений осуществляется по строгим
математическим правилам, используя норму, введённую в выбранном
конечномерном пространстве.
Процесс работы над заданием предполагает .выполнение
следующих этапов:
1) нахождение частного аналитического решения
квазилинейного уравнения теплопроводности и использование
поставленной математической задачи в качестве модельной при
отладке вычислительных программ;
2) построение следа в выбранном конечномерном пространстве;
спекто^ьную еТтпН-Ие разностной схемы на аппроксимацию и
коэффициентов; ичивость> используя принцип замороженных
Реали^Тва^Г1,ТмР1т,а РеШеНИЯ	задачи и его
высокого уровня; мы’ написанной на алгоритмическом языке
аналитаческогореще^^МЫ С учетом вычисленных значений следа
>	кладки ti/jiiijn .> мошрь
Ж
6)	проведение практических исследований вычислительных
свойств используемой разностной схемы на последовательно
удваиваемых сетках:
а)	соответствует ли порядок сходимости численного
Vi । Решения порядку аппроксимации на решении, установленном в
пункте 3;
б)	выяснение вопроса: насколько сильно зависят
? ! ; Численные результаты от выполнения условий устойчивости и
что будет происходить, когда они нарушаются?
в)	предложить и опытным путем проверить свои
варианты разностных схем.
СПИСОК ЛИТЕРАТУРЫ
1.	Демченко В.В. Уравнения и системы уравнений с
частными производными первого порядка: Учебное пособие. -
2- изд. - М.: МФТИ, 2004. - 116 с.
2.	Демченко В.В. Краевые задачи для обыкновенных
дифференциальных уравнений второго порядка. Метод
прогонки: Учебно-методическое пособие. - М.: МФТИ, 2004. -
56 с.
3.	Годунов С.К., Рябенький В.С. Разностные схемы:
Учебное пособие. — М.: Наука, 1973. — 400 с.
4.	Самарский А.А., Николаев Е.С Методы решения
сеточных уравнений. — М.: Наука, 1978. — 592 с.
5.	Марчук Г.И. Методы вычислительной математики-
М.: Наука, 1980. - 536 с.
183
Xi. ПРИ Л0ЖКШК • •	3 (феты
"	ад»~4» :«Ь.<к-ИИ!НЛал.У
Вариант№1 _ / ..,,>.f.. ,	;Т^. ..t? rfi - , 
Численно и аналитически решить одну из смешанных задач
чг
а
для квазилинейного уравнения теплопроводности с точностью
gslQ-4 и сравнить их значения в одиннадцати равноудаленных
Точках в момент времени Т = 1. j j . .. ы
i:;	.	ii U)IS ' -ЧЬ. Л-
Задание Ms 1	ч
к, Дифференциальная задача
ди д , ди. .	_, .	,
'	—=—(и—),0<t<l,0<x<l;
dt дх дх
и(0,х)=(1+х)77,0<х>1;
u(z,0) = l/(7-6r),0</^i;
и((,1)!=4/{7-:6г),0</$1.>
i‘ ' •	J-.O.’.q-Jt). ьЧх.ЦГ,:
>)! i - М f. .ft
F •«'л-
.M
X- .л‘>»йа
1 feint
fttHOfjii vjwajrRjm*
Задание Ms 2
* мкг Дифференциальная задача 4
«ио .' ; > •ди _ д 2 ди 
мим ,I
=2/75^4?,0<^l.
z	'	v’ Xi
-,»:/•.!	у	...»
Заданием 3	» • •; mi
. ! Дифференциальная задача
:лГГЙ^г*фф|гД;
’ -v -f. :цхкм лоцг;
.о dt
ычл •£.... ... «
:l< мж аеИНоИ£?
ШИВть t*,
	;.>).wr.v гыниетх;
д у2ди, л
& дх(и ar),0<Z£1’ Vx<1;
«М = (1 + х)4/121,0^х<1;
и(^0) = 1/(11-10О2,0<^1;’
«(М) = 16/(11-10/)2, 0<t ^1

Задание №4	% М
Дифференциальная задача жииии:	амин^хрффиД
~-=, 0 < WM < ><4;
dt дх дх
н(О,х) = (1 + х}* fl/25,0 $ х < 1;
н(/,О) = 1 /(5 -14//3)* О < t < 1;
а(М)=V16 /(5 -14//3)^; 0 4J t < I.
Задание №5
Дифференциальная задачи ®^ЕЙ: дан.»к®к5^.-5эффа&
- —=А(иМЭ^£)	Г
t dt дх дх	’	 ’
и{е,х) (1 + х)6/3375, О < х < 1;
u(t, 0) = 170 5 -14/)3,0 < / 51;
ЦМ) = 64/(15-14/М	1.
Задание №6	й»*/ <
Дифференциальная задача	ька-л'- «>;-<
dt г дг1 дг'
]	! w(O,r)==H/9,0iS7*^4;
J	и(/,О)=0,О</< 1;
[	- й ir(r,l)*W~8$j
Задание №7
Дифференциальная задача г-4»екзг
j 3/ г dr' dr	’
j :i u(Q,r)=r/y[7,0<r<l; :	. 1
I 1 «(/,0}= 0,0 <7'51;	;
Ш1
Задание №8
Дифференциальная задача
^=1—(п?2—), 0 < Г 1,0 < г < к
dt г дг дг
и(О>г) = г4/196, 0 £ х
и(/,О) = О,О</<1;
, и(М) = 1/(14-12г)2,0<г<1.	.
Задание М 9	?. •
Дифференциальная задача
^.-L^ru3—), О < t < 1,0 < г < 1;
dt г дг дг
и(0,г)=(г2/6)’/3, 0<г <1;
и(г,О) = О, О < / < 1;
w(/,l) = l/^/6-16z/3, О</< 1.
Задание/® 10	'• 3
Дифференциальная задача
ди 1 д щ ди
~Z) ’ 0 < 10 < г <
dt г дг дг
u(0,r) = r6/4913,0	< 1;
и(/,О) = О, О < t 1;
*0,1)=1/(17Ч6Оэ,Ок/£1.
Заданием 11
Дифференциальная задача
1 д j ди
и(0.г) = г2/11,0^г £1;
«(ЛО) = О,О<^1;
и(М) = 1/(11-1ОГХО</<1.
. <
188
Задание № 12	,!ь
Дифференциальная задача
9м 1 9 2 , ди. л
а7(г и'’0 < **0 <г < U
ot г or дг
u($,r} = rl3>0<r < 1;
и((,О) = О,О<Г <1;
м(/,1) = 1/79-&,0<г^1.
Заданием 13
Дифференциальная задача
ди 19.2irt ди. „
— = ~у—(ги' —), 0 < t < 1,0 < г < I;
dt г2 дг дг
м(0,г) = //225,0<г£1;
u(t,O) = O,Q<t£V,
u(r,l) = l/(15-14O2,0<f£l.
Задание № 14
Дифференциальная задача
ди \ д j уз gw
dt г2 дг дг
0<t S1,0 < г < 1;
w(0,r) = r‘/6859,0<r<l;
w(f,O) = O, О < Г < 1;
»(д1) = 1/(19-18О3, Q<tl 1.
Задание М 15
Дифференциальная задача
0<г<1;
dt г2 дг дг
и(О,г) = г*’/^Т,0<г^Г,
м(г,О) = О,О</<1;
«(г,1)=1/(9-2бг/3)*э,0</51.	(
К9
Вариант №2	С тлио»
Найти аналитическое и численное решения смешанной задачи
для квазилинейного уравнения теплопроводности согласно номеру
задания и сравнить их значения в момент времени Т = 1.
Задание №1
Дифференциальная задача . г у.
dt &V йс/ ЗК J
и(0,х,у) = (1 + х + у)2/13,0<;х,у<П;
u(t,0,y)=(I+y)2/(13-12O, 0<Т<И, 0<
и(М,у) = (2 + у)2 /(13-12/), 0<t< 1, 0< у<1;
и(?,х,0) = (1+х)2/(13-12г), 0<Г<1, 0< х<1;
иа,х,1)=(2 + х)2/(13-120, 0<Т<1, 0<х<1.
J SkiV.SibrX
 >  (У.' ФчЬ
Задание Ns 2
Дифференциальная задача
ди д ( -> ди
Т"=~ и' —
dt &(. дх
0<^1, 0<х,у<1;
и(0,х,у)=(1 + х + у)/3,0<х!у<1;	I
u(t,0,y) = (1 + у) / л/9 -8/, 0<Г <И,0£ у<А;
«(t,l,y)=(2 + y)/^9-8r,	0< ySl;
м(/,х,0) = (1 + x)/V9-8t, 0<Г<1,0< х<1;
w(l,x,l)=(2 + x)/>/9-8t,	Q,<x<l. '
Задание № ' 3 :	!	' w \ , —4- --
Дифференциальная задача' ’
Чэг;г;фс. ;К


Ш)
д (
ди д ( Ц2 ди
dt дх\ дх_
и(0,х,у) = (1 + х+у)4/44А,0^х,у^\-,	-
M0,0,j) = (l + y)4/(21-20r)2, 0</<1, 0< у<1;
и(1,1,у) = (2 + уУ/(21-201)2, 0<t<l, 0<
u(t, х, 0) = (1 + х)4 /(21 - 20/)2, 0 < / < 1, 0 < х < 1;
w(/,x,l) = (2 + x)4/(21-20/)2, Ь</<1, 0<х<1.
|/2 ди} л	М
U лГ ’ ®</^1»-0<х,у<1;

- .-?/
0</£1,0<х,у<1;

Задание №4	’	'
Дифференциальная задача
йи д( 3/2 ди} д( 3/2 ди}
— -— и*— +— и4 ------------
dt cbc^ дх) йу^ йу?
*Uo w(0,x,y) = (1+ х +y)^3/VioF, 0<х,у<1;
и(/,0,у) = (1+у)4/3/(10-28//3)2/3,0</<1, 0<у<1; •
«(/,1,У) = (2 + у)4/3/(10-28t/3)^, О</<1, 0<у<П;
u(t, х, 0) = (I + х)ф /(10 - 28/ / З)*3,0 <t < 1, 0 < х< 1;
м(/,х,1) = (2 + х)^ /(10-28//3)Л 0</$1> 0< х<\


0</<1, 0<^,у<;1;
Задание № 5
Дифференциальная задача
ди д ( уз ди} д ( |Л ди'
— =---- !/*’— | +-Нт ---
dt йг^ дх) йу^ ду)
w(0,x,y) = a + i+y)6/2700050^x,y<l; •
«(Л0,у) = (1 + у)6 /(30 -28/)3, 0<t< 1,0< у<1;
и(/,1,у) = (2+ у)6/(ЗО-28/)3, 0</<1, 0< у^1;
ы(/,х,О) = (1 + х)6/(ЗО-28/)3, 0</^1, 0< х<1;
u(t,x, 1) = (2 + х)6/(30 - 28/)3, 0</£1,0< х<1'
'..v. .<•!
л \г5
гг 16

*‘0 ,Т'': 1 > ОД - - £•'.?

У
ЗаданиеМб	г-J	“
Дифференциальная задача	'	,	;
5w =f»—\
Qt rfrl дг) г2д<р{ д<р) .	. . .O.v
и(0,г,<р) = г2со52<р/7,0^г<}^<^<я/2;
«(/,о,^=о, o<^i,o<w2;
«(/Др) = cos2p/(7-6t), 0</£1, 0<<Р<лг/2;	}
u(t, г,0) = г2 /(7 - 6/), 0 <t < 1, б < г < 1;

и(/,г,ж/2) = 0,0</<1,0<г<1.
s.'ww ишглшннзцуффкГ,
Заданием 7	. ч-,; $	,К; ц '> s
Дифференциальная задача '	 --*-[ ~	 ' гт”---
/	\	" . ?•>>	to
Эи д [ ди] д ди] п * .. '
^=_ГГМ^“ +~^~ «Т- кб<Ъ^1> 0<^<^/2; 0‘и -
dt rdr\ dr J г‘д<р\ dtp)	rl Л “
«(0,r,^) = r2sinV/7,0<г <1, 0<р<^/2; ; " П <’< C v,w
и(/,0,^) = 0,0<^1,0<^<лг/2; ' '	- *Х.ЬПи
	«(/,l,p) = sin2^/(7-6/), б<Г<1,0^^^лг/2;	"1 > и.Од. л)и
'	и(г,г,О) = О, 0</^1, 0<tr<Jir	•. •’
Задание № 8
Дифференциальная задача	. У
^-_Lfr,?auL 5 ( 2 ди }
dt гдг\ дг) г20^ др;
л
хб vS
,()д)и
. u(t,r,n/2) = r2/(7-6/), 0</^1, 0<г<1.
rf- ЖЬп.ЫЦ>.(Ц:иигрффпД
5 f «Л .Л;
*г + * — ' w
:	?•. > г;л
(0<^<я/2
H(0,r,p) = rcos^/>/5, O^r^l,	+ 11 л
и(Л0,9>) = 0,0</^1,0^^><я-/2;	?< 1 г-'/<*•>)«
и(М,Ф)^со8р/^-4/,	О^^л-Д- ’ !> <С) J
и(лг,о)=г/ХГ47>о</^1,о<г<1; Мл)м
и(/,г,^/2) = 0,0</^1,0<г<1.

Задание №9	г i ,
Дифференциальная задача	х<кЧ tiгЛ'!...й!И(’^..я.
д (	д ( , диУ	'
Д'"
w(0,r,^) = rsin^/V5, 0<г£1, 6&р^я/2;
u(t,Q,<p) = О, О <t < 1,0<р<д-/2;	- * А’
4(/,1,р) = sin^/75 -4/, О</< 1, 0<р£я/2;
;	, и(/,г,0) = 0, 0</<15 0<г<1;	’-‘ >
и(/,г,тт/2) = г/-J5 h,0<t<,],0<r<l.
' У< ;
— 1,0</,г<1,0<^<лг/2;r'
Задание Ns 10	RL.w;
f Дифференциальная задача
, du d f (M \	$ Г i/2
1 dt rdr\ dr) r2S^
j u(0,r,(p) = r4cos4^/121, 0<f < 1,
;c• - ; w(Z,O4p) = O, 0<t£1, 0£p<tr/2; n {.	>
a(r,l,p) = cos>/(ll-10r)2,0</^ls0<^^»’/2; Z
u(t,r,Q) = r* /(11 - 10Г)2 Л </<l, 0<r <1; 'b
1	1/(г,г,лг/2) = 0,0<г^1, d<r<l.
V.
v l

Заданием II
Дифференциальная задача
г2и—)+-г— ------fsintfec^- |,0<f^1» 0<r<|;
dr) r2sin6E0l de)

du
dt r2dr
!	0<6»<я-/2;
)	w(0,r,^Hcos25/7,0^r<l,0<^/2; ' *' u 
w(4o,fl) = O,u(V,«-/2) = o, 0<y<l, 0<г<1,0<|9<л-/2;
«(/,1,6»)=cos2^/(7-6t), 0<f<l, 0^в^я-/2; ru
w(/,r,0)=r2/^7 (hc/ib

Заданием 12	,
Дифференциальная задача	;
ди д _______________д----fsin0i3f|,O</<l, 0<г<1;	-:
dt r2cH dr) r2sm05O\ дВ)	v>
0<#<я/2;	, л.
M(O,r,0) = r2sin20/9, (Иг<1, О<05л-/2;	,
и((,0,^)=0,«(Лг,0) = 05 0</<1,0<г<1,0^<я:/2; {
«(/,1,0) = sin2 0/(9 - 8Г), 0<t<l, О<0< л-/2;
u(t,r,7t/2) = r2/(9-3t), 0<t<l, 0<г<1
 >.	. '..у. л.Ш.
Задание № 13
Дифференциальная задача
ди_ д (2 2	[ д ’f
dt r25r\ . dr J r2sin&30V
'	AV и
А» \ ’’rji d'i
sin0«2—0<r<l:
50J	>
O<0or/2;
M(O,r,0) = rcos0/V5, 0^г<1, О<0<яг/2; '
«(r,0,^) = 0,w(?,r,^/2) = 0,0</<l, 0<г<1, О£0<я/2;
u(t,l,0) = cos0/y/5-4t, O<t<l,O<0<ff/2;
w(t,r,0) = r/j5-4t,0<f,2150<r<l.
Задание №14
Дифференциальная задача	' 	’’
9 ( 2 гди\ д С	.Л,;‘ 'Ч
а’л’агГ" -лГййЖ"-* ав>°5».^о<г<Ц ,«
;	O<0<x/2t	'	' ' ’ V>
и(ЛО,0) = ОХГ,г,О) = О,О<Г<1,О<г<1,О^^^/2; *.	,
“(/’1>^) = sin0/V7^6F, 0<Г<1, 0<6»^ff/2;.
“(Аг,л-/2) = Я/>/7-6/,0</Sl,0<г<}(,' л.Нч
,194
Задание № 15
Дифференциальная задача
д ( г у/2 ди А д ( . У2 ди А
-— г2ич	—  sm(9w/2—	0<г<1;
г2аД dr) r2sm0de{ дв)
О<0<я-/2;
и(О,г,0) = r4sirt40/196,O<r <1,О<0<;я/2;
w(/,O,0) = 0,«(/,r,0) = 0,0<t< l,Q<r<\,Q^0<7r/2-,
U(/,l,^) = sin4<9/(14-12/)2,0</$l,0<to/2;
1/(/,г,я-/2) = г4/(14-12/)2,0</<1,0<г<1.	j
: П'ХЛ7П fl i iHdE.ТП1Г
ЗЯНТ/ШТАЬ’ NOriVJUtHH Oil
ПХ..:'i B"5MV«{	 tic
Oft-- -«MjiiT ?’ ,'C 
.ЛГЯ. Uhx
195
Учебное издание	г
ДЕМЧЕНКО' Владимир Владимирович
ВЫЧИСЛИТЕЛЬНЫЙ ПРАКТИКУМ
ПО ПРИКЛАДНОЙ МАТЕМАТИКЕ
Редактор И.А. Волкова
Корректор 0.77. Котова
Подписано в печать 07.09.2007. Формат 60 х 84 '/16. Бумага офсетная.
Печать офсетная. Усл. печ. л. 12,25. Уч,- изд. л. 10,75, Тираж 400 экз.
Заказ № 1231.
Государственное образовательное учреждение
высшего профессионального образования
Московский физико-технический институт
(государственный университет)
141700, Московская обл., г. Долгопрудный, Институтский пер., 9
Для заявок*, тел. (095) 408-58-22 rio@mail.mipt.ru
Отпечатано в полном соответствии с качеством
предоставленных диапозитивов
в ГУП МО «Орехово-Зуевскаятипография»
142603, Московская область, г. Орехово-Зуево, ул. Дзержинского, Д. 1
I
и