Вычислительная математика

a₀ + a₁ x₀ + a₂x + … + a_nx = y₀

a₀ + a₁ x₁ + a₂x + … + a_nx = y₁

a₀ + a₁ x₂ + a₂x + … + a_nx = y₂ (4.8)

.

a₀ + a₁ x_n + a₂x + … + a_nx = y_n

Вопрос о существовании и единственности интерполяционного многочлена решает следующая теорема:

Теорема 4.1. Существует единственный интерполяционный многочлен степени n, удовлетворяющий условиям (4.6).

Имеются различные формы записи интерполяционного многочлена. Широко распространенной формой записи является многочлен Лагранжа

L_n(x) = = . (4.9)

В частности, для линейной и квадратичной интерполяции по Лагранжу получим следующие интерполяционные многочлены:

L₁(x) = y₀+ y_1,

L₂(x) = y₀+ y₁+ y₂.

Пример 4.3.

Построим интерполяционный многочлен Лагранжа по следующим данным:

x	0	2	3	5
y	1	3	2	5

Степень многочлена Лагранжа для n +1 узла равна n. Для нашего примера многочлен Лагранжа имеет третью степень. В соответствии с (4.9)

L₃(x) = 1+3 + 2 + 5 = 1 + x - x² + x³.

Пример 4.4.

Рассмотрим пример использования интерполяционного многочлена Лагранжа для вычисления значения заданной функции в промежуточной точке. Эта задача возникает, например, когда заданы табличные значения функции с крупным шагом, а требуется составить таблицу значений с маленьким шагом.

Для функции y = sinx известны следующие данные.

x	0	/6	/3	/2
y	0	?		1

Вычислим y(0.25).

Найдем многочлен Лагранжа третьей степени:

L₃(x) = 0 + +

+ 1.

При x = 0.25 получим y(0.25) = sin 0.25 0.249.

Погрешность интерполяции. Пусть интерполяционный многочлен Лагранжа построен для известной функции f(x). Необходимо выяснить, насколько этот многочлен близок к функции в точках отрезка [a, b], отличных от узлов. Погрешность интерполяции равна |f(x) - P_n(x)|. Оценку погрешности можно получить на основании следующей теоремы.

Теорема 4.2. Пусть функция f(x) дифференцируема n +1 раз на отрезке [a, b], содержащем узлы интерполяции x_i [a, b], i = 0, 1, … , n. Тогда для погрешности интерполяции в точке x [a, b] справедлива оценка:

|f(x) - L_n(x)| |_n+1(x)|, (4.10)

где

M_n+1 = |f⁽ⁿ⁺¹⁾(x)|,

_n+1(x) = (x - x₀)(x - x₁)…. (x - x_n).

Для максимальной погрешности интерполяции на всем отрезке [a, b] справедлива оценка:

|f(x) - L_n(x)| |_n(x)| (4.11)

Пример 4.5.

Оценим погрешность приближения функции f(x) = в точке x = 116 и на всем отрезке [a, b], где a = 100, b = 144, с помощью интерполяционного много члена Лагранжа L₂(x) второй степени, построенного с узлами x₀ = 100, x₂ = 144.

Найдем первую, вторую и третью производные функции f(x):

f '(x)= x - ^1/2, f "(x)= - x -^3/2, f'''(x)= x -^5/2.

