Практичне заняття 1. Обрахунок степеневого поліному в системі matlab

Текст MATLAB-програми обрахунку значення степеневого поліному за формулами (1.2) - (1.3) і за MATLAB-функцією polyval такий.

%******************************************************************

%**** Polynom value calculation *******************************

%******************************************************************

clear;

n=10; a=[-1 2 -3 4 -5 6 -7 8 -9 20]; x=0.9;

%**** Calculation on definition *********************************

xn(1)=1; for i=2:n; xn(i)=xn(i-1)*x; end;

pol=a(1); i=2; while i<=n; pol=pol+a(i)*xn(i); i=i+1; end;

pol,

%**** Horner’s formula *****************************************

pol=a(n)*x; i=n-1; while i>1; pol=x*(a(i)+pol); i=i-1; end;

pol=pol+a(1),

%**** MATLAB function ****************************************

a=fliplr(a);

pol=polyval(a,x),

%**** Grafic of polynom ****************************************

x=[-1:0.05:1]; pol=polyval(a,x);

plot(x,pol); title('polynom value plot'); pause;

%**** 3D grafic of polynoms ***********************************

na=40; pol=zeros(na,length(x)); av=pol;

for i=1:na; pol(i,:)=polyval(a,x);

for j=1:length(x); av(i,j)=a(1); end;

a(1)=a(1)-1;

end;

plot3(x,av,pol); grid;

Назва програми і заголовки блоків є коментарями, що починаються символом % .

Перші оператори програми очищують робочу область пам’яті та задають числові значення порядку полінома n, його коефіцієнтів (вектор а) і аргумента х. Всі ці оператори завершуються символом ; , завдяки чому результат їх виконання не виводиться у командне вікно MATLAB.

Наступний блок програми обраховує значення полінома за формулою (1.2). Перший цикл за індексом і від 2 до n з кроком 1 записує у вектор xn цілі степені x від 1 до (n-1). Цикл починається заголовком зі службового слова for і закінчується службовим словом end. У другому циклі обчислено значення полінома pol як суми добутків коефіцієнтів, що містяться у векторі а, на степені аргумента х, накопичені у векторі xn. Цей цикл починається службовим словом while і закінчується також службовим словом end. Числове значення полінома виводиться у командне вікно MATLAB “Command window” оператором pol, який закінчується комою, а не крапкою з комою.

Наступний блок програми обраховує значення полінома за формулою Горнера (1.3). Підрахунок починається з внутрішніх дужок. Числове значення полінома також виводиться у командне вікно MATLAB.

Далі значення полінома обраховується і виводиться втретє за допомогою функції MATLAB polyval. Але для правильного використання цієї функції треба впорядкувати коефіцієнти полінома у векторі а, починаючи з коефіцієнта при найвищому степені аргумента х. Таку перестановку коефіцієтів робить функція MATLAB fliplr.

Графік залежності значення полінома від аргумента будується у наступному блоці програми. Для цього х переозначується як вектор, що містить числа від -1 до 1 з кроком 0.05. Мова MATLAB є матричною. Тому звертання до функції polyval, у якому другий аргумент – вектор х, породжує вектор pol тієї ж розмірності, що й вектор х. Далі функція plot будує і виводить у графічне вікно потрібний графік, який надписується функцією title (дивись рис. 1). Оператор pause зупиняє виконання програми і дозволяє роздивитись графік, поки не буде натиснена будь-яка клавіша на клавіатурі.

Останній блок програми будує тривимірну картинку, на якій зображено na=40 графіків полінома, що змінюються зі зміною першого коефіцієнта з кроком –1. Оскільки порядок коефіцієнтів обернено функцією fliplr, першому коефіцієнту відповідає останній у векторі а на початку програми. Числові значення полінома накопичуються у рядках матриці pol (цикл за індексом і), розмір якої початково задано функцією zeros, де функція length визначає довжину вектора х. Матриця av того ж розміру, що матриця pol. Кожний рядок матриці av містить однакові значення коефіцієнта а(1) (цикл за індексом j). Із зростанням номера рядка значення коефіцієнта а(1) зменшується на 1. Тривимірний графік будує функція plot3. За осями графіка відкладено: значення аргумента х; значення змінного коефіцієнта а(1); відповідні значення полінома (дивись рис. 2). На координатні площини накладено координатні сітки (функція grid).

Лабораторна робота 1. Обрахунок степеневого поліному в системі MATLAB

1.  Набрати та відлагодити текст програми обрахунку значення полінома.

2.  Отримати в командному вікні MATLAB (Command window) три однакові значення полінома.

3.  Отримати двовимірний та тривимірний графіки полінома (рис. 1, рис. 2) та описати їх.

4.  Пояснити всі оператори MATLAB-програми.

Контрольні завдання до розділу 1

1. Назвіть основні критерії вибору метода та алгоритма розв’язування.

