Министерствообразования и науки РФ
ГОУ ВПО “УГТУ-УПИ”
Курсовая работа
по “Вычислительной математике”
на тему: “Некоторые дополнительные вычислительныеметоды”
Семестр № 3
Преподаватель Кочнев В.П.
Студент гр. №р-23021д Логиновских М.А.
Номер зачетной книжки17309013
Екатеринбург
2004
_____________________________________________________________________________
Домашнее задание по ________________________________ №________________
№ записи в книге регистрации __________________ датарегистрации ___________200_г.
Преподаватель _________________________________________
Студент _________________________________________ группа №________________
Деканат ФДО _______________
СОДЕРЖАНИЕ
1. Решение систем линейных уравнений …………………………………………………… 3
а) Схема Халецкого ……………………………………………………………………… 3
б) Метод Зейделя и условия сходимости …………………………………………………5
2. Методы решения нелинейных уравнений …………………………………………………6
а) Метод хорд ………………………………………………………………………………. 7
б) Метод Ньютона (метод касательных) ………………………………………………….8
в) Метод итерации ………………………………………………………………………… 9
3. Интерполирование и экстраполирование …………………………………………………11
а) Интерполирование с помощью многочленов …………………………………………11
б) Интерполяционный многочлен Лагранжа …………………………………………….12
в) Интерполяционные многочлены Стирлинга и Бесселя ………………………………13
г) Тригонометрическое интерполирование …………………………………….…………15
д) Интерполяция сплайнами ……………………………………………………..………… 15
4. Численное дифференцирование и интегрирование ……………………………….………16
а) Постановка задачи численного интегрирования ………………………………………16
б) Составные квадратурные формулы ………………………………………………….…17
5. Приближенное вычисление обыкновенных дифференциальныхуравнений …………… 18
а) Метод Рунге-Кутта ……………………………………………………………………… 18
б) Экстраполяционные методы Адамса ……………………………………………………20
в) Метод Милна ……………………………………………………………………………. 20
г) Краевые задачи для обыкновенных дифференциальныхуравнений ………………… 21
6. Приближенные методы решения дифференциальных уравнений счастными производными ………………………………………………………………………………… 21
а) Классификация дифференциальных уравнений второгопорядка …………………… 22
б) Постановка краевых задач ……………………………………………………………… 23
в) Метод конечных разностей (метод сеток) ………………………………………………24
г) Разностные схемы для решения уравнениятеплопроводности ……………………… 25
д) Разностные схемы для решения уравнения колебанияструны ……………………… 26
7. Список литературы ………………………………………………………………………… 27
1. Решение систем линейных уравнений
Системы линейных уравнений (СЛУ) имеют в вычислениях оченьбольшое значение, так как к ним может быть приведено приближенное решениеширокого круга задач. Так, основными источниками возникновения СЛУ являютсятеория электрических цепей, уравнения балансов и сохранения в механике,гидравлике и т.д. Существует несколько способов решения таких систем, которые восновном делятся на два типа: 1) точные методы, представляющие собойконечные алгоритмы для вычисления корней системы, 2) итерационные методы,позволяющие получать корни системы с заданной точностью путем сходящихсябесконечных процессов. Заметим, что даже результаты точных методов являютсяприближенными из-за неизбежных округлений. Для итерационных процессов такжедобавляется погрешность метода.
Пример системы линейныхуравнений: />
Или вматричном виде: />,
где />матрица коэффициентов системы;
/> -вектор неизвестных; /> - векторсвободных членов.
Схема Халецкого
Запишем систему линейных уравнений в матричном виде: />,
где A=[aij] – квадратная матрица порядка n и
/>, /> - векторы-столбцы.
Представим матрицу A в виде произведения нижней треугольной матрицы B=[bij] и верхнейтреугольной матрицы C=[cij]с единичной диагональю />, где
/> и />.
Тогда элементы bij и cij определяютсяпо формулам
/> и />
Отсюда искомый вектор x можетбыть вычислен из уравнений /> и />.
Так как матрицы B и C – треугольные, то системы легко решаются:
/> и />
Из этих двух формул видно, что числа yi выгодно вычислять вместе с коэффициентами cij.Этот метод получил название схемы Халецкого. В схеме применяется обычныйконтроль с помощью сумм. Если матрица A – симметрическаяaij=aji, то/>
Пример. Решить систему />
Решение.
В первый раздел таблицы впишем матрицу коэффициентовсистемы, ее свободные члены и контрольные суммы. Далее так как /> />,то первый столбец из раздела 1 переносится в первый столбец раздела II. Чтобыполучить первую строку раздела II, делим все элементы первой строки раздела Iна элемент/>, в нашем случае на 3.
Имеем: />; />; />; />; />.
Переходим к заполнению второго столбца раздела II, начинаясо второй строки. Пользуясь формулами, определяем />:/>; />; />.
Далее определяя по формулам, заполняем вторую сетку дляраздела II:
/>/>
/>
/>
Затем переходим к третьему столбцу, вычисляя его элементы /> и /> по формулам и т.д., пока небудет заполнена вся таблица раздела II. Таким образом, заполнение раздела IIпроисходит способом “елочки”: столбец — строка, столбец — строка и т.д.
В разделе Ш, пользуясь формулами, определяем /> и />.
Текущий контроль осуществляется с помощью столбца ∑,над которым производятся те же действия, что и над столбцом свободных членов.
/>
/>
/>
/>
/>
/>
/>
/>
/>
/> I
/>
/>
/>
/>
/>
/> 3 1 -1 2 6 11 I
/>
/>
/>
/>
/>
/> -5 1 3 -4 -12 -17 I
/>
/>
/>
/>
/>
/> 2 1 1 1 3 I
/>
/>
/>
/>
/>
/> 1 -5 3 3 3 -1 II
/>
/>
/>
/>
/>
/> 3 0.333333 -0.333333 2 2 3.666667 II
/>
/>
/>
/>
/>
/> -5 2.666667 -0.25 0.25 -0.75 0.5 II
/>
/>
/>
/>
/>
/> 2 -0.666667 2 -1.25 -1.75 -2 II
/>
/>
/>
/>
/>
/> 1 -5.333333 6 2.5 3 4 III
/>
/> 2 1 III
/>
/> -0.75 -1 III
/>
/> -1.75 2 III
/>
/> 3 3
Метод Зейделя и условия сходимости
Этот метод представляет собоймодификацию метода простой итерации. Его смысл заключается в том, что привычислении (k+1)-го приближения неизвестной xi учитываются ужевычисленные ранее (k+1)-е приближения x1, x2, ..., xi-1.Пусть дана приведенная линейная система /> (i = 1, 2, …n). Выберем произвольноначальные приближения корней />,стараясь, чтобы они в какой-то мере соответствовали искомым неизвестным x1,x2, x3, ..., xn. Предположим, что k-еприближение /> корней известно, тогда всоответствии с идеей метода будем строить (k+1)–е приближение по следующимформулам: />
Обычно процесс Зейделя сходится быстрее, чем метод простойитерации. Бывает, что процесс Зейделя сходится, когда простая итерациярасходится и т.п. Правда, бывает и наоборот. Во всяком случае, достаточныеусловия сходимости для метода простой итерации достаточны и для сходимостиметода Зейделя. То есть процесс итерации сходится, если выполнено одно изусловий
1) /> или 2) />.
Пример. Методом Зейделя решить систему уравнений />
Решение. Приведем эту систему к виду, удобному для итерации,/>
В качестве нулевых приближений корней возьмем: />; />; />.
Применяя процесс Зейделя, последовательно получим:
/>
/> и т.д.
Результаты вычислений с точностью до четырех знаков помещеныв таблице:
/>
/>
/>
/> 1,2000 0,0000 0,0000 1 1,2000 1 ,0600 0,9480 2 0,9992 1,0054 0,9991 3 0,9996 1.0001 1,0001 4 1 ,0000 1,0000 1,0000 5 1 ,0000 1,0000 1,0000
Точные значения корней: />.
2. Методы решения нелинейных уравнений
Как известно, далеко не всякое уравнение f(x)=0 можно решить точно, т.е. не всегда можно найти число /> такое что f(/>)≡0. В первую очередьэто относится к трансцендентным уравнениям. Кроме того, даже для алгебраическихуравнений степени выше четвертой не существуют формулы, выражающей их решениячерез коэффициенты уравнения при помощи арифметических операций и извлечениекорней. Для уравнений третьей и четвертой степени формулы для отыскания корнейсуществуют, но они настолько сложны, что практически не применяются. Поэтомубольшое значение имеет приближенное вычисление корней уравнения f(x)=0. Для этого существуетмножество методов некоторые, из которых мы рассмотрим.
Метод хорд
Пусть дано уравнение f(x)=0, где функция f(x) определена инепрерывна на интервале
[a, b] и f(a)f(b)0. Разделим отрезок [a, b] в отношении — f(a):f(b). Это даст намприближенное значение корня x1 = a + h1, где
/>.
Далее этот прием применяем к одному из отрезков [a, x1]или [x1, b], на концах которого функция f(x) имеет противоположныезнаки. Аналогично находим второе приближение x2 и т.д. Геометрическиэтот способ эквивалентен замене кривой y = f(x) хордой, проходящей через точкиА[a, f(a)] и B[b, f(b)].
/>
f(b)
ξ f(a)
0 x ξ x3 x2 x1 b=x0 a=x0 x1 x2 b
/>/>/>
0 f(a) x />/>/>/> a
f(b)/> /> /> />
Действительно, уравнение хорды АВ имеет вид />
При х = х1 и y = 0, получим />
Полагая, что на отрезке [a, b] вторая производная f''(x)сохраняет постоянный знак, метод хорд сводится к двум различным вариантам.
Из рис. 1 видно, что конец а неподвижен и последовательныеприближения: x0=b;
/>
образуют ограниченную монотонно убывающуюпоследовательность, причем a
Из рис. 2 видно, что неподвижен конец b и последовательныеприближения: x0=a;
/>
образуют ограниченную монотонно возрастающуюпоследовательность, причем
x0
Таким образом, для вычисления корня уравнения имеем дверазличные вычислительные формулы. За неподвижный конец выбираем тот конец, длякоторого знак функции f(x) совпадаетсо знаком второй производной f''(x).
Пример. Найти положительный корень уравнения /> с точностью до 0,002.
Решение. Прежде всего отделяем корень. Так как /> и />, то искомый корень /> лежит в интервале />. Полученный интервал велик,поэтому разделим его пополам. Так как /> то/>. Последовательно применяяформулы, будем иметь:
/>
/>
/>
/>
/>
/>
Так как /> и при /> имеем />, то можно принять: />
Таким образом, />, где />. Заметим, что точный кореньуравнения есть />.
Метод Ньютона (метод касательных)
Пусть корень ξ уравнения f(x)=0, отделен на отрезке[a, b], причем первая и вторая производные f'(x) и f''(x) непрерывны исохраняют определенные знаки при />. Найдякакое-нибудь n-ое приближение корня />, мыможем уточнить его по методу Ньютона следующим образом. Пусть ξ=xn+hn,где hn — величина малая. Отсюда по формуле Тейлора получим: f(xn+ hn) ≈ f(xn)+hn f¢(xn)=0. Следовательно, />. Подставив полученноевыражение в формулу ξ=xn+hn, найдем следующеезначение корня: />
/>Графическое нахождение корняметодом Ньютона (рис. 3).
/>
Если в качестве начального приближения выбрать точку х0=В0, то процесс быстро сходится. Если же выбрать точку х0=А, то х1/>[a, b], и процесс нахождениякорня расходится. Рекомендуется: в качестве х0выбирать такую точку,где f(x0)f''(x0)>0.
Пример. Вычислить методом Ньютона отрицательный кореньуравнения
/> с пятью вернымизнаками.
Решение. Полагая в левой части уравнение /> получим />. Следовательно, искомый корень/> находится в интервале />. Сузим найденный интервал.Так как /> то />. В этом последнем интервале/> и />.Так как /> и />, то можем принять заначальное приближение />. Последовательныеприближения /> вычисляем по следующейсхеме:
/>
/>
/>
/>
/> -11 3453 -5183 0,7 1 -10,3 134,3 -4234 0,03 2 -10,27 37,8 -4196 0,009 3 -10,261 0,2 - -
Останавливаясь на />,проверяем знак значения />. Таккак />, то, /> и любое из этих чисел даетискомое приближение.
Метод итерации
Заменим уравнение f(x)=0 эквивалентным x=φ(x). Выберем некоторое начальное приближение x0и вычислим дальнейшие приближения по формулам x1=φ(x0), x2=φ(x1), …, xn=φ(xn-1). Еслипоследовательность xn имеет предел, тоитерационный процесс
xn= φ(xn-1) (n=1,2, …) называется сходящимся. Пусть функция φ(x)непрерывна. Переходя к пределу в равенстве xn=φ(xn-1), получим
/> Следовательно, /> является корнем уравнения x=φ(x) и может быть вычислен поформуле xn= φ(xn-1) (n=1, 2, …) с любой точностью. Для данного методасуществуют две теоремы:
Теорема 1. Пусть корень /> уравненияx=φ(x), а также егопоследовательные приближения x0, xn= φ(xn-1) (n=1, 2, …) содержатся в интервале [a,b] и /> на [a, b]. Тогда справедливы утверждения:итерационный процесс xn= φ(xn-1) сходится к корню уравнения />; />; />.
Следствие 1. Если требуется найти корень с точностью ε,то кончаем итерационный процесс тогда, когда />.
Следствие 2. Так как />=φ(/>) и /> =φ(xn-1), то />-xn= φ(/>)-φ(xn-1). По теоремеЛагранжа />. Из этого следует, что еслиφ'(x)>0 на (a, b), то последовательные приближения xn=φ(xn-1) (n=1, 2, …) сходятся к корню монотонно; если φ'(x)
Теорема 2. Если /> на [a, b], а корень /> и начальное приближение x0находятся на более узком отрезке [α, β],где />, то справедливы заключения теоремы1.
Привести уравнение f(x)=0 к виду x=φ(x) таким образом, чтобы получить сходящийся итерационныйпроцесс, можно различными способами. Рассмотрим два из них:
1) уравнение f(x)=0равносильно при λ≠0 уравнению λf(x)=0 и уравнению x= λf(x)+x. Обозначимλf(x)+xчерез φ(x), получим x= φ(x). Параметр λ подберем так, чтобы функция φ'(x)= λf'(x)+1на [a, b] была по модулю меньшеединицы.
2) если />, тоитерационный процесс расходится. Заменим уравнение x=φ(x) эквивалентным ему уравнением x=ψ(x), где ψ(x) – функция, обратнаяфункции φ(x). Так как />,то итерационный процесс xn=ψ(xn-1) будет сходящимся.
Пример. Методом итерации найти корень уравнения 5x-8lnx=8 с точностью 0,01.
Решение. Запишем уравнение в виде /> и построим соответствующиеграфики:
/>
Уравнение имеет два корня: />.За начальные приближения возьмем z0=0,5 и x0=3,5. Для уточнения запишем />. Здесь
/> Следовательно,итерационный процесс сходится. Погрешность оценим по формуле /> />,результаты вычислений приведены в таблице:n x 1+lnx
/>
/>
0
1
2
3
4
3,5
3,605
3,651
3,672
3,682
2,253
2,282
2,295
2,301
3,605
3,651
3,672
3,682
------
0,105
0,046
0,021
0,010
Так как φ’(z0)≈3>1,то итерационный процесс расходится. Найдем функцию />,обратную функции φ(x). Так как />, то итерационный процесс /> будет сходится. />, результаты вычисленийприведены в таблице:n
zn
/>
/>
/>
0
1
2
0,5
0,503
0,503
-0,688
-0,686
------
0,503
0,504
------
------
0,0015
0,0005
3. Интерполирование и экстраполирование
Задача интерполирования состоит в том, чтобы по значениямфункции f(x) в нескольких точках отрезка восстановить ее значения в остальныхточках данного отрезка. Разумеется, такая постановка задачи допускает скольугодно много решений. Задача интерполирования возникает, например, в томслучае, когда известны результаты измерений yk = f(xk)некоторой физической величины f(x) в точках xk, k = 0, 1,…, n итребуется определить ее значение в других точках. Интерполирование используетсятакже при необходимости сгущения таблиц, когда вычисление значений f(x) поточным формулам трудоемко. Иногда возникает необходимость приближенной замены (аппроксимации)данной функции (обычно заданной таблицей) другими функциями, которые легчевычислить. При обработке эмпирических (экспериментальных) зависимостей,результаты обычно представлены в табличном или графическом виде. Задачазаключается в аналитическом представлении искомой функциональной зависимости,то есть в подборе формулы, корректно описывающей экспериментальные данные.
Интерполирование с помощью многочленов
Пусть функциональная зависимость задана таблицей y0 = f(x0);…, y1= f(x1); …, yn= f(xn).Обычно задача интерполирования формулируется так: найти многочлен P(x) = Pn(x) степени не выше n, значениякоторого в точках xi (i = 0, 1 2,…,n) совпадают со значениями данной функции, то есть P(xi) = yi. Геометрически это означает, что нужнонайти алгебраическую кривую вида /> проходящую через заданную систему точек Мi(xi, yi) (см. рис. 4). Многочлен Р(х) называетсяинтерполяционным многочленом. Точки xi (i = 0, 1, 2,…, n) называются узлами интерполяции.
/>
Для любой непрерывной функции f(x) сформулированная задачаимеет единственное решение. Действительно, для отыскания коэффициентов а0,а1, а2 ,…, аn получаем систему линейных уравнений/> определитель которойотличен от нуля, если среди точек xi (i = 0, 1, 2,…, n) нетсовпадающих. Решение системы можно записать различным образом. Однако наиболееупотребительна запись интерполяционного многочлена в форме Лагранжа или в формеНьютона.
Интерполяционныймногочлен Лагранжа
Пусть на отрезке [a,b] некоторая функцияf(x)задана лишь в некоторых точках />,т.е. известны ее значения />,которые, собирают в таблицу:
x
x0
x1 ...
xn
f(x)
y0
y1 ...
yn
Кроме того, пусть задана некоторая точка />. Построим по таблицеследующий многочлен: />.
Этот многочлен называется многочленом Лагранжа.
Его основные свойства:
1) это — многочлен степени />;
2) />, т.е. многочлен Лагранжаимеет в точках /> те жезначения, что и функция />;
3) если фиксировать любое число /> то окажется выполненнымнеравенство />
где /> научастке />, т.е.число /> ограничиваетпроизводную />гопорядка функции />.
Сказанное означает, что если функция /> задана своей таблицей итребуется найти значение /> где-тов промежуточной точке c, то можно по таблице построить многочленЛагранжа и его значение в этой точке принять за значение функции. Отысканиепромежуточного значения функции называется интерполяцией; когда этоделается с помощью многочлена Лагранжа, то говорят об интерполяционноммногочлене Лагранжа или об интерполяции по Лагранжу.
Пример. Построить интерполяционный многочлен Лагранжа дляфункции заданной таблицейx 1 2 3 5 y 1 5 14 81
И найти значение функции при x=4.
Решение. Используя формулу Лагранжа найдем: />
После некоторыхпреобразований получим /> Тогда f(4)≈L3(4)=36,5.
Интерполяционные многочлены Стирлинга и Бесселя
Взяв среднее арифметическое первой и второй интерполяционныхформул Гаусса /> />
/> и /> />
/>, получим формулуСтирлинга />
где />.
Легко видеть, что /> при />.
Кроме формулы Стирлинга, часто употребляется формулаБесселя. Для вывода этой формулы воспользуемся второй интерполяционной формулойГаусса /> />
/>.
Возьмем /> равностоящихузлов интерполирования /> с шагом/>, и пусть /> — заданные значения функции/>.
Если выбрать за начальные значения /> и />, то, используя узлы />, будем иметь:
/>/>.
Примем теперь за начальные значения /> и /> и используем узлы />. Тогда />, причем соответственноиндексы всех разностей в правой части предыдущей формулы возрастут на единицу.Заменив в правой части этой формулы /> на /> и увеличив индексы всехразностей на 1, получим вспомогательную интерполяционную формулу:
/>/>.
Взяв среднее арифметическое формул, после несложныхпреобразований получим интерполяционную формулу Бесселя
/>
/>
где />.
Интерполяционная формула Бесселя, как следует из способаполучения ее, представляет собой полином, совпадающий с данной функцией /> в /> точках />.
Тригонометрическое интерполирование
Пустьфункция f(х) представлена на некотором отрезке [0, 2p]таблицей значений f(хi) в
равноотстоящихузлах хi =2p(i-1)/(2N+1),i =1, 2, ...,2N+1. Тогда тригонометрическиминтерполирующим многочленом назовем многочлен степени mвида:
/>.
Задачатригонометрической интерполяции состоит в построении тригонометрическогополинома, который бы наиболее полно удовлетворял условиям Рm(хi)= f(хi) для любогоi=1,2, ..., 2N+1.
Можно показать, что решением этой задачи являетсяполином именно того вида, коэффициенты которого вычисляют по следующимформулам:
/>;
/>;
/>.
Интерполяция сплайнами
Пусть отрезок [a, b] разбит на n равных частей [xi,xi+1], где xi=a+ih, i=0, ..., n, xn=b,
h=(b-a)/n.
Сплайном называется функция, которая вместе с несколькимипроизводными непрерывна на всем заданном отрезке [a, b], а на каждом частичномотрезке [xi, xi+1] в отдельности является некоторымалгебраическим многочленом.
Максимальная по всем частичным отрезкам степень многочленовназывается степенью сплайна, а разность между степенью сплайна и порядкомнаивысшей непрерывной на [a,b] производной — дефектом сплайна.
На практике широкое распространение получили сплайны третьейстепени, имеющие на [a, b] непрерывную, по крайней мере, первую производную.Эти сплайны называются кубическими и обозначаются S3(x).
Пусть на отрезке [a, b] в узлах сетки D заданы значения некоторой функции
fi =f(xi), i=0, ..., n.
Интерполяционным кубическим сплайном S3(x)называется сплайн
S3(x)=аi0 +аi1(x — xi)+аi2(x — xi)2+аi3(x — xi)3, xÎ[xi,xi+1], удовлетворяющий условиям
S3(xi)=f(xi), i=0, ..., n.
Данный сплайн на каждом из отрезков [xi, xi+1],i=0, ..., n-1 определяется четырьмя коэффициентами, и поэтому для егопостроения на всем промежутке [a, b] необходимо определить 4n коэффициентов.Для их однозначного определения необходимо задать 4n уравнений.
Условие S3(xi)=f(xi), i=0, ...,n дает 2n уравнений, при этом функция S3(xi), удовлетворяющаяэтим условиям, будет непрерывна во всех внутренних узлах.
Условие непрерывности производных сплайна />, r=1,2 во всех внутреннихузлах xi, i=1, ..., n-1 сетки Dдает 2(n-1) равенств.
Вместе получается 4N-2 уравнений.
Два дополнительных условия обычно задаются в видеограничений на значение производных сплайна на концах промежутка [a, b] иназываются краевыми условиями.
Наиболее употребительны следующие типы краевых условий:
а) S'3(а)=f'(а),S'(b)=f'(b);
б) S"3(а)=f"(а),S"(b)=f"(b);
в) />;
г) S'''3(xp+0)=S'''3(xp-0), р =1, n-1.
4. Численное дифференцирование и интегрирование
Если функция f(x)заданна аналитически ее первообразная F(x) является элементарной функцией, то /> вычисляется по формулеНьютона-Лейбница: /> В тех случаях, когдафункция f(x) заданааналитически, но ее первообразная не является элементарной функцией илиотыскать ее сложно, а также в случае, когда функция f(x) задана графически или таблично, для вычисления /> применяются приближенныеметоды.
Постановка задачи численного интегрирования
Задача численного интегрирования функции заключается ввычислении определенного интеграла на основании ряда значений подынтегральнойфункции. Численное вычисление однократного интеграла называется механическойквадратурой. Обычный прием механической квадратуры состоит в том, чтоданную функцию f(x) нарассматриваемом отрезке [a, b] заменяютинтерполирующей или аппроксимирующей функцией φ(x)простого вида, а затем приближенно полагают: /> Функцияφ(x) должна быть такова, чтобы интеграл /> вычислялся непосредственно.Если функция f(x) заданнааналитически, то ставится вопрос об оценке погрешности. Пусть для функции y=f(x) известныв n+1 точках x0, x1, x2, …, xn отрезка [a, b] соответствующие значения f(xi)=yi (i=0, 1, 2, …, n). Требуетсяприближенно найти /> По заданнымзначениям yi построим полином Лагранжа />, где
Пn+1(x)=(x-x0)(x-x1)…(x-xn),причем Ln(xi)=yi (i=0, 1, 2, …, n). Заменяя функцию f(x) полиномом Ln(x), получим равенство /> гдеRn[f] – ошибкаквадратурной формулы. Отсюда получаем приближенную квадратурную формулу
/> где /> (i=0,1, 2, …, n). Для вычисления Ai заметим, что
1) коэффициенты Ai при данном расположении узлов не зависят от выбора функции f(x);
2) для полинома степени n полученная формула – точная, так как тогда Ln(x)=f(x); следовательно,формула /> - точная при y=xk (k=0,1, 2, …, n), т.е. Rn[xk]=0 при k=0, 1, …, n. Полагая y=xk (k=0, 1, 2, …, n), получимлинейную систему из n+1 уравнений /> - где /> (k=0,1, …, n), из которой можно определить коэффициенты A0, A1, …, An.
Составные квадратурные формулы
Приведем ряд простейших квадратурных формул, используемых впрактике численного интегрирования функции f(x) на некотором интервале [a, b], разбитого на n равных отрезков точками a0=a, a1=a+h, a2=a+2h, …, an=a+nh+b, где n=0,1, …, k и/> Положим f(xn)=yn=f(a+nh).
Формула прямоугольников: />
Погрешность формулы определяется выражением
/> где />
Формула трапеций: />
Погрешность формулы определяется выражением
/> где />
Формула Симпсона: /> где />
Погрешность формулы определяется выражением
/> где />
Если длина интервала [a, b] велика для применения простейших квадратурных формул, топоступают следующим образом:
1) интервал [a, b]разбивают точками xi, /> на n интервалов по некоторому правилу;
2) на каждом частичном интервале [xi,xi+1] применяют простейшуюквадратурную формулу, находят приближенное значение интеграла />
3) из полученных выражений Qi составляют (отсюда и название составная формула)квадратурную формулу для всего интервала [a, b];
4) абсолютную погрешность R составной формулы находят суммированием погрешностей Ri на каждом частичноминтервале.
5. Приближенное вычисление обыкновенныхдифференциальных уравнений
Обыкновенным дифференциальным уравнением называетсяравенство />, в котором/> - независимаяпеременная, изменяющаяся в некотором отрезке />, а /> - неизвестная функция от />, которую и надо найти.Различают два типа обыкновенных дифференциальных уравнений — уравнения безначальных условий и уравнения с начальными условиями. Уравнения безначальных условий — это как раз то, что было только что определено. А уравнениес начальными условиями — это записанное выше уравнение относительно функции/>, но в которомтребуется найти лишь такую функцию />,которая удовлетворяет при некотором /> следующим условиям:
/>, т.е. вточке /> функция />и ее первые /> производных принимаютнаперед заданные значения. В этой ситуации число />называется порядкомуравнения.
Метод Рунге-Кутта
Изложим идею метода на примере: />
Интегрируя это уравнение в пределах от x до x + h (0 котороепосредством последнего интеграла связывает значения решения рассматриваемогоуравнения в двух точках, удаленных друг от друга на расстояние шага h. Дляудобства записи данного выражения используем обозначение
∆y=y(x+h)–y(x) и замену переменной интегрирования t=x+ah. Окончательно получим:
/> Указавэффективный метод приближенного вычисления интеграла в выражении />, мы получим при этом одноиз правил численного интегрирования уравнения />
Постараемся составить линейную комбинацию величин ji, i = 0, 1, ..., q, котораябудет являться аналогом квадратурной суммы и позволит вычислить приближенноезначение приращения Dy: /> где
/>
Метод четвертого порядка для q= 3, имеет вид /> где
/>
Особо широко известно другоевычислительное правило Рунге-Кутта четвертого порядка точности: /> где
/>
Метод Рунге-Кутта имеетпогрешность четвертого порядка (~ h4 ).
Правило Рунге. Если приближенный метод имеет порядокпогрешности m, то погрешность можно приближенно оценить по формуле />
В формуле O(xi) – главный член погрешности, />и /> - приближенные решения вточке xi, найденные с шагом h и 2h соответственно.
Экстраполяционные методы Адамса
Широко распространенным семейством многошаговых методовявляются методы Адамса. Простейший из них, получающийся при />, совпадает с рассмотреннымранее методом Эйлера первого порядка точности. В практических расчетах чащевсего используется вариант метода Адамса, имеющий четвертый порядок точности ииспользующий на каждом шаге результаты предыдущих четырех. Именно его иназывают обычно методом Адамса. Рассмотрим этот метод.
Пусть найдены значения /> вчетырех последовательных узлах />. Приэтом имеются также вычисленные ранее значения правой части />. В качествеинтерполяционного многочлена /> можновзять многочлен Ньютона. В случае постоянного шага /> конечныеразности для правой части в узле /> имеютвид />.
Тогда разностная схема четвертого порядка метода Адамсазапишется в виде />.
Сравнивая метод Адамса с методом Рунге — Кутта той жеточности, отмечаем его экономичность, поскольку он требует вычисления лишьодного значения правой части на каждом шаге. Но метод Адамса неудобен тем, чтоневозможно начать счет по одному лишь известному значению />. Расчет может быть начатлишь с узла />. Значения /> необходимые для вычисления />, нужно получить каким-либодругим способом, что существенно усложняет алгоритм. Кроме того, метод Адамсане позволяет изменить шаг /> впроцессе счета; этого недостатка лишены одношаговые методы.
Метод Милна
Пусть на отрезке [a, b] требуется найти численное решение дифференциальногоуравнения /> с начальным условием />. Разобьем отрезок [a, b] на n равных частей точками />,где h=(b-a)/n – шаг интегрирования. Используя начальные данные, находимкаким-либо способом последовательные значения /> искомойфункции y(x). Таким образом,становится известным />. Приближения /> и /> для следующих значений /> последовательно находятсяпо формулам Милна
/> – где />.
Абсолютная погрешность значения /> приближенноравна />.
Пример. Дано дифференциальное уравнение y’=y-x, удовлетворяющие начальномуусловию x0=0, y(x0)=1,5. Вычислить с точность до 0,01 значениерешения этого уравнения при x=1,5.
Решение. Выберем начальный шаг вычисления. Из условия h4i
xi
yi
y’i=f(xi, yi)=yi-xi
/>
y'i= f(xi, yi)=yi-xi
/>
εi
1
2
3
4
5
6
0,25
0,50
0,75
1,00
1,25
1,50
1,5000
1,8920
2,3243
2,8084
1,5000
1,6420
1,8243
2,0584
3,3588
3,9947
4,7402
2,3588
2,7447
3,2402
3,3590
3,9950
4,7406
7*10-5
10-5
1,4*10-5
Получаем ответ y=(1,5)=4,74.
Краевые задачи для обыкновенных дифференциальныхуравнений
На практике приходится часто решать задачи, когда условиязадаются при двух значениях независимой переменной (на концах рассматриваемогоотрезка). Такие задачи, называемые краевыми, получаются при решении уравненийвысших порядков или систем уравнений. Стандартная постановка краевой задачидля обыкновенных дифференциальных уравнений выглядит следующим образом
/>, а дополнительныеусловия ставятся более, чем в одной точке отрезка интегрирования уравнений (вэтом случае порядок системы не может быть меньше второго): />, />, />.
Общая классификация методов решения краевых задач:существуют точные, приближенные и численные методы.
6. Приближенные методы решения дифференциальныхуравнений с частными производными
Кроме обычных дифференциальных уравнений существуют такназываемые дифференциальные уравнения с частными производными. Далее они будутрассмотрены более подробно.
Классификация дифференциальных уравнений второгопорядка
Рассмотрим уравнение второго порядка />, где /> - функции /> и />. Говорят, что указанноеуравнение в области /> принадлежитгиперболическому типу, если в этой области />.Если />, то уравнение в области /> принадлежит параболическомутипу. Если />, то уравнение принадлежитэллиптическому типу.
Уравнение /> называетсяканоническим уравнением гиперболического типа.
Уравнение /> называетсяканоническим уравнением параболического типа.
Уравнение />называетсяканоническим уравнением эллиптического типа.
Дифференциальное уравнение /> называетсяуравнением характеристик уравнения />.
Если последнее уравнение гиперболического типа, то уравнениехарактеристик имеет два интеграла: /> т.е.существуют два семейства вещественных характеристик.
С помощью замены переменных />,/> дифференциальное уравнение /> приводится к каноническомувиду: />. Для уравненияпараболического типа оба семейства характеристик совпадают, т.е. уравнениехарактеристик дает лишь один интеграл />.
В этом случае осуществляем замену переменных />, />, где /> — какая-нибудь функция, длякоторой />. После замены переменныхполучаем уравнение />. Для уравненияэллиптического типа интегралы уравнения характеристик имеют вид />, где /> и /> — вещественные функции.
Полагая /> и />, приводим уравнение /> к виду />.
Постановка краевых задач
Классическим решением краевой задачи называются всякафункция, удовлетворяющая дифференциальному уравнению в каждой точке внутриобласти задания этого уравнения и непрерывная в рассматриваемой области,включая границу. Соответствующую постановку краевой задачи называютклассической. Существует несколько таких задач:Задача Коши для бесконечной области. Рассмотрим эту задачу на примере
уравнения колебания струны и уравнения теплопроводности.
Рассмотрим процесс колебания тонкой бесконечной струны поддействием непрерывно распределенной внешней силы с плотностью f.Предположим, что сила действует в одной плоскости – плоскости колебания струны(x, u), а струна являетсягибкой упругой нитью. Пусть величина натяжения, возникающая в струне вследствиеее изгиба, подчиняется закону Гука, а сами колебания достаточно малы. Тогдавеличина смещения u (x, t) удовлетворяет уравнению колебания струны: />. Для однозначности процессанеобходимо задать еще начальное смещение и начальное распределение скоростей.Математически это соответствует заданию начальных условий: />. Требуется найтиклассическое решение уравнения, удовлетворяющие начальным условиям.Сформулированная таким образом задача называется задачей Коши длягиперболического уравнения.
Исследуем теперь процесс распределения температуры в тонкомбесконечном стержне. Предполагается, что тепловой поток подчиняется законуФурье, а изменение температуры тела пропорционально количеству теплоты,сообщаемой телу. Предположим, что внутри стержня может выделяться и поглощатьсятеплота, характеризуемая плотностью тепловых источников f.Тогда распределение температуры в стержне описывается уравнениемтеплопроводности: />. Для однозначногозадания процесса необходимо указать начальное распределение температуры. Этосоответствует заданию начального условия: />.Требуется найти классическое решение уравнения, удовлетворяющие начальнымусловиям. Сформулированная таким образом задача называется задачей Коши дляпараболического уравнения.Стационарная задача (задача без начальных данных). Рассмотрим установившийся
режим распределения температуры в ограниченной тонкойпластине произвольной формы с гладкой границей. Пусть функция u(x, y) выражает температуру каждойточки пластины. При обычных законах распространения тепла функция u(x, y) удовлетворяетуравнению Пуассона: />, где функция азадает плотность тепловых источников пластины. В случае отсутствия источника (f=0) данное уравнение называется уравнением Лапласа: />. Для однозначного описанияпроцесса необходимо задать тепловой режим на границе пластины. Это может бытьсделано с помощью задания распределения температуры на границе илираспределения теплового потока. Возможен также режим теплового равновесияизлучающего тела с окружающей средой.
В зависимости от теплового режима на границе получаются триграничных условия для функции u(x,y). Пусть Г – граница рассматриваемой области D – определения уравнения Лапласа. Математическаяформулировка граничных условий может быть задана в следующем виде:
граничное условие Iрода: />;
граничное условие II рода: />;
граничное условие III рода: />.
Производная берется по внешней нормали к кривой Г; λ>0– коэффициент теплопроводности; φ0, φ1, φ2– заданные на Г функции, причем φ2 есть произведениекоэффициента теплопроводности на температуру внешней среды, соприкасающейся стелом.
Таким образом, краевая задача заключается в том, чтобы найтиклассическое решение уравнения Пуассона или Лапласа, удовлетворяющее одному изграничных условий.Смешанная краевая задача. Рассмотрим задачу распространения тепла в тонком
стержне единичной длины. Поместим один из концов в точку x=0, а другой – в точку x=1. Распределениетемпературы в таком стержне в течение некоторого интервала времени 0, с начальным условием />, а для единственностирешения в этом случае необходимо еще задать температурный режим на концахстержня. Это можно сделать с помощью граничных условий, аналогичных тем,которые были сформулированы для уравнений Пуассона и Лапласа.
Граничное условие Iрода (на конце стержня x=0 заданна температура):/>.
Граничное условие II рода (на конце стержня x=0 задантепловой поток): />.
Граничное условие III рода: />.
Для другого конца стержня x=1 правыечасти граничных условий заменяются соответственно на ψ0(t), ψ1(t), ψ2(t). Заметим, что начальное и граничное условия должны удовлетворятьтак называемым условиям сопряжения, т.е. при условии I рода u0(0)=φ0(0),при условии II рода u0x(0)=φ1(0),при условии III рода -u0x(0)+λu0(0)=φ2(0). Аналогичные условиясопряжения должны выполнятся и на другом конце стержня x=1.
Сформулируем одну из возможных краевых задач. Найтиклассическое решение уравнения />,удовлетворяющее начальному условию /> иследующим граничным условиям />. Этазадача обычно называется первой краевой задачей для уравнения теплопроводности.Соответственно краевые задачи с граничными условиями II роди или III называютсявторой и третьей краевой задачей для уравнения теплопроводности.
Метод конечных разностей (метод сеток)
Численные методы, основанные на разностной аппроксимациипроизводных называется разностным методом, методом конечных разностей илиметодом сеток.
Пусть заданно линейное дифференциальное уравнение,записанное в символическом виде: />. Здесь u – искомое решение уравнения; L – некоторыйдифференциальный оператор, сокращенно обозначающий соответствующуюдифференциальную операцию; f – правая часть уравнения(заданная функция).
Для единственного решения данного уравнения к немунеобходимо присоединить краевые условия: />.
Разностный метод решения этих двух задач можно представить ввиде двух этапов:построение разностной схемы, аппроксимирующей данную непрерывную задачу; получение решения разностной задачи и оценка погрешности этого решения.
Для построения разностной схемы первым шагом является заменаобласти /> непрерывного измененияаргументов областью дискретного их изменения – сеточной областью />, т.е.множеством точек (xn, ym), называемыхузлами сетки. Для квадрата />сеточнуюобласть можно построить следующим образом. Проведем прямые />. Множество точекпересечения этих прямых и составит сеточную область, а сами точки образуют узлысетки. Всякая функция />, определенная нассеке />, называется сеточнойфункцией и обозначается />.
Второй шаг в построении разностной схемы состоит ваппроксимации дифференциального выражения Lu некоторым разностным выражением, а функцию непрерывногоаргумента f – сеточной функцией, т.е. в построениенекоторого разностного аналога для данного уравнения, при данных краевыхусловиях.
Такая аппроксимация приводит к системе алгебраическихуравнений относительно значений некоторой сеточной функции />. Эту систему можно записатьв следующем виде:
/>
Где Lh иφh – разностные операторы,аппроксимирующие соответственно L иl; υh – искомаясеточная функция, аппроксимирующая решение u; fh, φh – заданныесеточные функции, аппроксимирующие f иφ.
Совокупность разносных уравнений, аппроксимирующих исходнуюзадачу – есть разностная схема. Рассмотрим их подробнее на примерах уравнениятеплопроводности и колебания струны.
Разностные схемы для решения уравнениятеплопроводности (параболический тип)
Рассмотрим первую краевую задачу для уравнениятеплопроводности в прямоугольнике />.Требуется найти непрерывное в /> решениезадачи:
/>
В области /> введемпрямоугольную равномерную сетку />{xn, tk} с шагомh=1/N покоординате x и с шагом τ=T/M покоординате t:
/>.
Производные левой части уравнения /> аппроксимируем следующимразностными выражениями: />
В соответствии с данной аппроксимацией построим дваразностных аналога уравнения /> снеизвестной сеточной функцией υhτ:/>
Здесь /> - значениенекоторой сеточной функции fhτ, соответствующейправой части уравнения />. Дляпервой разностной схемы />, а длявторой — /> .
Начальное и граничное условия для первой краевой задачи аппроксимируютсяточно: />
Для второй и третьей краевых задач граничные условияаппроксимируются на основе разностных выражений.
Полагая r=τ/h2 получим /> -для первой разностной схемы, /> - — длявторой разностной схемы.
Анализ показывает, что погрешность аппроксимации схем есть />.
Разностные схемы для решения уравнения колебанияструны (гиперболический тип)
Рассмотрим первую краевую задачу для уравнения колебанияструны в прямоугольнике />.Требуется найти непрерывное в /> решениезадачи:
/>
Применение метода конечных разностей к решению задачи посуществу мало чем отличается от его применения к уравнению теплопроводности. Область/> покрывается сеткой />. Отличие заключается вприближении второй производной по переменной t:
/>.
Разностная аппроксимация принимает вид
/>.
Начальные условия аппроксимируются следующим образом: />.
Граничные условия аппроксимируются точно так же, как и дляуравнения теплопроводности: />.
Значение /> являетсяфиктивным неизвестным, которое можно определить по формуле: />, где γ=τ/h.
Анализ показывает, что погрешность аппроксимации схем есть />.
Список литературыДемидович Б.П., Марон И.А. Основы вычислительной математики. Наука, 1970.Минкова Р.М., Вайсбурд Р.А. Методы вычислительной математики. УПИ, 1981.Боглаев Ю.П. Вычислительная математика и программирование. Высшая школа, 1990.Кацман Ю.Я. Прикладная математика. Численные методы. ТПУ, 2000.