Учебное пособие: Вычислительная математика

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

Pm(x)=a0 + a1x + a2x2+...+amxm.                                                  (4.14)

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

S =

Условия минимума S можно записать, приравнивая нулю частные производные S по  всем переменным a0,  a1,  a2, … , am. Получим систему уравнений

 =  –= 0, или

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

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

a0+ a1+ … +am= , k = 0, 1, … , m     (4.16)

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

ck = , bk = .

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

a0ck + a1ck+1 + … + ck+mam = bk, k = 0, 1, … , m.                                   (4.17)

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

  c0a0 + c1a1 + c2a2… + cmam = b0

  c1a0 + c2a1 + c3a2… + cm+1am  = b1

(4.18)

cma0 + cm+1a1 + cm+2a2… + c2mam  = bm

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

Ca = b.                                                                                          (4.19)

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

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

D = .                                                                (4.20)

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

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

P1(x) = a0 + a1x.

c0 = = n + 1; c1 = = ; c2 = ;                         (4.21)

b0 = = ; b1 = = .                                  (4.22)

           c0           c1                         n+1        

C =                             =                                      ,

           c1               c2                                   

b = (b0, b1)T = (,)T.

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

a0 = , a1 = ,

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

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

a0 = ,  a1 = .                                                    (4.23)

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

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

Шаг 2. Вычислить коэффициенты c0,  c1, b0, b1 по формулам (4.21), (4.22).

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

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

D1 = .                                                          (4.24)

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

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

P2(x) = a0 + a1x + a2x2.

c0 == n+1; c1 ==; c2 =; c3 =; c4 =.  (4.25)

b0 ==; b1 ==; b2 = .                      (4.26)

           c0    c1    c2                 

C =     c1     c2    c3    .

           c2    c3    c4

b = (b0, b1, b2)T .

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

ai = , i = 0, 1,

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

úCú = c0c2c4 + 2c1c2c3  –  c –  сc4  – cc0.                                       (4.27)

               b0    c1  c2 

úC1ú =   b1    c2   c3  = b0c2c4 + b2c1c3 + b1c2c3 – b2c– b1c1c4 – b0c.    (4.28)

              b2    c3   c4

             c0    b0   c2   

úC2ú =  c1    b1    c3 = b1c0c4 + b0c2c3 + b2c1c2 – b1c– b0c1c4 – b2c0c3.   (4.29)

             c2    b2   c4

  c0    c1    b0 

úC3ú =  c1    c2    b1  = b2c0c2 + b1c1c2 + b0c1c3 – b0c– b2c – b1c0c3.    (4.30)

c2    c3    b2

a0 = , a1 = , a2 = .                                                                 (4.31)

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

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

Шаг 2. Вычислить коэффициенты c0,  c1, c2, c3, c4, b0, b1, b2, по формулам (4.25), (4.26).

Шаг 3. Вычислить úCú, úC1ú, úC2ú, úC3ú по формулам (4.27) – (4.30).

Шаг 4. Вычислить a0, a1, a2 по формулам (4.31).

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

D2 = .                                              (4.32)

Шаг 5. Вывести на экран результаты : аппроксимирующую квадратичную функцию P2(x) = a0 + a1x + a2x2  и величину погрешности D2.

Пример 4.6.

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

Таблица 4.1

Вычислим коэффициенты c0,  c1, c2, c3, c4, b0, b1, b2, по формулам (4.25), (4.26):

c0 = 5; c1 = 15; c2 = 55; c3 = 225; c4 = 979;

b0 = 12;  b1 = 53;  b2 = 235.

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

Система уравнений для определения коэффициентов a0 и a1 многочлена первой степени P2(x) = a0 + a1x + a2x2  имеет вид

5a0 + 15a1 = 12

15a0 + 55a1 = 53

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

a0 =  » –2.7,  a1 = » 1.7.    

P1(x) = a0 + a1x =  –2.7 + 1.7x.

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