2. Назвіть основні ознаки алгоритмічних мов.

3. Назвіть всі етапи розв’язування задачі на ЕОМ.

4. Чому необхідний етап інтерпретації результатів?

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

Розв’язування систем лінійних алгебричних рівнянь є фундаментальною обчислювальною задачею, бо до неї зводиться більшість розрахункових задач.

Системою лінійних алгебричних рівнянь (СЛАР) називають таку систему m рівнянь з n невідомими:

(2.1)

Тут коефіцієнти a_ij; (i=1,…,m; j=1,…,n) та праві частини рівнянь b_i; (i=1,…,m) є заданими, а x_j, (j=1,…,n) – невідомі. Розв’язати СЛАР – значить знайти такі x_j, що задовольняють всі рівняння системи (2.1).

Систему (2.1) зручно записувати у матричному вигляді:

AX=B, (2.2)

де

Кількість рівнянь СЛАР m у загальному випадку не співпадає з кількістю невідомих n. Матриця системи A тоді є прямокутною. Якщо m=n, то матриця A є квадратною, а систему (2.2) називають СЛАР n-го порядку.

Дамо графічну інтерпретацію розв’язування СЛАР 2-го порядку:

(2.3)

Систему рівнянь (2.3) можна зобразити в площині (x₁, x₂) двома прямими (дивись рис. 3). Нехай пряма 1 є геометричним місцем точок, координати яких задовольняють перше рівняння системи, тобто координати кожної точки прямої 1 є розв’язком рівняння 1. Те саме стосується прямої 2 відносно рівняння 2. Тоді координати точки перетину цих прямих (x_1р, x_2р) є розв’язком одночасно обох рівнянь, тобто системи (2.3). Очевидно, якщо точка перетину існує, то вона єдина.

Можливі вироджені випадки, коли прямі паралельні (система (2.3) не має розв’язків), і коли прямі співпадають (розв’язків безмежно багато). У першому випадку виконуються такі співвідношення між коефіцієнтами системи (2.3): а₁₁/а₂₁=а₁₂/а₂₂≠b₁/b₂. У другому – а₁₁/а₂₁=а₁₂/а₂₂=b₁/b₂. Легко переконатись, що в обох випадках визначник матриці А рівний нулю: detA=0. Натомість прямі перетинаються при а₁₁/а₂₁≠а₁₂/а₂₂, а необхідна і достатня умова існування єдиного розв’язку така: detA≠0. Ця умова поширюється на СЛАР будь-якого порядку.

Якщо прямі на рис. 3 майже паралельні, то detA≈0, а розв’язок (x_1р, x_2р) сильно залежить від коефіцієнтів системи (2.3). Ці ознаки свідчать про погану обумовленість системи. Такі системи важко розв’язувати, якщо кількість значущих цифр у числах є обмеженою.

Одним із методів розв’язування квадратних СЛАР є метод Крамера (відомий швейцарський математик 18-го сторіччя). Розв’язок системи (2.3) за цим методом є відношеннями визначників:

x_1р= ; x_2р= ; де D=detA; D₁=det ; D₂=det .

За методом Крамера можна розв’язувати СЛАР будь-якого порядку n. Однак із зростанням n метод стає дуже громіздким, бо кількість арифметичних операцій для обчислення визначника приблизно рівна (n+1)!, а визначників є (n+1).

Можна знаходити розв’язок СЛАР (2.2) при m=n з використанням оберненої матриці А^-1. Домноживши зліва рівняння (2.2) на А^-1, одержуємо Х=А^-1В. Однак за кількістю арифметичних дій обертання матриці не простіше за формули Крамера.

Найпоширенішим обчислювальним методом розв’язування СЛАР є метод, запропонований ще в кінці 18-го сторіччя великим німецьким матема­тиком, фізиком та геодезистом Карлом Фрідріхом Гауссом (Carl Friedrich Gauss). Походячи з майже неграмотної бідної родини, Гаусс досяг слави найвизначнішого математика та обчислювача свого часу. Він заклав наукові основи геодезії та картографії. За кілька років до Морзе він вперше створив електричний телеграф та передав ним повідомлення. На його честь одиниця напруженості магнітного поля названа його ім’ям.

Метод Гаусса грунтується на еквівалентних перетвореннях, які не змінюють розв’язок СЛАР. Розв’язок не зміниться, якщо будь-яке рівняння СЛАР домножити на дійсну константу, а також при додаванні будь-яких рівнянь СЛАР. За методом Гаусса СЛАР еквівалентно перетворюють у систему з верхньою трикутною матрицею (прямий хід), з якої послідовно обчислюють всі складові розв’язку (зворотній хід). Верхньою трикутною є матриця, елементи якої нижче головної діагоналі є нульовими. Метод Гаусса та його модифікації вигідно вирізняються меншою кількістю арифметичних дій, приблизно рівною n³. Разом з тим, обчислення прямого та зворотнього ходу з обмеженою кількістю значущих цифр можуть накопичувати похибки, які спотворюють розв’язок великих погано обумовлених СЛАР.

