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

МЕТОДИ РОЗВ’ЯЗАННЯ СИСТЕМ ЛІНІЙНИХ АЛГЕБРАЇЧНИХ РІВНЯНЬ.

Розглянемо чисельні методи розв’язання систем лінійних алгебраїчних рівнянь

Ax=f T (1)

де A- матриця m*m, x=(x1,x2,...,xm)- шуканийвектор,

f=(f1, f2, ... , fm) -заданий вектор.

Припускаємо, що та визначник матриці А відмінний від нуля, так що існує єдиний розв’язок х. З курсу алгебри відомо, що систему (1) можна розв’язати за формулами Крамера*. Для великих m цей спосіб практично нереалізований тому, що потребує порядку m! aрифметичних дій. Тому широко використовуються інші методи розв’язання, наприклад, метод Гаусса**, який потребує дій.

Методи чисельного розв’язання системи (1) поділяються на дві групи:

-прямі методи;

-ітераційні методи.

У прямих (або точних) методах розв’язок x системи (1) відшукується за скінченну кількість арифметичних дій. Внаслідок похибок заокруглення прямі методи насправді не приводять до точного розв’язку системи (1) і назвати їх точними можливо лише залишаючи осторонь похибки заокруглення.

Ітераційні методи (їх також називають методами послідовних наближень)полягають у тому, що розв’язок x системи (1) відшукується як границя при послідовних наближень де n- номер ітерації. Як правило, за скінченну кількість ітерацій ця границя не досягається.

______________________

* Крамер Габрієль (1704-1752)- швейцарський математик.

** Гаус Карл Фридрих (1777-1855)- німецький математик, астроном, фізик, геодезист, професор Гетінгенського університету.

МЕТОД ГАУССА .

Запишемо систему (1) у розгорнутому вигляді:

а11x1+a12x2+...+a1mxm=f1 ,

a21x1+a22x2+...+a2mxm =f2 , (2)

......................................

am1x1+am2x2+...+ammxm =fm .

Метод Гаусса розв’язання системи (2) полягає у послідовному вилученні невідомих x1, x2, ..., xm-1 з цієї системи.

Припустимо, що a110 . Поділив перше рівняння на a11, одержимо

x1+c12x2 +...+c1m xm =y1 , (3)

де : c1j=a1j /a11 ; j=2,m ; y1=f1/a11 .

Розглянемо тепер рівняння системи (2), що залишилися

ai1x1+ai2x2+...+aimxm=fi ; i= 2,m . (4)

Помножимо (3) на ai1 та віднімемо одержане рівняння з і-го рівняння системи (4), i=2,m.

У результаті одержимо наступну систему рівнянь:

x1+c12x2+...+c1jxj+...+c1mxm =y1 ,

(1) (1) (1) (1)

a22x2+... +a2jxj+...+a2mxm=f2 ,

............................................ (5)

(1) (1) (1) (1)

am2x2+...+amjxj+...+ammxm=fm .

Tут позначено:

(1) (1)

aij=aij-c1jai1; fi=fi -y1ai1; i,j=2,m . (6)

Матриця системи (5) має вигляд:

Матриці такої стуктури заведено позначати так:

де хрестиками позначені ненульові елементи.

У системі (5) невідоме х міститься тільки в першому рівнянні, тому у подальшому достатньо мати справу із скороченою системою рівнянь:

(1) (1) (1) (1)

a22x2 +...+a2jxj +...+a2mxm =f2 ,

.............................................. (7)

(1) (1) (1) (1)

am2x2 +...+amjxj +...+ammxm =fm .

Тим самим ми здійснили перший крок методу Гаусса . Коли , то з системи (7) зовсім аналогічно можна вилучити невідоме x2 і прийти до системи, еквівалентній (2),що має матрицю такої структури:

При цьому перше рівняння системи (5) залишається без зміни.

Вилучая таким же чином невідомі х 3, х4 ,... ,xm-1,приходимо остаточно до системи рівнянь виду:

x1 +c12x2 +...+c1,m-1xm-1+c1mxm =y1,

x2 +...+c2,m-1xm-1+c2mxm =y2 ,

................................

xm-1+cm-1,mxm=ym-1,

xm=ym ,

що еквівалентна початковій системі (2) .

Матриця цієї системи

містить нулі усюди нижче головної діагоналі. Матриці такого виду називаються верхніми трикутними матрицями. Нижньою трикутною матрицею називається така матриця, у якої дорівнюють нулю усі елементи, що містяться вище головної діагоналі.

Побудова системи (8) складає прямий хід методу Гаусса. Зворотнiй хід полягає у відшуканні невідомих х1, ... ,хmз системи (8). Тому що матриця системи має трикутний вигляд, можна послідовно, починаючи з хm, відшукати всі невідомі. Дійсно, xm=ym,

