Файл: n1.pdf

Даже для простых дифференциальных уравнений первого порядка не всегда удается получить аналитическое решение. Поэтому большое значение имеют численные методы решения. Численные методы позволяют определить приближенные значения искомого решения y(t) на некоторой выбранной сетке значений аргумента ti, (i = 0, 1, …). Точки ti называются узлами сетки, а величина hi = ti+1 – ti – шагом сетки. Часто рассматривают равномерные сетки, для которых шаг hi постоянен,

hi = h = T −n t0 . При этом решение получается в ви-

де таблицы, в которой каждому узлу сетки ti соответствуют приближенные значения функции y(t) в

узлах сетки yi ≈ y(ti).

Численные методы не позволяют найти решение в общем виде, зато они применимы к широкому классу дифференциальных уравнений.

Сходимость численных методов решения задачи Коши

Пусть y(t) – решение задачи Коши. Назовем

глобальной погрешностью (или просто погрешно-

стью) численного метода функцию εi= y(ti) – yi, заданную в узлах сетки ti. В качестве абсолютной

погрешности примем величину R = max | y(ti) – yi|

0≤i≤n

Численный метод решения задачи Коши называется сходящимся, если для него R → 0 при

104

Вычитая из второго, третьего и четвертого уравнений системы (3.10) первое уравнение, умноженное соответственно на m12 , m13 , m14 , получим новую

систему:

Вычитая из третьего и четвертого уравнений системы (3.11) второе уравнение, умноженное соответственно на m 32 и m 24 , приходим к системе:

3-ий шаг. Вычислим множитель:

41

Вычитая из четвертого уравнения системы (3.12) третье, умноженное на m 34 , приведем систему к треугольному виду:

Обратный ход. Из последнего уравнения системы (3.13) находим x4 = 1.000. Подставляя значение x4 в третье уравнение, получим x3 = 2.000. Подставляя найденные значения x4 и x3 во второе уравнение, найдем x2 = 3.000. Наконец, из первого уравнения, подставив в него найденные значения x4, x3 и x2, вычислим x1 = – 1.000.

Итак система (3.10) имеет следующее реше-

ние:

x1 = 1.000, x2 = 2.000, x3 = 3.000, x4 = – 1.000.

3.3. Метод исключения Гаусса с выбором главного элемента по столбцу

Хотя метод Гаусса является точным методом, ошибки округления могут привести к существенным погрешностям результата. Кроме того, исключение по формулам (3.7) нельзя проводить, если элемент главной диагонали a kkk−1 равен нулю.

42

Производную y'(t) в каждой точке (t, y) можно геометрически интерпретировать как тангенс угла α наклона касательной к графику решения, проходящего через эту точку, т е.: k = tgα = f (t, y).

Уравнение (6.1) определяет целое семейство решений. Чтобы выделить одно решение, задают

начальное условие:

где t0 – некоторое заданное значение аргумента t, а y0 – начальное значение функции.

Задача Коши заключается в отыскании функции y = y(t), удовлетворяющей уравнению (6.1) и начальному условию (6.2). Обычно определяют решение задачи Коши на отрезке, расположенном справа от начального значения t0, т. е. для

t [t0, T].

Разрешимость задачи Коши определяет следующая теорема.

Теорема 6.1. Пусть функция f (t, y) определена и непрерывна при t0 ≤ t ≤ T, –∞ < y < ∞ и удовлетворяет условию Липшица:

| f (t, y1) – f (t, y2)| ≤ L| y1 – y2|,

где L некоторая постоянная, а y1 , y2 – произвольные значения.

Тогда для каждого начального значения y0 существует единственное решение y(t) задачи Коши для t [t0, T].

103

= 13 (0.74678581 – 0.74667084) ≈ 4 10–5.

Поскольку |ε3| < ε, требуемая точность достигнута и I ≈ 0.7468 ± 0.0001.

Тема 6. Численное решение дифференциальных уравнений

6.1. Постановка задачи Коши

Известно, что обыкновенное дифференциальное уравнение первого порядка имеет вид:

y' (t) = f(t, y(t)). (6.1)

Решением уравнения (6.1) является дифференцируемая функция y(t), которая при подстановке в уравнение (6.1) обращает его в тождество. На рис. 6.1 приведен график решения дифференциального уравнения (6.1). График решения дифференциального уравнения называется интегральной кривой.