M₃ = | f'''(x)| = 100 -^5/2 = 10 -⁵.

В соответствии с (4.9) получим оценку погрешности в точке x = 116:

| - L₂(116)| |(116 - 100)(116 - 121)(116 - 144)| = 10 -⁵16528 = 1.410 ^{- 3}.

Оценим погрешность приближения функции f(x) = на всем отрезке в соответствии с (4.11):

| - L₂(x)| |(x - 100)(x - 121)(x -144)| 2.510-³.

4.4 Аппроксимация функций. Метод наименьших квадратов

В инженерной деятельности часто возникает необходимость описать в виде функциональной зависимости связь между величинами, заданными таблично или в виде набора точек с координатами (x_i, y_i), i = 0, 1, 2,... , n, где n - общее количество точек. Как правило, эти табличные данные получены экспериментально и имеют погрешности (рис. 2.5)

Рис.4.2

При аппроксимации желательно получить относительно простую функциональную зависимость (например, многочлен), которая позволила бы "сгладить" экспериментальные погрешности, вычислять значения функции в точках, не содержащихся в исходной таблице.

Эта функциональная зависимость должна с достаточной точностью соответствовать исходной табличной зависимости. В качестве критерия точности чаще всего используют критерий наименьших квадратов, т.е. определяют такую функциональную зависимость f(x), при которой

S =, (4.12)

обращается в минимум.

Погрешность приближения оценивается величиной среднеквадратического уклонения

= . (4.13)

В качестве функциональной зависимости рассмотрим многочлен

P_m(x)=a₀ + a₁x + a₂x²+...+a_mx^m. (4.14)

Формула (4.12) примет вид

S =

Условия минимума S можно записать, приравнивая нулю частные производные S по всем переменным a₀, a₁, a₂, … , a_m. Получим систему уравнений

= -= 0, или

= 0, k = 0, 1, … , m. (4.15)

Систему уравнений (4.15) перепишем в следующем виде:

a₀+ a₁+ … +a_m= , k = 0, 1, … , m (4.16)

Введем обозначения:

c_k = , b_k = .

Система (4.16) может быть записана так:

a₀c_k + a₁c_k+1 + … + c_k+ma_m = b_k, k = 0, 1, … , m. (4.17)

Перепишем систему (4.17) в развернутом виде:

c₀a₀ + c₁a₁ + c₂a₂… + c_ma_m = b₀

c₁a₀ + c₂a₁ + c₃a₂… + c_m+1a_m = b₁

(4.18)

c_ma₀ + c_m+1a₁ + c_m+2a₂… + c_2ma_m = b_m

Матричная запись системы (4.18) имеет следующий вид:

Ca = b. (4.19)

Для определения коэффициентов a_k, k = 0, 1, … , m, и, следовательно, искомого многочлена (4.14) необходимо вычислить суммы c_k, b_k и решить систему уравнений (4.18). Матрица C системы (4.19) называется матрицей Грама и является симметричной и положительно определенной. Эти полезные свойства используются при решении.

Погрешность приближения в соответствии с формулой (4.13) составит

 = . (4.20)

Рассмотрим частные случаи m =1 и m = 2.

1. Линейная аппроксимация (m = 1).

P₁(x) = a₀ + a₁x.

c₀ = = n + 1; c₁ = = ; c₂ = ; (4.21)

b₀ = = ; b₁ = = . (4.22)

c₀ c₁ n+1

C = = ,

c₁ c₂

b = (b₀, b₁)^T = (,)^T.

Решение системы уравнений Ca = b найдем по правилу Крамера:

a₀ = , a₁ = ,

где C - определитель матрицы C, аC_i - определитель матрицы C_i, полученной из матрицы C заменой i-го столбца столбцом свободных членов b, i = 1, 2.

Таким образом,

a₀ = , a₁ = . (4.23)

Алгоритм 4.1 (Алгоритм метода наименьших квадратов. Линейная аппроксимация).

Шаг 1. Ввести исходные данные: x_i, y_i, i=0, 1, 2, ... , n.

Шаг 2. Вычислить коэффициенты c₀, c₁, b₀, b₁ по формулам (4.21), (4.22).

Шаг 3. Вычислить a₀, a₁ по формулам (4.23).

Шаг 4. Вычислить величину погрешности

₁ = . (4.24)

Шаг 5. Вывести на экран результаты: аппроксимирующую линейную функцию P₁(x) = a₀ + a₁x и величину погрешности ₁.

2. Квадратичная аппроксимация (m = 2).

P₂(x) = a₀ + a₁x + a₂x².

c₀ == n+1; c₁ ==; c₂ =; c₃ =; c₄ =. (4.25)

b₀ ==; b₁ ==; b₂ = . (4.26)

c₀ c₁ c₂

C = c₁ c₂ c₃ .

c₂ c₃ c₄

b = (b₀, b₁, b₂)^T .

Решение системы уравнений Ca = b найдем по правилу Крамера:

a_i = , i = 0, 1,

где C - определитель матрицы C, аC_i - определитель матрицы C_i, полученной из матрицы C заменой i-го столбца столбцом свободных членов b.

C = c₀c₂c₄ + 2c₁c₂c₃ - c - сc₄ - cc₀. (4.27)

b₀ c₁ c₂

C₁ = b₁ c₂ c₃ = b₀c₂c₄ + b₂c₁c₃ + b₁c₂c₃ - b₂c- b₁c₁c₄ - b₀c. (4.28)

b₂ c₃ c₄

c₀ b₀ c₂

C₂ = c₁ b₁ c₃ = b₁c₀c₄ + b₀c₂c₃ + b₂c₁c₂ - b₁c- b₀c₁c₄ - b₂c₀c₃. (4.29)

c₂ b₂ c₄

c₀ c₁ b₀

C₃ = c₁ c₂ b₁ = b₂c₀c₂ + b₁c₁c₂ + b₀c₁c₃ - b₀c- b₂c - b₁c₀c₃. (4.30)

c₂ c₃ b₂

a₀ = , a₁ = , a₂ = . (4.31)

Алгоритм 4.2 (Алгоритм метода наименьших квадратов. Квадратичная аппроксимация).

Шаг 1. Ввести исходные данные: x_i, y_i, i=0, 1, 2, ... , n.

Шаг 2. Вычислить коэффициенты c₀, c₁, c₂, c₃, c₄, b₀, b₁, b₂, по формулам (4.25), (4.26).

Шаг 3. Вычислить C, C₁, C₂, C₃ по формулам (4.27) - (4.30).

Шаг 4. Вычислить a₀, a₁, a₂ по формулам (4.31).

Шаг 5. Вычислить величину погрешности

₂ = . (4.32)

Шаг 5. Вывести на экран результаты : аппроксимирующую квадратичную функцию P₂(x) = a₀ + a₁x + a₂x² и величину погрешности ₂.

Пример 4.6.

Построим по методу наименьших квадратов многочлены первой и второй степени и оценим степень приближения. Значения y_i в точках x_i , i =0, 1, 2, 3, 4 приведены в таблице 2.3.

Таблица 4.1

i	0	1	2	3	4
xi	1	2	3	4	5
yi	-1	1	2	4	6

Вычислим коэффициенты c₀, c₁, c₂, c₃, c₄, b₀, b₁, b₂, по формулам (4.25), (4.26):

c₀ = 5; c₁ = 15; c₂ = 55; c₃ = 225; c₄ = 979;

b₀ = 12; b₁ = 53; b₂ = 235.

1. Линейная аппроксимация (m =1).

Система уравнений для определения коэффициентов a₀ и a₁ многочлена первой степени P₂(x) = a₀ + a₁x + a₂x² имеет вид

5a₀ + 15a₁ = 12

15a₀ + 55a₁ = 53

По формулам (4.23) найдем коэффициенты a₀ и a₁:

a₀ = -2.7, a₁ = 1.7.

P₁(x) = a₀ + a₁x = -2.7 + 1.7x.

2. Квадратичная аппроксимация (m =2).

Система уравнений для определения коэффициентов a₀, a₁ и a₂ многочлена второй степени P₂(x) = a₀ + a₁x + a₂x² имеет вид

5a₀ + 15a₁ + 55a₂ = 12

15a₀ + 55a₁ + 225a₂ = 53

55a₀ + 225a₁ + 979a₂ = 235

По формулам (4.31) найдем коэффициенты a₀, a₁ и a₂:

a₀ -2.20, a₁ 1.27, a₂ 0.07.

P₂(x) = a₀ + a₁x + a₂x² = -2.20 + 1.27x + 0.07x².

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

Таблица 4.2

i	0	1	2	3	4
xi	1	2	3	4	5
yi	-1	1	2	4	6
P1(xi)	-1	0.7	2.4	4.1	5.8
P2(xi)	-1	0.62	2.24	4	6.9

Погрешность приближения в соответствии с формулами (4.24) и (4.32) составит

₁ = = 0.245.

₂ = = 0.084.

Тема 5. Численное интегрирование функций одной переменной

5.1 Постановка задачи численного интегрирования

Далеко не все интегралы можно вычислить по известной из математического анализа формуле Ньютона - Лейбница:

I == F(b) - F(a), (5.1)

где F(x) - первообразная функции f(x). Например, в элементарных функциях не выражается интеграл . Но даже в тех случаях, когда удается выразить первообразную функцию F(x) через элементарные функции, она может оказаться очень сложной для вычислений. Кроме того, точное значение интеграла по формуле (5.1) нельзя получить, если функция f(x) задается таблицей. В этих случаях обращаются к методам численного интегрирования.

Суть численного интегрирования заключается в том, что подынтегральную функцию f(x) заменяют другой приближенной функцией, так, чтобы, во-первых, она была близка к f(x) и, во вторых, интеграл от нее легко вычислялся. Например, можно заменить подынтегральную функцию интерполяционным многочленом. Широко используют квадратурные формулы:

, (5.2)

где x_i - некоторые точки на отрезке [a, b],называемые узлами квадратурной формулы, A_i - числовые коэффициенты, называемые весами квадратурной формулы, n 0 - целое число.

5.2 Метод прямоугольников

Формулу прямоугольников можно получить из геометрической интерпретации интеграла. Будем интерпретировать интеграл как площадь криволинейной трапеции, ограниченной графиком функции y = f(x), осью абсцисс и прямыми x = a и x = b (рис. 5.1).

Рис. 5.1

Разобьем отрезок [a, b] на n равных частей длиной h, так, что h = . При этом получим точки a = x₀ < x₁< x₂ < … < x_n = b и x_i+1 = x_i + h, i = 0, 1, … , n - 1 (рис. 5.2)

Рис. 5.2

Заменим приближенно площадь криволинейной трапеции площадью ступенчатой фигуры, изображенной на рис. 5.3.

Рис. 5.3

Эта фигура состоит из n прямоугольников. Основание i-го прямоугольника образует отрезок [x_i, x_i+1] длины h, а высота основания равна значению функции в середине отрезка [x_i, x_i+1], т е. f(рис. 5.4).

Рис. 5.4

Тогда получим квадратурную формулу средних прямоугольников:

I = I_пр = (5.3)

Формулу (5.3) называют также формулой средних прямоугольников. Иногда используют формулы

I I = , (5.4)

I I = , (5.5)

которые называют соответственно квадратурными формулами левых и правых прямоугольников.

Геометрические иллюстрации этих формул приведены на рис. 5.5 и 5.6.

Рис. 5.5

Рис. 5. 6

Оценка погрешности. Для оценки погрешности формулы прямоугольников воспользуемся следующей теоремой .

Теорема 5.1. Пусть функция f дважды непрерывно дифференцируема на отрезке [a, b]. Тогда для формулы прямоугольников справедлива следующая оценка погрешности:

| I - I_пр | h², (5.6)

где M₂ = |f "(x)|

Пример 5.1.

Вычислим значение интеграла по формуле средних прямоугольников (5.3) с шагом h = 0.1.

Составим таблицу значений функции e(табл. 5.1):

Таблица 5.1

x	e	x	e
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50	1.0000000 0.9975031 0.9900498 0.9777512 0.9607894 0.9394131 0.9139312 0.8847059 0.8521438 0.8166865 0.7788008	0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00	0.7389685 0.6976763 0.6554063 0.6126264 0.5697828 0.5272924 0.4855369 0.4448581 0.4055545 0.3678794

Производя вычисления по формуле (5.3), получим:

I_пр = 0.74713088.

Оценим погрешность полученного значения. Имеем:

f "(x) = (e)" = (4x² - 2) e.

Нетрудно убедиться, что | f "(x)| M₂ = 2. Поэтому по формуле(5.4)

| I - I_пр | (0.1)² 0.84 10^-3.

5.3 Метод трапеций

Выведем формулу трапеций так же, как и формулу прямоугольников, из геометрических соображений. Заменим график функции y = f(x) (рис.5.1) ломаной линией (рис.5.7), полученной следующим образом. Из точек a = x₀, x₁, x₂,…, x_n = b проведем ординаты до пересечения с кривой y = f(x). Концы ординат соединим прямолинейными отрезками.

Рис. 5.7

Тогда площадь криволинейной трапеции приближенно можно считать равной площади фигуры, составленной из трапеций. Так как площадь трапеции, построенной на отрезке [x_i, x_i+1] длины h = , равна h , то, пользуясь этой формулой для i = 0, 2, … , n - 1, получим квадратурную формулу трапеций:

I=I_тр =h= (5.7)

Оценка погрешности. Для оценки погрешности формулы трапеций воспользуемся следующей теоремой.

Теорема 5.2. Пусть функция f дважды непрерывно дифференцируема на отрезке [a, b]. Тогда для формулы трапеций справедлива следующая оценка погрешности:

| I - I_тр | h², (5.8)

где M₂ = |f "(x)|.

Пример 5.2.

Вычислим значение интеграла по формуле трапеций (5.7) и сравним полученный результат с результатом примера 5.1.

Используя таблицу значений функции eиз примера 5.1 и производя вычисления по формуле трапеций (5.7), получим: I_тр = 0.74621079.

Оценим погрешность полученного значения. В примере (5.1) получили оценку: | f "(x)| M₂ = 2. Поэтому по формуле (5.8)

I - I_тр | (0.1)² 1.7 10^-3.

Сравнивая результаты примеров 5.1 и 5.2, видим, что метод средних прямоугольников имеет меньшую погрешность, т.е. он более точный.

5.4 Метод Симпсона (метод парабол)

Заменим график функции y = f(x) на отрезке [x_i, x_i+1], i = 0, 2, … , n - 1, параболой, проведенной через точки (x_i, f(x_i)), (x,f(x)), (x_i+1, f(x_i+1)), где x - середина отрезка [x_i, x_i+1]. Эта парабола есть интерполяционный многочлен второй степени L₂(x) с узлами x_i, x, x_i+1. Нетрудно убедиться, что уравнение этой параболы имеет вид:

y = L₂(x) =

f(x) + (x - x) + (x - x)², (5.9)

где h = .

Проинтегрировав функцию (5.9) на отрезке [x_i, x_i+1], получим

I_i = = ( f(x_i) + 4f(x) + f(x_i+1)). (5.10)

Суммируя выражение (5.10) по i = 0, 1, 2, … , n - 1, получим квадратурную формулу Симпсона (или формулу парабол):

I = I_С = ( f(x₀) + f(x_n) + 4 + 2). (5.11)

Оценка погрешности. Для оценки погрешности формулы Симпсона воспользуемся следующей теоремой.

Теорема 5.2. Пусть функция f имеет на отрезке [a, b] непрерывную производную четвертого порядка f ⁽⁴⁾(x). Тогда для формулы Симпсона (5.9) справедлива следующая оценка погрешности:

| I - I_С | h⁴, (5.12)

где M₄ = | f ⁽⁴⁾(x)|.

Замечание. Если число элементарных отрезков, на которые делится отрезок [a, b], четно , т.е. n = 2m, то параболы можно проводить через узлы с целыми индексами, и вместо элементарного отрезка [x_i, x_i+1] длины h рассматривать отрезок [x_2i, x_2i+2] длины 2h. Тогда формула Симпсона примет вид:

I (f(x₀) + f(x_2m) + 4 + 2), (5.13)

а вместо оценки (5.10) будет справедлива следующая оценка погрешности:

| I - I_С | h⁴, (5.14)

Пример 5.3.

Вычислим значение интеграла по формуле Симпсона (5.11) и сравним полученный результат с результатами примеров 5.1 и 5.2.

Используя таблицу значений функции eиз примера 5.1 и производя вычисления по формуле Симпсона (5.11) , получим:

I_С = 0.74682418.

Оценим погрешность полученного значения. Вычислим четвертую производную f ⁽⁴⁾(x).

f ⁽⁴⁾(x) = (16x⁴ - 48x² + 12) e, | f ⁽⁴⁾(x)| 12.

Поэтому

| I - I_С | (0.1)⁴ 0.42 10^-6.

Сравнивая результаты примеров 5.1, 5.2 и 5.3, видим , что метод Симпсона имеет меньшую погрешность, чем метод средних прямоугольников и метод трапеций.

5.5 Правило Рунге практической оценки погрешности

Оценки погрешности по формулам (5.4), (5.8) и (5.12) являются априорными. Они зависят от длины элементарного отрезка h, и при достаточно малом h справедливо приближенное равенство:

I - I^h Ch^k, (5.15)

где I^h приближенное значение интеграла, вычисленное по одной из формул (5.3), (5.5), (5.9), C 0 и k > 0 - величины, не зависящие от h.

Если уменьшить шаг h в два раза, то, в соответствии с (5.15) получим:

I - I^h/2 Ch^k ( I - I^h). (5.16)

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

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

I^h/2 - I^h Ch^k(2^k - 1). (5.17)

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

I - I^h/2 . (5.18)

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

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

I - I_пр , (5.19)

I - I_тр , (5.20)

I - I_С . (5.21)

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

Пример 5.4.

Найдем значение интеграла с точностью = 10^-4, используя формулу трапеций и применяя вышеизложенную процедуру дробления шага. В примере 5.2 было получено значение I при h₁ = 0.1, I^h =0.74621079. Уменьшим шаг вдвое: h₂ = 0.05 и вычислим I= 0.74667084, ₂ = ( I- I) = (0.74667084 - 0.74621079) 1.510^-4. Так как |₂| > , то снова дробим шаг: h₃ = 0.025, вычисляем I= 0.74678581, ₂ = ( I- I) = (0.74678581 - 0.74667084) 410^-5. Поскольку |₃| < , требуемая точность достигнута и 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

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

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

y(t₀ ) = y₀, (6.2)

где t₀ - некоторое заданное значение аргумента t, а y₀ - начальное значение функции.

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

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

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

| f(t, y₁) - f(t, y₂)| L| y₁ - y₂|,

где L некоторая постоянная, а y₁ , y₂ - произвольные значения.

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

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

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

Сходимость численных методов решения задачи Коши. Пусть y(t) - решение задачи Коши. Назовем глобальной погрешностью (или просто погрешностью) численного метода функцию _i = y(t_i) - y_i , заданную в узлах сетки t_i. В качестве абсолютной погрешности примем величину R = | y(t_i) - y_i|

Численный метод решения задачи Коши называется сходящимся, если для него R 0 при h 0. Говорят, что метод имеет p-ый порядок точности, если для погрешности справедлива оценка R Ch^p, p > 0, C - константа, C 0.

6.2 Метод Эйлера

Простейшим методом решения задачи Коши является метод Эйлера.

Будем решать задачу Коши

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

y(t₀ ) = y₀,

на отрезке [t₀, T]. Выберем шаг h = , и построим сетку с системой узлов

t_i = t₀ + ih, i = 0, 1, …, n.

В методе Эйлера вычисляются приближенные значения функции y(t) в узлах сетки :y_i y(t_i).

Заменив производную y' (t) конечными разностями на отрезках [t_i, t_i+1], i = 0, 1, …, n - 1, получим приближенное равенство:

= f(t_i, y_i), i = 0, 1, …, n - 1,

которое можно переписать так:

y_i+1 = y_i + h f(t_i, y_i), i = 0, 1, …, n - 1. (6.3)

Формулы (6.3) и начальное условие (6.2) являются расчетными формулами метода Эйлера.

Геометрическая интерпретация одного шага метода Эйлера заключается в том, что решение на отрезке [t_i, t_i+1] заменяется касательной y = y' (t_i)( t - t_i), проведенной в точке (t_i, y(t_i)) к интегральной кривой, проходящей через эту точку. После выполнения n шагов неизвестная интегральная кривая заменяется ломаной линией (ломаной Эйлера).

Оценка погрешности. Для оценки погрешности метода Эйлера воспользуемся следующей теоремой.

Теорема 6.2. Пусть функция f удовлетворяет условиям:

K, = L. (6.4)

Тогда для метода Эйлера справедлива следующая оценка погрешности:

R = | y(t_i) - y_i| = ,

где l - длина отрезка [t₀, T]. Мы видим , что метод Эйлера имеет первый порядок точности.

Оценка погрешности метода Эйлера часто бывает затруднительна, так как требует вычисления производных функции f(t, y(t)). Грубую оценку погрешности дает правило Рунге (правило двойного пересчета), которое используется для различных одношаговых методов, имеющих p-ый порядок точности. Правило Рунге заключается в следующем. Пусть y - приближения, полученные с шагом , а y - приближения, полученные с шагом h. Тогда справедливо приближенное равенство:

|y- y(t_i)| |y- y| . (6.5)

Таким образом, чтобы оценить погрешность одношагового метода с шагом , нужно найти то же решение с шагом h и вычислить величину, стоящую справа в формуле (6.5), т е.

R |y- y| (6.6)

Так как метод Эйлера имеет первый порядок точности, т. е. p = 1, то приближенное равенство (6.6) примет вид

R |y- y| (6.7)

Используя правило Рунге, можно построить процедуру приближенного вычисления решения задачи Коши с заданной точностью . Нужно, начав вычисления с некоторого значения шага h, последовательно уменьшать это значение в два раза, каждый раз вычисляя приближенное значение y, i = 0, 1, …, n. Вычисления прекращаются тогда, когда будет выполнено условие:

R |y- y| < . (6.8)

Для метода Эйлера условие (6.8) примет вид

R |y- y| < (6.9)

Приближенным решением будут значения y, i = 0, 1, …, n.

Пример 6.1.

Найдем решение на отрезке [0, 1] следующей задачи Коши:

y' (t) = y - , (6.10)

y(0) = 1.

Возьмем шаг h = 0.2. Тогда n = = 5.

В соответствии с (6.3) получим расчетную формулу метода Эйлера:

y_i+1 = y_i + 0.2 , y₀ = 1, i = 0, 1, 2, 3, 4, 5.

Решение представим в виде таблицы 6.1:

Таблица 6.1

i	0	1	2	3	4	5
ti	0	0.2	0.4	0.6	0.8	1.0
yi	1.0000	1.2000	1.3733	1.5294	1. 6786	1.8237

Уравнение (6.10) есть уравнение Бернулли. Его решение можно найти в явном виде:

y = . (6.11)

Для сравнения точного и приближенного решений представим точное решение (6.11) в виде таблицы 6.2:

Таблица 6.2

i	0	1	2	3	4	5
ti	0	0.2	0.4	0.6	0.8	1.0
y(ti)	1.0000	1.1832	1.3416	1.4832	1. 6124	1.7320

Из таблицы видно, что погрешность составляет R = | y(t_i) - y_i| = 0.0917.

6.3 Модифицированные методы Эйлера

Первый модифицированный метод Эйлера. Суть этого метода состоит в следующем. Сначала вычисляются вспомогательные значения искомой функции y в точках t = t_i + с помощью формулы:

y = y_i + f_i = y_i +f(t_i, y_i).

Затем находится значение правой части уравнения (6.1) в средней точке

f = f(t, y)

и затем полагается

y_i+1 = y_i + h f, i = 0, 1, …, n - 1. (6.12)

Формулы (6.12) являются расчетными формулами первого модифицированного метода Эйлера.

Первый модифицированный метод Эйлера является одношаговым методом со вторым порядком точности

Второй модифицированный метод Эйлера - Коши. Суть этого метода состоит в следующем. Сначала вычисляются вспомогательные значения

= y_i + h f(t_i, y_i). (6.13)

Затем приближения искомого решения находятся по формуле:

y_i+1 = y_i + [f(t_i, y_i) + f(t_i+1, )], i = 0, 1, …, n - 1. (6.14)

Формулы (6.14) являются расчетными формулами второго модифицированного метода Эйлера - Коши.

Второй модифицированный метод Эйлера - Коши, так же, как и первый, является одношаговым методом со вторым порядком точности.

Оценка погрешности. Приближенная оценка погрешности модифицированных методов Эйлера осуществляется как и для простого метода Эйлера с использованием правила Рунге (см. предыдущий раздел 6.2). Так как оба модифицированных метода Эйлера имеют второй порядок точности, т. е. p = 2, то оценка погрешности (6.6) примет вид

R |y- y|. (6.15)

Используя правило Рунге, можно построить процедуру приближенного вычисления решения задачи Коши модифицированными методами Эйлера с заданной точностью . Нужно, начав вычисления с некоторого значения шага h, последовательно уменьшать это значение в два раза, каждый раз вычисляя приближенное значение y, i = 0, 1, …, n. Вычисления прекращаются тогда, когда будет выполнено условие:

R |y- y| < . (6.16)

Приближенным решением будут значения y, i = 0, 1, …, n.

Пример 6.2.

Применим первый модифицированный метод Эйлера для решения задачи Коши

y' (t) = y - , y(0) = 1,

рассмотренной ранее в примере 6.1.

Возьмем шаг h = 0.2. Тогда n = = 5.

В соответствии с (6.3) получим расчетную формулу первого модифицированного метода Эйлера:

y_i+1 = y_i + h f = y_i + 0.2 f, где

f = f(t, y) = y - ,

t = t_i + = t_i + 0.1,

y = y_i +f(t_i, y_i) = y_i +0.1,

t₀ = 0, y₀ = 1, i = 0, 1, …, 4.

Решение представим в виде таблицы 6.3:

Таблица 6.3

i	ti	yi	f(ti, yi)	t	y	h f
0 1 2 3 4 5	0 0.2 0.4 0.6 0.8 1.0	1 1.1836 1.3426 1.4850 1.6152 1.7362	0.1 0.0850 0.0747 0.0677 0.0625	0.1 0.3 0.5 0.7 0.9	1.1 1.2682 1.4173 1.5527 1.6777	0.1836 0.1590 0.1424 0.1302 0.1210

Третий столбец таблицы 6.3 содержит приближенное решение y_i, i = 0, 1, …, 5.

Сравним полученное приближенное решение с точным решением (6.11), представленном в таблице 6.2. Виднм, что погрешность составляет R = | y(t_i) - y_i| = 0.0042.

Пример 6.3.

Применим второй модифицированный метод Эйлера - Коши для решения задачи Коши

y' (t) = y - , y(0) = 1,

рассмотренной ранее в примерах 6.1 и 6.2. Так же, как и ранее, зададим шаг h = 0.2. Тогда n = = 5.

В соответствии с (6.14) получим расчетную формулу метода Эйлера - Коши:

y_i+1 = y_i + [f(t_i, y_i) + f(t_i+1, )] = y_i + 0.1[f(t_i, y_i) + f(t_i+1, )],

где

f(t_i, y_i) = y_i -

= y_i + h f(t_i, y_i) = y_i + 0.1

t₀ = 0, y₀ = 1, i = 0, 1, …, 4.

Решение представим в виде таблицы 6.4:

Таблица 6.4

i	ti	yi	f(ti, yi)	ti+1

Подобные документы

• Численные методы решения типовых математических задач

  Решение систем линейных алгебраических уравнений методом простой итерации. Полиномиальная интерполяция функции методом Ньютона с разделенными разностями. Среднеквадратическое приближение функции. Численное интегрирование функций методом Гаусса.
  
  курсовая работа [2,4 M], добавлен 14.04.2009
  

• Решение систем линейных алгебраических уравнений

  Метод Зейделя как модификация метода простой итерации. Особенности решения систем линейных алгебраических уравнений. Анализ способов построения графика функций. Основное назначение формул Симпсона. Характеристика модифицированного метода Эйлера.
  
  контрольная работа [191,3 K], добавлен 30.01.2014
  

• MatLab в задачах вычислительной математики

  Решение систем линейных алгебраических уравнений методом исключения Гаусса. Табулирование и аппроксимация функций. Численное решение обыкновенных дифференциальных уравнений. Приближенное вычисление определенных интегралов. Решение оптимизационных задач.
  
  курсовая работа [1,6 M], добавлен 21.11.2013
  

• Решение систем линейных алгебраических уравнений

  Задачи вычислительной линейной алгебры. Математическое моделирование разнообразных процессов. Решение систем линейных алгебраических уравнений большой размерности. Метод обратной матрицы и метод Гаусса. Критерии совместности и определенности системы.
  
  курсовая работа [220,0 K], добавлен 21.10.2011
  

• Решение систем уравнений различными способами

  Формирование системы их пяти уравнений по заданным параметрам, ее решение методом Гаусса с выбором главного элемента. Интерполяционный многочлен Ньютона. Численное интегрирование. Решение нелинейных уравнений. Метод Рунге-Кутта четвертого порядка.
  
  контрольная работа [115,5 K], добавлен 27.05.2013
  

• Метод Ньютона (метод касательных). Решение систем нелинейных алгебраических уравнений

  Геометрическая интерпретация методов Ньютона, итерации и спуска. Определение корня уравнения с заданной степенью точности. Решение систем нелинейных алгебраических уравнений. Нахождение эквивалентного преобразования для выполнения условия сходимости.
  
  курсовая работа [371,6 K], добавлен 14.01.2015
  

• Численные методы

  Решение системы линейных уравнений с неизвестными методами Гаусса, Зейделя и простой итерации. Вычисление корня уравнения методами дихотомии, хорды и простой итерации. Нахождение приближённого значения интеграла с точностью до 0,001 методом Симпсона.
  
  контрольная работа [1,7 M], добавлен 05.07.2014
  

• Численные методы

  Приближенные числа и действия над ними. Решение систем линейных алгебраических уравнений. Интерполирование и экстраполирование функций. Численное решение обыкновенных дифференциальных уравнений. Отделение корня уравнения. Поиск погрешности результата.
  
  контрольная работа [604,7 K], добавлен 18.10.2012
  

• Системы линейных уравнений и неравенств

  Основные понятия теории систем уравнений. Метод Гаусса — метод последовательного исключения переменных. Формулы Крамера. Решение систем линейных уравнений методом обратной матрицы. Теорема Кронекер–Капелли. Совместность систем однородных уравнений.
  
  лекция [24,2 K], добавлен 14.12.2010
  

• Численные методы анализа

  Численные методы решения систем линейных уравнений: Гаусса, простой итерации, Зейделя. Методы аппроксимации и интерполяции функций: неопределенных коэффициентов, наименьших квадратов. Решения нелинейных уравнений и вычисление определенных интегралов.
  
  курсовая работа [322,7 K], добавлен 27.04.2011
  

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