x m-1 =ym-1 -cm-1,m x m i т. д.

Загальні форми зворотнього ходу мають вигляд:

xi =yi - е cijxj ; i=m-1,1; xm =ym . (10)

j=i+1

При реалізації на ЕОМ прямого ходу методу Гаусса немає необхідності діяти із змінними x1 ,x2 ,... ,xm. Досить вказати алгоритм,за яким початкова матриця А перетворюється до трикутного вигляду (9), та вказати відповідне перетворення правих частин системи.

Одержимо ці загальні формули.

Нехай вже здійснені перші к-1 кроків, тобто вже вилучені змінні

x1 , x2,..., xk-1 .

Тоді маємо систему:

x1+c12 x2 +...+c1,k-1xk-1+ c1kxk+....+c1mxm =y1 ,

x2 +...+c2,k-1xk-1+ c2kxk+....+c2mxm =y2 ,

..............................................

xk-1+ck-1,kxk+...+ck-1,mxm=yk-1 , (11)

(k-1) (k-1) (k-1)

akkxk+...+akmxm =fk ,

............................

(k-1) (k-1) (k-1)

amkxk+...+ammxm =fm .

Розглянемо К-те рівняння цієї системи:

(k-1) (k-1) (k-1)

akkxk+...+akmxm=fk ,

та припустимо, що . Поділивши обидві частини цього рівняння на , одержимо

xk+ck,k+1xk+1+...+ckmxm=yk , (12)

(k-1) (k-1) (k-1) (k-1)

де ckj=akj / akk ; j=k+1,m ; yk=fk / akk .

Далі помножимо рівняння (12) на та віднімемо одержане співвідношення з i-го рівняння системи (11). У результаті остання група рівнянь системи (11) набуває наступного вигляду:

x k+ck,k+1xk+1 +...+ ckmxm=yk,

(k) (k) (k)

ak+1,k+1xk+1+...+ ak+1,mxm=fk+1,

.......................................

(k) (k) (k)

am,k+1xk+1+... + ammxm=fm ,

(k) (k-1) (k-1) (k) (k-1) (k-1)

де: aij =aij - aikckj ; i,j=k+1,m ; fi= fi - aikyk ; i=k+1,m .

Таким чином, у прямому ході методу Гаусса коефіцієнти рівнянь перетворюються за наступним правилом

(0)

akj =akj ; k,j=1,m ;

(k-1) (k-1)

ckj=akj /akk ; j=k+1,m ; k=1,m ; (13)

(k) (k-1) (k-1)

aij =aij - aikckj ; i,j=k+1,m ; k=1,m . (14)

Обчислення правих частин системи (8) здійснюється за формулами:

(0) (k-1) (k-1)

fk=fk ; yk = fk / akk ; k=1,m ; (15)

(k) (k-1) (k-1)

fi = fi - aikyk ; k=1,m . (16)

Коефіціенти cij і праві частини yi ; i=1,m ; j=i+1,m зберігаються у пам’яті ЕОМ і використовуються при здійсненні зворотнього ходу за формулами (10).

Основним обмеженням методу є припущення, що усі елементи , на які здійснюється ділення, відрізняються від нуля. Число називається провідним елементом на К-му кроці вилучення. Навіть, якщо деякий провідний елемент не дорівнює нулеві, а просто є близьким до нуля, в процесі обчислень може мати місце нагромадження похибок. Вихід з цієї ситуації полягає в тому, шо як провідний елемент вибирається не , а інше число ( тобто на К-му кроці вилучається не xk, а інша змінна xj , ) . Така стратегія вибору провідних елементів здійснюється в методі Гаусса з вибором головного елементу, який буде розглянуто пізніш.

А тепер підрахуємо число арифметичних дій, що необхідні для розв’язання системи за допомогою методу Гаусса. Оскільки виконання операцій множення і ділення на ЕОМ потребує набагато більше часу, ніж виконання додавання і віднімання, обмежимось підрахуванням числа множень і ділень.

1.Обчислення коефіцієнтів за формулами (13) потребує:

(m-k)=1+2+...+(m-1)= ділень .

2.Обчислення усіх коефіцієнтів за формулами (14) потребує

множень (тут ми використовуємо легко перевіряєму за індукцією рівність ).

Таким чином, обчислення ненульових елементів трикутної матриці С потребує

операцій множення і ділення.

3.Обчислення правих частин yk за формулами (15) потребує m

ділень, а відшукання за формулами (16)

множень. Разом, обчислення правих частин перетвореної системи (8) потребує

дій множення і ділення.

Усього для реалізації прямого ходу методу Гаусса необхідно виконати

дій.

4.Для реалізації зворотнього ходу методу Гаусса за формулами (10) необхідно

множень.

Всього, для реалізації методу Гаусса необхідно виконати