Катастрофічні похибки метода Гаусса виявились лише після появи ЕОМ, коли стало можливим розв’язувати великі системи. Внаслідок грунтовних досліджень цієї проблеми був створений метод, що майже усунув ці похибки. Цей метод називають QR-розкладом (QR-factorization).

За сучасною матричною термінологією метод Гаусса полягає у розкладі матриці СЛАР на добуток двох матриць, нижньої трикутної L та верхньої трикутної U: A=LU. Завдяки цьому систему AX=B розв’язують у два етапи: спочатку розкладають матрицю А на добуток LU (прямий хід метода Гаусса), а потім послідовно розв’язують системи LY=B та UX=Y (зворотній хід). Метод Гаусса часто називають LU-розкладом (LU-factorization).

Відповідно QR-розклад полягає у розкладі матриці СЛАР на добуток дійсної ортогональної матриці Q і верхньої трикутної R. Поетапний розв’язок системи AX=B виконують так: спочатку обчислюють Y=Q'B, а потім розв’язують систему RX=Y. Верхній штрих означає транспонування матриці.

QR-розклад складніший, ніж LU-розклад, і погано обумовлені системи не так часто зустрічаються. Але треба знати про існування надійнішого метода, коли виникне необхідність.

Перейдемо до розв’язування СЛАР з прямокутними матрицями, коли в системі (2.2) m≠n. В практично цікавих випадках кількість рівнянь СЛАР перевищує кількість невідомих: m>n. Очевидно, тоді не існує розв’язок, що задовольняє всі рівняння, і треба шукати якийсь компроміс.

Розглянемо графічну інтерпретацію розв’язування системи трьох рівнянь з двома невідомими

(2.4)

К ожне з рівнянь зображене на рис. 4 прямою (дивись рис. 3). У загаль­ному випадку три прямі не перетинаються в одній точці, а утворюють три точки перетину (дивись рис. 4). Координати кожної з точок перетину є розв’язком певних двох рівнянь, але не третього. Треба знайти точку S, координати якої не є точним розв’язком всіх рівнянь, а в певному сенсі є компромісом між трьома частковими розв’язками. Напрошується вибрати таку точку, щоб загальна похибка трьох рівнянь системи була найменшою.

Якщо похибка і-го рівняння дорівнює (a_i1x₁+a_i2x₂–b_i)², то наша задача розв’язується за методом найменших квадратів, вперше запропонованим великим Гауссом, а її математичний запис може бути таким:

. (2.5)

Якщо ж похибка і-го рівняння дорівнює |a_i1x₁+a_i2x₂–b_i|, то координати компромісної точки S можна знайти за методом мінімаксу:

. (2.6)

Розв’яжемо систему (2.4) за методом найменших квадратів (2.5). Шуканий мінімум є сумісним розв’язком двох рівнянь, кожне з яких є умовою мінімуму (2.5) за однією зі змінних:

Виконаємо нескладні перетворення:

Ми отримали СЛАР другого порядку, розв’язок якої є розв’язком (2.5). Отже, задача (2.5) зводиться до розв’язку СЛАР другого порядку.

В загальному випадку розв’язок СЛАР з m рівнянь та n невідомих за методом найменших квадратів є задачею

, (2.7) що зводиться до розв’язку СЛАР n-го порядку.

Мінімаксна задача (2.6) в загальному випадку формулюється так:

(2.8) і зводиться до розв’язку задачі лінійного програмування, яка розглядається в розділі 8. Розв’язок (2.8) має деякі переваги порівняно з (2.7). Але задача (2.8) складніша, ніж (2.7), тому СЛАР з прямокутними матрицями найчастіше розв’язують методом найменших квадратів.

Великий Гаусс мав блискучі результати в багатьох галузях знань. Зокрема, він розробив ефективні методи розрахунку орбіт небесних тіл за результатами астрономічних спостережень, де треба було багатократно розв’язувати СЛАР. Саме цими методами користувались молоді математики й астрономи англієць Джон Адамс (Adams) та француз Урбен Левер’є (Le Verrier), коли в середині 19-го сторіччя майже одночасно і незалежно один від одного взялись за розрахунок орбіти невідомої на той час планети Нептун, що мала бути далі від орбіти планети Уран і спотворювала її рух. Адамс завершив дворічні розрахунки раніше за Левер’є, але англійські астро­номи, яким від передав результати, не поспішали підтвердити їх безпосеред­нім спостереженням. Тому гучна честь відкриття „на кінчику пера” нової планети дісталась Левер’є, який згодом став директором Паризької обсерваторії. Адамс зайняв посаду директора Кембріджської обсерваторії і став президентом Лондонського астрономічного товариства. Між бувшими суперниками встановились дружні стосунки, що тривали решту їх життя.