Система уравнений для определения коэффициентов a0, a1 и a2 многочлена второй степени P2(x) = a0 + a1x + a2x2  имеет вид

5a0 +   15a1 +   55a2 =   12

15a0 +   55a1 + 225a2 =   53

55a0 + 225a1 + 979a2 = 235

По формулам (4.31) найдем  коэффициенты a0, a1 и a2:

a0 » –2.20, a1 » 1.27, a2 » 0.07.

P2(x) = a0 + a1x + a2x2 = –2.20 + 1.27x + 0.07x2.

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

Таблица 4.2

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

D1 =  = 0.245.

D2 =  = 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)

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

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

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

Рис. 5.1

Разобьем отрезок [a, b] на n равных частей длиной h, так, что h = . При этом получим точки a = x0 < x1< x2 < … < xn = b и  xi+1 = xi  + h, i = 0, 1, … , n – 1 (рис. 5.2)

Рис. 5.2

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

Рис. 5.3

Эта фигура состоит из n прямоугольников. Основание i-го прямоугольника образует отрезок [xi, xi+1] длины h, а высота основания равна значению функции в середине отрезка [xi, xi+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пр | £  h2,                                                                          (5.6)

где M2  = |f "(x)|

Пример 5.1.

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

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

Таблица 5.1

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

Iпр  = 0.74713088.

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

f "(x) = (e)" = (4x2  –  2) e.

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

| I  –  Iпр | £ (0.1)2  » 0.84× 10-3.

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

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

Рис. 5.7

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

I=»Iтр =h= (5.7)

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

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

| I  –  Iтр | £  h2,                                                                         (5.8)

где M2  = |f "(x)|.

Пример 5.2.

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

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

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

 I  –  Iтр | £ (0.1)2  » 1.7× 10-3.

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

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

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

y = L2(x) =

f(x) + (x – x) + (x - x)2,    (5.9)

где h = .

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

Ii = »  = ( f(xi) + 4f(x) + f(xi+1)).                     (5.10)

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

 I =» IС = ( f(x0) + f(xn) + 4 + 2).           (5.11)

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

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

| I  –  IС | £ h4,                                                                      (5.12)

где M4 = | f (4)(x)|.

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

I » (f(x0) + f(x2m) + 4 + 2),                                   (5.13)

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

| I  –  IС | £  h4,                                                                    (5.14)

Пример 5.3.

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

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

IС = 0.74682418.

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

f (4)(x) = (16x4 – 48x2 + 12) e, | f (4)(x)| £ 12.

Поэтому

| I  –  IС | £ (0.1)4 » 0.42 × 10-6.

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

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

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

 I  –  Ih  » Chk,                                                                                     (5.15)

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

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

 I  –  Ih/2 » Chk » ( I  –  Ih).                                                              (5.16)

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

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

Ih/2  –  Ih  » Chk(2k – 1).                                                                      (5.17)

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

I  –  Ih/2 »  .                                                                                (5.18)

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

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

I  –  Iпр  » ,                                                                           (5.19)

I  –  Iтр  » ,                                                                            (5.20)

I  –  IС  » .                                                                             (5.21)

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

Пример 5.4.

Найдем значение интеграла   с точностью e = 10-4, используя формулу трапеций и применяя вышеизложенную процедуру дробления шага. В примере 5.2 было получено значение I при h1 = 0.1, Ih =0.74621079. Уменьшим шаг вдвое: h2 = 0.05 и вычислим I= 0.74667084, e2  = ( I- I) = (0.74667084 – 0.74621079) » 1.5×10-4. Так как |e2| > e, то снова дробим шаг: h3 = 0.025, вычисляем I= 0.74678581, e2  = ( I- I) = (0.74678581 – 0.74667084) » 4×10-5. Поскольку  |e3| < e,  требуемая точность достигнута и 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). График решения дифференциального уравнения называется интегральной кривой.

Страницы: 1, 2, 3, 4, 5