дій множення і ділення, причому основний час витрачається на прямий хід. Для великих mчисло дій множення і ділення у методі Гаусса близьке до За витратами часу та необхідній машинній пам’яті метод Гаусса придатний для розв’язання систем рівнянь (2) загального вигляду з кількістю змінних m порядку 100.

ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ .

Пусть имеется функция которую необходимо продифференцировать несколько раз и найти эту производную в некоторой точке.

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

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

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

Будем предполагать, что функция задана в равностоящих узлах

Ее значения и значения производных в узлах будем обозначать

Пусть функция задана в двух точках и ее значения

Посстроим интерполяционный многочлен первой степени

Производная равна

Производную функцию в точке приближенно заменяем производной интерполяционного многочлена

(1)

Величина называется первой разностной производной.

Пусть задана в трех точках

Интерполяционный многочлен Ньютона второй степени имеет вид

Берем производную

В точке она равна

Получаем приближенную формулу

(2)

Величина называется центральной разностной производной.

Наконец, если взять вторую производную

получаем приближенную формулу.

(3)

Величина называется второй разностной производной.

Формулы (1)-(3) называются формулами численного дифференцирования.

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

В дальнейшем нам понадобится следующая лемма.

Лемма 1. Пусть произвольные точки, Тогда существует такая точка что

Доказательство. Очевидно неравенство

По теореме Больцано-Коши о промежуточных значениях непрерывной функции на замкнутом отрезке она принимает все значения между и Значит существует такая точка что выполняет указанное в лемме равенство.

Погрешности формул численного дифференцирования дает следующая лемма.

Лемма 2.

1.Предположим, что Тогда существует такая точка , что

(4)

Если то существует такая точка , что

(5)

Когда то существует такая, что

(6) Доказательство. По формуле Тейлора

откуда следует (4).

Если то по формуле Тейлора

(7)

где

Подставим (7) в Получаем

Заменяя в соответствии с леммою 1

получаем

Откуда и следует (6).

Равенство (5) доказывается аналогично ( доказательство провести самостоятельно).

Формулы (4)-(6) называются формулами численного дифференцирования с остаточными членами.

Погрешности формул (1)-(3) оцениваются с помощью следующих неравенств, которые вытекают из соотношений (4)-(6):

Говорят, что погрешность формулы (1) имеет первый порядок относительно (или порядка ), а погрешность формул (2) и (3) имеет второй порядок относительно (или порядка ). Также говорят, что формула численного дифференцирования (1) первого порядка точности (относительно ), а формулы (2) и (3) имеют второй порядок точности.

Указанным способом можно получать формулы численного дифференцирования для более старших производных и для большего количества узлов интерполирования.

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

(8)

Пусть в некоторой окрестности точки производные, через которые выражаются остаточные члены в формулах (5), (6), непрерывны и удовлетворяют неравенствам

(9)

где - некоторые числа. Тогда полная погрешность формул (2), (3) (без учета погрешностей округления) в соответствии с (5), (6), (8), (9)не превосходит соответственно величин

Минимизация по этих величин приводит к следующим значениям :

(12)

при этом

(13)

Если при выбранном для какой-либо из формул (2), (3) значении отрезок не выходит за пределы окрестности точки , в которой выполняется соответствующее неравенство (9), то найденное есть оптимальным и полная погрешность численного дифференцирования оценивается соответствующей величиной (13).

ИНТЕРПОЛИРОВАНИЕ СПЛАЙНАМИ.

Интерполирование многочленом Лагранжа или Ньютона на отрезке с использованием большого числа узлов интерполяции часто приводит к плохому приближению, что объясняется сильным накоплением погрешностей в процессе вычислений. Кроме того из-за расходимости процесса интерполяции увеличение числа узлов не обязано приводить к повышению точности. Для того, чтобы избежать больших погрешностей, весь отрезок разбивают на частичные отрезки и на каждом из частичных отрезков приближенно заменяют функцию многочленом невысокой степени ( так называемая кусочно-полиномиальная интерполяция).

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

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

Преимущество сплайнов перед обычной интерполяцией является, во-первых, их сходимость, и, во-вторых, устойчивость процесса вычислений.

Рассмотрим частный, но распространенный в вычислительной практике случай, когда сплайн определяется с помощью многочленов третьей степени ( кубический сплайн).

Пусть на задана непрерывная функция. Введем узлы ( сетку):

и обозначим

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

а) на кождом сегменте функция является многочленом третьей степени;

б) функция , а так же ее первая и вторая производные непрерывны на ;

в)

Последнее условие называется условием интерполирования.

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

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

(1)

где - коэффициенты, подлежащие определению. Выясним смысл введенных коэффициентов. Имеем

поэтому

Из условий интерполирования получаем, что

Доопределим , кроме того , .

Далее , требование непрерывности функции приводит к условиям