Рис. 6.1.

102

Если элемент a kkk−1 мал, то велики ошибки округле-

ния при делении на этот элемент. Для уменьшения ошибок округления применяют метод исключения Гаусса с выбором главного элемента по столбцу.

Прямой ход так же, как и для схемы единственного деления, состоит из n – 1 шагов. На первом шаге прежде, чем исключать переменную x1, уравнения переставляются так, чтобы в левом верхнем углу был наибольший по модулю коэффициент ai1, i = 1, 2, …, n. В дальнейшем, на k-м шаге, прежде чем исключать переменную xk, уравнения переставляются так, чтобы в левом верхнем углу был наибольший по модулю коэффициент aik, i = k, k +1, …, n. После этой перестановки исключение переменной xk производят, как в схеме единственного деления.

Трудоемкость метода. Дополнительные действия по выбору главных элементов требуют примерно n2 операций, что практически не влияет на общую трудоемкость метода.

Пример 3.2.

Применим метод исключения Гаусса с выбором главного элемента по столбцу для решения системы уравнений (3.10) из примера 3.1.

Прямой ход. 1-ый шаг. Так как коэффициент a11 = 2.0 наибольший из коэффициентов первого столбца, перестановки строк не требуется и 1-ый

43

шаг полностью совпадает с 1-ым шагом примера 3.1. Из второго, третьего и четвертого уравнений исключается переменная x1 и система приводится к виду (3.11).

2-ой шаг. Наибольший по модулю коэффициент при x2 в системе (3.11)

a132 = – 1.15.

Поэтому переставим уравнения следующим образом:

Вычислим множители:

Вычитая из третьего и четвертого уравнений системы (3.14) второе уравнение, умноженное соответственно на m 32 и m 24 , приходим к системе:

4.28478x3 – 7.38261x4 = 20.23696

2.28522x3 – 2.81739x4 = 9.67305

44

Используя правило Рунге, можно построить процедуру приближенного вычисления интеграла с заданной точностью ε. Нужно, начав вычисления с некоторого значения шага h, последовательно уменьшать это значения в два раза, каждый раз вычисляя приближенное значение I hi . Вычисления прекращаются тогда, когда результаты двух последующих вычислений будут различаться меньше, чем на ε.

Пример 5.4.

1

Найдем значение интеграла ∫e−x2 dx с точно-

0

стью ε = 10-4, используя формулу трапеций и применяя вышеизложенную процедуру дробления шага. В примере 5.2 было получено значение I h1 при

h1 = 0.1, Ih =0.74621079.

Уменьшим шаг вдвое: h2 = 0.05 и вычислим

= 13 (0.74667084 – 0.74621079) ≈ 1.5 10– 4.

Так как |ε2| > ε, то снова дробим шаг: h3 = 0.025, вычисляем

I h3 = 0.74678581, ε2 = 13 ( I h3 – I h2 ) = 101

Непосредственное использование оценок погрешности (5.4), (5.8) и (5.12) неудобно, так как при этом требуется вычисление производных функции f (x). В вычислительной практике используются другие оценки.

Вычтем из равенства (5.15) равенство (5.16):

Учитывая приближенное равенство (5.16), получим следующее приближенное равенство:

Приближенное равенство (5.18) дает апостериорную оценку погрешности. Вычисление этой оценки называется правилом Рунге. Правило Рунге

– это эмпирический способ оценки погрешности, основанный на сравнении результатов вычислений , проводимых с разными шагами h.

Для формул прямоугольников и трапеций k = 2, а для формулы Симпсона k = 4. Поэтому для этих формул приближенное равенство (5.18) принимает вид:

3-ий шаг. Вычислим множитель:

Вычитая из четвертого уравнения системы (3.15) третье, умноженное на m 34 , приведем систему к треугольному виду:

Обратный ход. Обратный ход полностью совпадает с обратным ходом примера 3.1. Решение системы имеет вид:

x1 = 1.000, x2 = 2.000, x3 = 3.000, x4 = – 1.000.

3.4. Вычисление определителя методом исключения Гаусса

Из курса линейной алгебры известно, что определитель треугольной матрицы равен произведению диагональных элементов. В результате метода исключений Гаусса система линейных уравнений (3.2) с квадратной матрицей A приводится к эквивалентной ей системе (3.8) с треугольной матрицей An. Поэтому

det A = (–1)s det An,