Отсюда,учитывая выражения для функций получаем при уравнения

Обозначая перепишем эти уравнения в виде

(2)

Условия непрерывности первой производной

приводят к уравнениям

(3)

Из условий непрерывности второй производной получаем уравнения

. (4)

Объединяя (2) -(4) , получим систему уравнений относительно неизвестных

Два недостающих условия получают, задавая те или иные граничные условия для Предположим, например, что функция удовлетворяет условиям Тогда естественно требовать, чтобы Отсюда получаем

т.е.

Заметим, что условие совпадает с уравнением (4) при . Таким образом, приходим к замкнутой системе уравнений для определения коэффициентов кубического сплайна:

Убедимся в том, что эта система имеет единственное решение. Исключим из (5)- (7) переменные и получим систему, содержащую только Для этого рассмотрим два соседних уравнения (7) :

и вычтем второе уравнение из первого. Тогда получим

Подставляя найденное выражение для в правую часть уравнения (6), получим

(8)

Далее, из уравнения (5) получаем

И подставляя эти выражения в (8) , приходим к уравнению

Окончательно для определения коэффициентов получаем систему уравнений

(9)

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

(10)

Таким образом, доказано, что существует единственный кубический сплайн, определяемый условиями а)-в) и граничными условиями Заметим , что можно рассматривать и другие граничные условия.

ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ.

На практике редко удается вычислить точно определенный интеграл. Например, в элементарных функциях не вычисляется функция Лапласа

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

Рассмотрим некотрые широко используемые приемы приближенного вычисления определенных интегралов.

Квадратурные формулы.

Введем понятие квадратурные формулы. Пусть дан определенный интеграл

(1)

от непрерывной на отрезке функции . Приближенное неравенство

(2)

где - некоторые числа, - некотрые точки отрезка , называется квадратурной формулой, определяемой весами и узлами .

Говорят, что квадратурная формула точна для многочленов степени , если при замене на произвольный алгебраический многочлен степени приближенное равенство (2) становится точным.

Рассмотрим наиболее простые квадратурные формулы.

Формула прямоугольников. Допустим, что . Положим приближенно

(3)

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

Найдем остаточный член , т.е. погрешность формулы (3) .

Пусть

(4)

Так как

то согласно формуле Тейлора с остаточным членом в форме Лагранжа имеем

(5)

где -некоторые точки ,

Функция является первообразной для Поэтому для интеграла, стоящего в левой части приближенного равенства (3), из формулы Ньтона - Лейбница с расчетом (5) вытекает следующее соотношеие

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

(6)

Формула трапеций. Пусть Полагаем

(7)

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

Найдем остаточный член, т.е. погрешность формулы (7). Выразим і где - функция (4), по формуле Тейлора с остаточным членом в интегральной форме :

(8)

(9)

Согласно (8) имеем

(10)

Отделив в правой части (9) слагаемое и заменив его выражением (10), с учетом того, что находим

Преобразуем теперь второе слагаемое в правой части, используя обобщенную теорему о среднем.

* Формула Тейлора с остаточным членом в интегральной форме

Теорема 1 (обобщенная теорема о среднем). Пусть причем на Тогда существует такая точка что

Доказательство. Положим

(11)

Тогд, так как то

и, следовательно,

Если то и в качестве можн взять любую точку из

Если то вытекает существование такого числа с, удовлетворяющего неравенствам ( для этого делим все части на ):

(12)

что

(13)

По теореме о промежуточных значениях непрерывной функции в силу (11) , (12) найдется точка , в которой что вместе с равенством (13) доказывает теорему .

Теперь, так как то по доказанной теоремою

где - некоторая точка . Подставляя полученное в , приходим к формуле трапеций с остаточным членом :

(14)

Формула Симпсона . Предположим, что Интеграл приближенного заменяем площадью заштрихованной криволинейной трапеции, ограниченной сверху параболой, проходящей через точки де

Указанная парабола задается уравнением

в чем нетрудно убедиться, положив поочередно (ее можно также получить, построив интерполяционный многочлен второй степени и приводя подобные ) Отсюда находи ( проверить самостоятельно)

Таким образом , формула Симпсона , называемая также формулой парабол , имеет вид

(15)

Положим где -функция (4). Поскольку

то согласно формул Тейлора с остаточным членом в интегральной форме имеем

Отсюда получаем

(16)

т.к. остальные члены взаимно уничтожаются.

Поскольку то применяя к интегралу (16) теорему 1 , а затем к полученному результату лемму, находим

(17)

где нектрые точки.

Принимая во внимание, что из (16), (17) приходим к формуле

(18) т.е. к формуле Симпсона с остаточным членом.

Рассмотрим квадратурные формулы прямоугольников (3), трапеций (7) и Симпсона (15) называются каноничными.

Усложненные квадратурные формулы.

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