Matlab 2

[ виправити ] текст може містити помилки, будь ласка перевіряйте перш ніж використовувати.

скачати

Міністерство освіти Республіки Білорусь

Установа освіти

«Гомельський державний університет ім. Ф. Скорини »

Математичний факультет

Кафедра МПУ

MATLAB

Реферат

Виконавець:

Студентка групи М-53

Гумарева Л.С.

Гомель 2004

Введення

MATLAB - матрична лабораторія - найбільш розвинена система програмування для науково-технічних розрахунків, доповнена до теперішнього часу кількома десятками більш приватних застосувань, що відносяться до обчислювальної математики, обралботке інформації, конструювання електронних приладів, економіці та ряду інших розділів прикладної науки. Вивчення MATLAB'а за фірмовою документації, яка тепер додається на інсталяційному компакт-диску, займає у починаючих користувачів занадто багато часу не тільки із-за необхідності читати її англійською мовою зі специфічним сленгом, але, головним чином, з огляду на неминучого для таких настанов послідовного і досить формального викладення великого обсягу інформації, а наявні російською мовою допомоги в основному йдуть цьому стереотипу. Навіть для досвідченого фахівця з розрахунками на комп'ютері таке вивчення пов'язане з невиправдано великими витратами праці.

MATLAB призначений насамперед для програмування чисельних алгоритмів. Він розробляється вже більше 15 років і виник на основі більш ранніх прикладних пакетів LINPACK і EIGPACK, створених в 1970-і рр.. в США, і в свою чергу вплинув на появу таких систем, як MathCad, MAPLE і Mathematica. Удосконалення системи MATLAB відбувалося як у зв'язку з досягненнями в обчислювальній математиці, так і у зв'язку із змінами в архітектурі персональних комп'ютерів і розвитком загальносистемних засобів. З часом MATLAB був доповнений цілим рядом вже згадуваних додатків (toolboxes), далеко раздвинувший межі його застосовності. Далі мова піде лише про ядро MATLAB'а, яке ми будемо називати системою, і конкретно про її версії 5.2, випущеної фірмою MathWorks в січні 1998 р.

MATLAB - система програмування високого рівня, що працює як інтерпретатор і включає великий набір інструкцій (команд) для виконання найрізноманітніших обчислень, завдання структур даних та графічного представлення інформації. Команди ці розбиті на тематичні групи, розташовані в різних директоріях системи. Тепер у системі близько 800 команд, і приблизно половина з них цілком доступна починаючому користувачеві. Команди з великим можливим обсягом обчислень написані на С, але багато і таких команд, які представлені в термінах цих перших. Тому система виявляється майже відкритою для користувача. Є великі можливості для виведення двовимірної і тривимірної графіки і засоби управління нею. Користувач може без особливих труднощів додавати свої команди і писати програми в термінах вже існуючих команд; дещо складніше робити це в рамках Фортрану і С. Можна обмінюватися даними з програмами на цих мовах, а з них звертатися до системи. Стислість і наочність програмування і виняткові можливості візуалізації результатів роблять систему дуже ефективною при пошуках і апробації нових алгоритмів, при проведенні разових розрахунків і в навчальному процесі, оскільки її можна освоювати без попереднього знайомства з основами програмування й виконувати такі складні приклади, які неможливо робити з використанням інших систем.

Документація по системі і її додатків містить багато тисяч сторінок, і тому природно постає питання про те, як її освоювати. Робота з системою вимагає певної математичної підготовки, так що навчання можна починати на другому курсі вузу. Основні відомості про систему викладені у посібнику / 1 / - / 2 /: / 1 / - це підручник з описом обчислювальних можливостей і архітектури системи, / 2 / - опис її графічних можливостей. Звичайно, можна читати підряд / 1 /, / 2 / і при необхідності звертатися за уточненнями до команди help або довідником / 3 /, в якому описані майже всі команди. Але набагато більш ефективним, на наш погляд, є виклад основних обчислювальних процедур за допомогою найбільш уживаних команд системи. Саме так ми і познайомимося з MATLAB'ом, а точніше приблизно з 30-40 його командами. Після цих занять користуватися документацією / 1 / і / 2 / буде набагато легше.

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

Для роботи з системою досить мати комп'ютер PC 486 з оперативною пам'яттю хоча б 16 Mb і з встановленими на ньому системами Windows 95 і MATLAB 5.2. У дійсності MATLAB може працювати і з друогімі операційними системами, такими, наприклад, як Macintosh, Unix і OS / 2.

За кордоном вийшло вже досить багато навчальних посібників за системою, але на російську мову жодне з них поки не переводилося і навіть в центральних бібліотеках їх тепер немає через скорочення фінансування. Видані у нас допомоги (наприклад, / 4 / - / 12 /) в основному йдуть посібникам / 1 / - / 3 /, тоді як нам представляється корисним дати менш формальне введення в предмет, спираючись насамперед на інтуїцію слухача.

1. Змінні

Змінні можуть бути числовими, текстовими та інших типів. У нас будуть тільки числові (це у всіх деталях) і текстові (зовсім небагато). Назва змінної починається з латинської літери, далі можуть бути букви і цифри (не більше 31 символу). Малі та великі літери тут різняться.

1. Числові змінні. Це числа, вектори, матриці та багатовимірні масиви. У комп'ютері всі числа представлені приблизно з 16 десятковими знаками, під кожне дійсне число відводиться 8 байтів, під комплексний - 16.

1.1 Введення чисел

Цілі числа. У системі вони не виділяються явно. Наберемо і виконаємо окремо кожну команду:

a = 2 a = 2.0 a = 2; a = 1:6 b = 1:20 c = 10: -2:5

Командна вікно. Командний рядок. Редагування командного рядка. Буфер виконаних команд. Як вибирати інформацію з командного вікна і з буфера виконаних командних рядків. Не можна допускати збігу імені змінної з ім'ям якої-небудь команди.

Дійсні числа. Виконаємо окремо наступні команди:

d = 0.5:0.3:2.5 d =. 5: .3:2.5 d =. 5 +1: .3 -. 1:2.5 * 2 ​​length (d)

d (end) d (end-2) d (1) d (0) d (2:7) d (7: -1:2) d (150)

f = linspace (1.5,30,143); length (f)

Індекси завжди починаються зі значення 1. Команди набираються на малому латинською регістрі. Можлива многопараметрічность команд.

Діапазон дійсних чисел:

realmax realmin

Інші константи MATLAB'а:

pi i j eps

Їх не слід псувати.

Комплексні числа:

q = 1 +2 * iq = 1 +2 i real (q) imag (q) abs (q) conj (q) s = angle (q) (тут - p <s <= p).

q = 1 +2 * i; r = 3; fi = 0: .01: pi; z = q + r * exp (i * fi); plot (z) Це верхня півколо.

1.2 Введення векторів

Вектори-рядки:

a = 1:6 linspace (1,6,10)

Вектори-стовпці:

a = (1:6) 'linspace (1,6,10)'

Оператори. 'І':

y1 = linspace (1,6,4) '; y2 = y1; y = y1 + i * y2; y.' y '

Команди linspace і: застосовні для завдання тільки речових векторів.

Введення матриць. A (i, j) - елемент з i-го рядка і j-го стовпця. A (k) - k-й елемент таблиці, витягнутої в стовпець.

A = [1,2; 3,4] A = [1; 2,3; 4] A (2,2) A (3) A (5) size (A) A (3,4) = 10 size ( A)

A (5) = 6 size (A) A (22) = 3 A = A (:) A (22) = 3 size (A) [m, n] = size (A)

A = reshape (1:24,4,6) size (A) A ([1, end],:) = [] A (:, [1, end]) = [] size (A)

1. 3 Деякі спеціальні матриці

m = 3; n = 4; eye (m, n) eye (m) eye (n) ones (m, n) ones (m) ones (n) zeros (m, n)

rand (m, n) rand (m, n) rand ('state', 0) rand (m, n) rand (m) Це рівномірний розподіл на інтервалі (0, 1).

randn (m, n) randn ('state', 100) Це нормальний розподіл, у нього мат.ожіданіе = 0, дисперсія = 1

v1 = 1:4 v2 = 7:12 toeplitz (v1, v2) toeplitz (v1)

1. 4 Деякі прості команди

A = reshape (1:24,4,6) triu (A) triu (A, 0) triu (A, 2) triu (A, -1) tril (A)

v = 1:5 diag (v) diag (v, 2) diag (v, -1)

diag (A) diag (A, 2) diag (A, -1)

A = reshape (1:24,4,6) rot90 (A) rot90 (A, 2)

Видачі на екран. Команда format з різними можливостями.

У звичайному форматі (forrmat short) видається 5 знаків, для цілих чисел 9 знаків, порядки змінюються від -308 до +308. У повному форматі (format long e) 16 знаків.

a = 2 a =. 001 a = 1e-3 a = 1e-5 a = 123456789 a = 1234567891 a = 1 +3 * i

format long e, 2 ^ .5, format short

Опція format short e дозволяє отримувати рівні стовпці.

Вони беруться в лапки (на букві е. латинською регістрі), символ займає 2 байти. Використовуються для завдання заголовків в числових ведучих і на графіках, для завдання формул і т.д. Можна перекладати текстові змінні в числові і навпаки. Виконаємо в командному рядку

t = 'Moscow - столиця Росії' (після дефіса потрібно перейти на російський шрифт і потім не забути знову повернутися на латинський).

Інші типи змінних - комірки і структури.

Система help.

help видає список директорій системи;

help <ім'я директорії> видає список команд директорії;

help <ім'я команди> видає опис команди.

type <ім'я команди> видає текст команди або програми користувача, якщо він складений у термінах MATLAB 'а.

2. Елементи xy-графіки

1.Як відкривати графічне вікно:

figure whitebg zoom on

Тепер побудуємо графік функціпі y = sin (2 p x), 0 <= x <= 5, виконавши рядок

x = 0:1 e-3: 5; y = sin (2 * pi * x); plot (y) plot (x, y), grid

Використання режиму zoom:

k = 100; y = sin (2 * pi * k * x); plot (y)

2.Автоматіческое чергування квітів. Тепер будемо, як правило, нумерувати рядка.

1; x = linspace (0,1,20); k =. 1: .1: .8; y = k '* x; plot (x, y)

Тут визначається вектор-рядок x = 0:20, потім вектор-рядок k з 8 кутових коефіцієнтів, далі виходить матриця y = k '* x як добуток вектора-стовпця k' на вектор-рядок x. Рядки цієї матриці складаються з точок відповідних прямолінійних відрізків. Нарешті, будуються графіки цих відрізків як функцій від x - перша нижня лінія (вона жовта) відповідає k =. 1, остання, теж жовта, - для k =. 8. Ми бачимо, що кольори, яких всього 7, чергуються циклічно в такому порядку (під російськими англійські назви):

жовтий фіолетовий блакитний червоний зелений синій білий

yellow magenta cyan red green blue white

Викличемо рядок 1 і відредагуємо в ній команду plot:

1; x = linspace (0,1,20); k =. 1: .1: .8; y = k '* x; plot (x, y,' g. ')

тобто додамо там третій (текстової, бо він в апострофа) аргумент. Всі криві на малюнку стануть зеленими (green), а лінії будуть зображуватися окремими точками. Аналогічно вживаються і інші кольори з цього списку - за першою літерою. У текстовому аргументі може бути до трьох символів. Для зображення точок графіка крім. вживаються ще: - -. * Xo + і деякі інші символи.

3.Графікі в полярних координатах:

x = 1: .01:3; nx = length (x); r = x. ^ 2; fi = linspace (0,5 * pi, nx); polar (fi, r)

4.Еще один приклад - легко будуються багатозначні функції:

x = 0: .1:6 * pi; y = cos (x); plot (x, y) plot (y, x)

5. Управління осями:

axis off axis on axis ([-10,10, -5,20]) axis auto axis equal axis square

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

3. Прості приклади, що ілюструють ефективність MATLAB 'а

1. Підсумовування. Знайдемо при заданому n часткову суму ряду s (n) = 1 / k ^ 2, k = 1: n. Для цього виконаємо рядок

1; n = 100; k = 1: n; f = k. ^ (-2); plot (cumsum (f)), [sum (f), pi ^ 2 / 6] = 1000

Команда cumsum (f) підраховує всі часткові суми s (k) від f (1: k) для кожного k від 1 до n, так що на графіку можна спостерігати процес формування потрібної нам величини. У кінці рядка видається чисельний і точний результати:

ans = 1.6350 1.6449.

Вважаючи n = 1000, отримаємо

ans = 1.6439 1.6449,

тобто помилку в 1 одиницю 4-й значущої цифри.

Збіжність не завжди настільки очевидна, як на цьому графіку. Щоб у цьому переконатися, ускладнимо наш приклад: при заданих m> 1 і n знайдемо часткову суму ряду s (m, n) = sum (1 / k ^ m), k = 1: n (при m = 1 виходить вже розходиться гармонійний ряд). Для проведення обчислень відредагуємо рядок 1:

2; m = 2; n = 1000; k = 1: n; f = k. ^ (- M); plot (cumsum (f)), sum (f)

= 1.5 = 1 e 4

= 1.2

і спочатку для перевірки отримаємо свій старий результат. Але вже при m = 1.5 у нас, дивлячись на графік, немає повної впевненості у досягненні збіжності. Це тим більше так при m = 1.2: для n = 1000 ans = 4.3358, а для n = 1e4 ans = 4.7991. Факт збіжності ряду при m = 1.01 не можна встановити чисельно з-за низької швидкості його збіжності.

Щоб краще запам'ятати дію команди cumsum, обчислимо ò (x / sin (x)) dx, x Î [0, 3]. Підінтегральна функція f = x / sin (x) не має в нулі особливості, і тому досить виконати рядок

3; n = 100; h = 3 / n; x = h / 2: h :3-h / 2; f = x. / sin (x); plot (h * cumsum (f)), grid, sum ( h * f) = 1000

тобто апроксимувати f в серединах інтервалів (ці точки x називають напівцілим на відміну від кінців рахункових інтервалів - цілих точок). Порівняння відповіді ans = 8.4495 і графіка наводить на думку про те, що поки збіжність ще не досягнуто, але при n = 1e3 ans = 8.4552, так що при n = 1e2 зі збіжністю в дійсності все гаразд, а зростання функції h * cumsum ( f) на правому кінці відбувається через зростання там функції f - це можна побачити, виконавши

4; plot (f)

Для матриці A команди sum і cumsum працюють вздовж стовпців (значить, по першому індексу), а для вектора - вздовж нього незалежно від того, рядок це або стовпець. Щоб підсумувати для матриці A уздовж її рядків, потрібно виконати sum (A, 2), тобто вказати для виконання команди другої індекс. Це правило відноситься до багатьох командам MATLAB'a і до багатовимірних матриць теж - за умовчанням мається на увазі перший індекс, а в іншому випадку потрібно завжди вказувати, за яким індексом повинна працювати команда, і це вказівка ​​не зберігається для наступних команд.

2. Твори. Аналогічно підсумовування за допомогою команд prod і cumprod обчислюються і обробляються твори. Наприклад, знайдемо Õ (1-1 / k ^ 2), k = 2:100 (при k ® ¥ Õ ® 1 / 2), виконавши рядок

1; n = 100; k 2 = (2: n). # 2; a = 1-1. / K 2; cp = cumprod (a); cp (end), plot (cp / .5), grid

Результат cp (end) = 0.5050 говорить про те, що збіжність тут не дуже швидка. Це видно і з графіка, на якому представлена ​​відносна помилка результату. Зверніть увагу на назви змінних k2 = k ^ 2 і cp = cumprod (..): при виборі імен змінних завжди потрібно прагнути до того, щоб ці імена хоч якось відображали суть справи (це особливо важливо при написанні великих програм, де багато змінних).

При обчисленні творів можна вийти за числову шкалу. Знайдемо, наприклад, для яких k можна знайти k!. Ясно, що максимально допустимий km навряд чи більше 200, так що рядок

2; n = 200; k = 1: n; kf = cumprod (k); plot (kf)

повинна дати відповідь на наше питання. Через швидке зростання kf і обмеженої розв'язності дисплея (це не більше 0.5% від максимального значення на графіку) ми бачимо лише одну точку kf (km), перед якою, як нам помилково здається, йдуть нулі і за якою йдуть числа inf (infinity ), взагалі ніяк не представлені на малюнку. Точно так само графіка обходиться і змінною NaN (not a number), і ця обставина може бути іноді корисним. Змінна NaN виникає в таких ситуаціях:

0 / 0 inf - inf inf / inf

Змінні inf і NaN (вони виходять зі знаком) можна використовувати в програмах. Для визначення km виконаємо рядок

3; sum (isinf (kf))

в якій isinf (kf) видасть 1 на тих позиціях вектора розмірів kf, де елементи kf є inf, і 0 на інших позиціях. Оскільки ans = 30, km = n-30 = 170, що можна було б отримати і відразу, виконавши рядок

4; km = sum (isfinite (kf))

де isfinite зазначає ті елементи числової змінної, які відмінні від inf і NaN. При виході твору за числову шкалу для співмножників можна використовувати команди

log (взяття натурального логарифма),

log10 (взяття десяткового логарифма),

abs (взяття модуля),

sign (взяття знака, видає 1, 0 і -1).

3. Логічні завдання. Зазвичай при освоєнні програмування логічні дії даються важче арифметичних. Наведемо тут два простих прикладу завдань логічного характеру.

1. Напишемо рядок для знаходження спільних елементів двох векторів:

x = 1:20; y = 15:30; [X, Y] = meshgrid (x, y); v = X (X == Y)

2. Другий приклад трохи складніше, і початківці вивчати MATLAB зазвичай намагаються вирішити його за допомогою циклів for-end, що абсолютно неправильно. Взявши на сторонах одиничного квадрата по 200 інтервалів, визначимо, скільки точок вийшла таким чином сітки потрапляє всередину вписаною в нього кола. Потрібна програма має вид

1; tic, x = 0:1 / 200:1; [X, Y] = meshgrid (x); M = abs (X + i * Y-.5-i *. 5) <1 / 2; s = sum (M (:)), t1 = toc

і дасть відповідь s = 31397 точок, t1 = 0.16 сек, тоді як рядок для циклів for-end

2; tic, s = 0; w = 1:201; for I = w, for J = w, if norm ([x (I), x (J)] -. 5) <.5, s = s + 1; end, end, end, s, t2 = toc

дає те ж саме s і t2 = 7.47 сек, так що t2/t1 = 46. Це зайвий раз говорить про те, що потрібно розумно підходити до використання операторів мови програмування.

4. Графічний спосіб розв'язання рівнянь

1. Простий приклад: знайти корені рівняння x * sin (x ^ 2) = 0 на відрізку [0,3]. Програма:

1; x = 0: .01:3; f = x .* sin (x. ^ 2); plot (x, [f; 0 * f]), grid

2; ginput

У команді ginput точка знімається натисканням лівої клавіші миші, Enter - вихід з ginput.

Перевіримо це іншим способом:

3; nx = length (x); w = 1: nx -1; x (find (f (w) .* f (w +1) <0 | f (w) == 0)) Відп: 0, 1.77 , 2.5.

Цей рядок можна спростити:

4; nx = length (x); w = 1: nx -1; x (f (w) .* f (w +1) <0 | f (w) == 0)

Матриці та вектори з елементами 0-1.

2. Складний приклад - неявні функції. Побудуємо графік неявної функції f (x, y) º x 3 y-2xy 2 + y-0.2 = 0, x, y = [0, 1]. Це виконає програма

1; h =. 02; x = 0: h: 1; [X, Y] = meshgrid (x); f = X. ^ 3 .* Y -2 * X .* Y. ^ 2 + Y -. 2 ;

2; v = [0,0]; contour (x, x, f, v), grid

На графіку зелена лінія (праворуч вона двозначна) представляє шуканий результат. Область у першому квадраті між цими кривими позначимо через G. Це завдання зовсім непросто зробити в інших системах програмування перш за все тому, що обчислення утворюють лінії рівнів точок - у загальному випадку дуже складна процедура.

З'ясуємо, який знак має f в області G, для чого виконаємо

3; mesh (x, x, f .* (f> 0))

Це приклад тривимірної, тобто xyz-графіки. У ній колір використовується для зображення амплітуди (значення z),

змінюючись з зростанням z від темно-синьому через блакитний, зелений і жовтий до темнокрасно.

Обчислимо площу S цій галузі:

4; S = h ^ 2 * sum (f (:)>= 0) (S = 0.7296).

Для h = 0.01 виконаємо рядок 1, потім рядок 4 і отримаємо S = 0.7204, а для h = 0.005 знайдемо S = 0.7152. При інтегруванні завжди природно робити такі перевірки.

З'ясуємо, який обсяг укладений між поверхнею f (x, y) і областю G, де f (x, y)> = 0. Для цього знову візьмемо в рядку 1 h = 0.02 і обчислимо

5; V = h ^ 2 * sum (f (f> = 0)) (V = 0.1268)

Для h = 0.01 V = 0.1235, а для h = 0.005 V = 0.1219. Тепер не потрібно писати f (:), оскільки f (f> = 0) є вектор.

Звичайно, ці результати наближені (з точністю до 1 - 2%), але відзначте, як швидко і просто вони були отримані. Такі прийоми можна застосовувати для вирішення досить широкого кола завдань.

Виконаємо рядок

6; C = contour (x, x, f); clabel (C)

яка зашле числову інформацію про графік у матрицю C і побудує графік, вибравши значення рівнів автоматично. З матриці C можна послідовно вибирати все криві.

Узагальнення. Графічним способом можна розв'язувати системи рівнянь та рівняння в комплексній площині. Команда contour 3 будує лінії рівнів для функцій f (x, y, z), при цьому сітки з аргументів завжди повинні бути прямокутними.

5. Поліноми

За ступенем придатності, за різноманітністю і якості відповідних команд скалярні поліноми - наступні за матрицями математичні об'єкти в MATLAB'е. Поліном

p (x) = a n x n + a n-1 x n-1 +...+ a 0 задається вектором-рядком p з чисел a n, a n-1, ... , A 0,

тобто коефіцієнтами, розташованими в порядку убування показника ступеня. Його ступінь n ставити не треба, оскільки n = length (p) -1; поліном може бути і константою - тоді n = 0; коефіцієнти a k - будь-які комплексні числа. Вектор p інтерпретується системою як поліном тільки тоді, коли він задається як параметр для однієї з команд, які виробляють обчислення з поліномами. Так як в цих командах не перевіряється умова a n ¹ 0, треба намагатися самим дотримуватися його, оскільки іноді це може служити джерелом помилок.

Основні команди для дій з поліномами такі:

conv (p, q) - твір поліномів p і q. Назва команди походить від слова convolution (згортка), оскільки коефіцієнти твори дійсно виходять як компоненти згортки векторів p і q.

[Q, r] = deconv (b, a) - приватне (q) і залишок (r) від ділення b на a, так що conv (a, q) + r = b.

residue (b, a) - розкладання раціональної функції b (x) / a (x) на елементарні дроби над полем комплексних чисел з виділенням цілої частини. Якщо a (x) має кратні або близькі один до одного коріння, результати можуть бути невірними, оскільки таке завдання погано обумовлена. Погана обумовленість, тобто вкрай сильна залежність результату від коефіцієнтів, ілюструється заключним прикладом з цієї теми.

p = poly (r) - побудова полінома по корінню, заданим у векторі-стовпці r. Для квадратної матриці r поліном p буде її характеристичним многочленом.

polyval (p, x) - поелементне обчислення значень полінома p на безлічі x, де x може бути як вектором, так і матрицею. Розміри результату збігаються з size (x).

polyder (p) - похідна від p.

roots (p) - вектор-стовпець, що містить всі корені полінома. Їх порядок не визначений. Сказане з приводу нестійкості результатів для команди residue точно так само відноситься і до команди roots. Коріння полінома обчислюються як власні значення деякої матриці A (p) того ж порядку.

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

1. Побудуємо графік полінома p (x) = x 3-x +2 на відрізку -1 <= x <= 1. Це виглядає так:

p = [1,0, -1,2]; x =- 1: .01:1; f = polyval (p, x); plot (x, f), grid

Поліном не має коренів на заданому відрізку. Це підтверджує і команда

roots (p) '

яка дасть

ans = -1.5214 0.7607 - 0.8579 i 0.7607 + 0.8579 i.

2.Разделім попередній поліном на x-3:

[Q, r] = deconv (p, [1, -3])

Тоді

q = 1 3 8, r = 0 0 0 26.

Іншими словами, приватне q (x) = x 2 +3 x +8, а залишок r = 26.

3. Розкладемо функцію (x-3) / p (x) на елементарні дроби:

1; [r, s, k] = residue ([1, -3], p); r ', s', k '

Для r 'отримаємо вектор з трьох компонент r 1, r 2, r 3:

-0.7607 0.3803 - 0.3803 0.4289i + 0.4289i,

для s '- також вектор з трьох компонент s 1, s 2, s 3:

-1.5214 0.7607 - 0.7607 0.8579i + 0.8579i

і k = [] (це означає, що цілої частини в розкладанні немає - дійсно, у чисельника перша, а у знааменателя третього ступеня). Компоненти векторів r і s означають, що

(X -3) / p (x) = sum (r i / (x - s i), i = 1:3.

Команда residue працює і у зворотний бік:

2; [q, p] = residue (r, s, k)

відновить вихідні чисельник і знаменник:

q = 0 1 -3 (він вийшов точно),

p = 1.0000 -0.0000 -1.0000 2.0000 (тут вже позначилися помилки округлення і старший коефіцієнт не дорівнює 1 автоматично).

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

Перейдемо тепер до опису прикладу. Нехай vn = 1: n, де n> 1 - цілочисельний параметр, і pn = poly (v ') - поліном з корінням 1: n, які добре відокремлені один від одного, а wn = roots (pn) - вектор-стовпець з обчисленими корінням полінома pn. Проведемо порівняння vn 'і wn для різних n. Почнемо з n = 2:

1; n = 2; vn = 1: n; pn = poly (vn '); wn = roots (pn); [vn', wn]

і отримаємо ans = 1 2 2 1

звідки видно, що елементи в wn потрібно впорядкувати. Виконуючи при n = 3 відредаговану рядок 1

2; n = 3; vn = 1: n; pn = poly (vn '); wn = roots (pn); R = [vn', sort (wn)]

знайдемо R = 1.0000 1.0000

2.0000 2.0000

3.0000 3.0000.

З'явилися в першому стовпці R цифри 0 як би "наведені" значеннями з другого стовпця, і таких дрібних шорсткостей у команди format чимало. А цифри 0 у другому стовпці R говорять про те, що вже з'явилася похибка у визначенні коренів. Вона, звичайно, ще дуже мала:

3; ((R (:, 2) - R (:, 1)). / R (:, 1)) '

дає відносну помилку для коріння

ans = 1e-14 * (0.1110 -0.0444 -0.1184)

- Це означає, що в чисельному результаті вірні приблизно 14 знаків.

Для n = 10 виконання рядки (її можна створити з рядків 2 і 3)

4; n = 10; vn = 1: n; pn = poly (vn '); wn = roots (pn); R = [vn', sort (wn)]; R1 = (R (:, 2)-R (:, 1)). / R (:, 1)

говорить про те, що точність результату поступово падає. Рядок

5; me = max (abs (R1))

дає me = 4.2009e-10, тобто тепер для всіх коренів вірні тільки 9 знаків (m e - max. error). Але коріння ще залишаються речовими, оскільки

6; iwn = sum (abs (imag (wn)))

дає iwn = 0.

Виконаємо рядок 4 для n = 20 і потім рядком 5 знайдемо максимальну відносну помилку me = 0.0073, так що тепер для всіх коренів вірні тільки 2 знаки, але результат поки що залишається речовим, оскільки після рядка 6 отримаємо iwn = 0. Порівняємо точні і обчислені коріння графічно: графік

7; plot (R), grid

показує, що похибка для деяких коренів вже видно на око - для коренів 10:17 жовта і фіолетова лінії трохи розрізняються.

Тепер приступимо до самого цікавого в даному прикладі. При n = 20 pn (2) = -210 (це коефіцієнт при x 19). Додамо до нього 1e-7, тобто внесемо до нього мале обурення приблизно в 10-му знаку, і повторимо розрахунок з відредагованою рядком 4:

8; n = 20; vn = 1: n; pn = poly (vn '); pn (2) = pn (2) +1 e-7; wn = roots (pn); R = [vn', wn], plot (R), grid

Незважаючи на таке мале обурення в коефіцієнті pn (2), деякі корені стали комплексними. Це видно, по-перше, з видачі на екрані (їх уявні частині сягають за модулем 2.7), по-друге, з того, що тепер рядком 6 отримаємо iwn = 18.67, і з графіка. Якщо побудувати графік

9; plot (R ,'.'), grid

тобто прибрати з'єднання між точками, результат буде виглядати більш рельєфно. На таких графіках режим zoom працює.

З'ясуємо, чому так сильно змінилися результати при внесенні такої малої обурення. Позначимо через p (x) наш незбурений поліном pn при n = 20 і через a його другий коефіцієнт:

p (x) = prod (xk), k = 1:20, або p (x) = x 20 + ax 19 + ... +20! , A =- 210.

Тоді для коренів x = 1:20 за теоремі про похідної неявної функції будемо мати

p / x * dx / da + p / a = 0,

звідки

dx / da = - p / a / p / x.

У нас p / a = x 19, а поліном p / x знаходиться як polyder (pn) (див. початок цієї теми). Тому для обчислення dxda на безлічі vn наших коренів спочатку виконаємо відредаговану рядок 4 з n = 20

10; n = 20; vn = 1: n; pn = poly (vn ');

а потім рядок (на графіку значення функції визначені тільки для x = 1:20)

11; dpn = polyder (pn); dxda =- (vn. ^ 19). / Polyval (dpn, vn); plot (log10 (abs (dxda )),'.'), grid

і побачимо, що вже при x = 8 dx / da = 10 5 і буде ще більше з зростанням x. Тому внесення в коефіцієнт a обурення 10 -7 має в обов'язковому порядку помітно позначитися на значеннях деяких коренів, яким би методом вони не знаходилися. Більш того, якщо ці необхідні зміни ніяк не проявилися, метод слід забракувати.

6. Ітерації

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

x k +1 = F (x k), k = 0, 1, 2, ... , X 0 задано,

де F - задана перетворення y = F (x), x 0 - якось вибраного початкове наближення, x k - значення змінної x на k-й ітерації, а сама змінна x може бути будь-який - числом, вектором або матрицею. Межа ітерацій, якщо він існує, будемо позначати через X, і тоді рівняння

X = F (X) (1)

має обов'язково виконуватися для ітераційного перетворення F. Ця обставина допомагає вибирати різні варіанти для F. Всі рішення X цього рівняння називаються нерухомими точками перетворення F (x). Всі такі x 0, для яких послідовність x k сходиться, утворюють область збіжності ітераційного перетворення F. Швидкість збіжності ітерацій x k характеризується поведінкою числової послідовності

v k = norm (x k +1-x k) / norm (x k-x k-1),

де норма може вибиратися по-різному. Для збіжності x k вже не обов'язково, щоб існував межа V послідовності v k, але дуже часто для збіжних ітерацій він існує, і якщо це так, то нехай a = abs (V). Тоді при a = 1 збіжність x k буде вкрай повільною, при 0 <a <1 це буде збіжність за геометричній прогресії зі знаменником q = V, а при a = 0 послідовність x k буде збігатися швидше будь-якої геометричної прогресії. Покажемо тепер на найпростіших прикладах, як все це виглядає на ділі. Щоб мати можливість різноманітності варіантів, завдання візьмемо нелінійної. Рвссмотрім рівняння

(X-2) (x-3) = 0 або f (x) º x 2-5x +6 = 0 c корінням x 1 = 2, x 2 = 3, (2)

побудуємо для знаходження його рішень x 1 = 2, x 2 = 3 кілька ітераційних перетворень чи схем F і проаналізуємо їх роботу.

1.Пусть для рівняння (2)

x = F (x), де y = F (x) = (x 2 +6) / 5.

Обов'язкова умова (1) для перетворення F виконано, однак при цьому переході з'явилося додаткове рішення x = inf (μ). Втрата деяких або придбання нових рішень часто трапляється при переході від вихідного рівняння до його ітераційної формі. Переходячи до потрібної нам індексного запису, будемо мати

x (k +1) = F (x (k)), k = 1: n,

де початкове наближення x (1) і число ітерацій n повинні бути задані. Будемо накопичувати значення x (k) у змінній x, а поточне значення x (k) позначимо через xt. Потрібні обчислення реалізує рядок

1; xt = 0; n = 100; x = xt; TF = '(xt ^ 2 +6) / 5 "; for k = 1: n, xt = eval (TF); x = [x, xt]; end, plot (x), grid

Тут записані текстова змінна TF і команда eval (вона інтерпретує TF як вираз xt ^ 2 +6) / 5 і виконує його). Після виконання рядка 1 на графіку відобразиться 101 значення x (k), включаючи початкове наближення. Ітерації зійшлися до px = 2. Рядок

2; xt = 0; n = 100; x = xt; TF = 'xt = (xt ^ 2 +6) / 5;'; for k = 1: n, eval (TF), x = [x, xt]; end, plot (x), grid

справить ті ж обчислення, але зверніть увагу на те, як у ній записані TF і eval.

Хоча зовні збіжність x (k) до x 1 = 2 не викликає сумнівів, це насправді завжди потрібно перевіряти ретельніше. Так як y '= F' (x) = 2x / 5, то F '(2) = 0.8 <1, а F' (3) = 1.2> 1, але ітерації, як відомо, можуть зійтися тільки до такої нерухому точку X перетворення F, для якої | F '(X) | £ 1. Звідси ясно, чому X = 2. Проте далеко не завжди можна знайти F '(x) аналітично, і тому в загальному випадку доводиться використовувати обчислювальний підхід. При відомих вже x (k) його реалізує рядок

3; w = 2: n; v = (x (w +1) - x (w ))./( x (w) - x (w -1)); plot (v), grid

З графіка видно, що v (k) дійсно сходяться до a = 0.8, як і має бути згідно теорії, яка тут виглядає очевидною.

Тепер будемо варіювати початкове наближення, виконуючи рядка 1 і 3.

(A) xt = 1.5, 1.9, 1.99, 1.9999. При останніх двох значеннях xt на графіку з рядка 3 з'являться осциляції справа, так як тут різниці x (k +1)-x (k) і тим самим значення v (k) втрачають точність: x (k) вже зійшлися з багатьма знаками.

(A ') xt = 2. Графік з рядка 2 взагалі порожній, бо тоді всі v (k) = 0 / 0 = NaN.

(B) xt = 2.01, 2.5, 2.9, 2.99, 2.9999. В останньому випадку x (k) спочатку досить довго (до k = 20) затримуються в районі нерухомої точки x 2 = 3 (вони весь час йдуть від неї, але на графіку цього не видно), а потім приблизно за 60 ітерацій монотонно рухаються до x 1 = 2.

(B ') xt = 3. Всі x (k) = 3, а графік з рядка 3 порожній. Це сталося тому, що при xt = 3 F (xt) = 15 / 5 = 3 виходить без помилок округлення.

(B'') xt = 3.01. Межами для x (k) і v (k) буде inf, так що x 2 = 3 є нестійкою нерухомою точкою перетворення F: при найменшому зсуві x 0 з x 2 в межі ітерацій вийде яка стійка нерухома точка x 1, або inf - набута нерухома точка. Проварьіруем цей зрушення:

(С) xt = 3 +1 e-15, 3 +1 e-14, 3 +1 e-8. У 1-му варіанті догляду xt не відбувається, у 2-му тенденція відходу вже зародилася і він неминуче станеться зі збільшенням числа ітерацій, в 3-м догляд виявився повною мірою вже на 100 ітераціях.

(З ') xt =- 2.5, -2.9, -2.99, -3, -3.01, -100, 100. При xt =- 3 знову отримаємо x 2 за одну ітерацію, а при xt = -3.01 і далі отримаємо inf. Таким чином, при матеріальному x 0 ітерації сходяться до стійкої нерухому точку x 1 = 2 при -3 <x 0 <3 і до inf при | x 0 |> 3. Прорахуємо той випадок, коли | x 0 | = 3. Цей розрахунок виконується рядком (редагуємо рядок 1)

4; n = 100; fi =- pi: pi/20: pi; xt = 3 * exp (i * fi); TF = '(xt. ^ 2 +6) / 5 "; for k = 1: n, xt = eval (TF); end, plot (xt ,'.')

На графіку (він комплекснозначний) видно 4 точки: точка z = 3, точки z = 3 ± s * i з малим s> 0 і точка z = 2. Щоб розібратися в результаті, виконаємо рядок 4 з n = 1000 і, зробивши вікно MATLAB'а повним, видамо

5; imag (xt ')

Образи першої та останньої точки початковій окружності злегка відрізняються від z = 2, а її 21-а точка (вона відповідає fi = 0) є z = 3. Це сталося тому, що sin (pi) та sin (-pi) не дорівнюють нулю в точності. Знову зробимо вікно MATLAB'а невеликим і виконаємо рядок 4 з радіусом 3.01, а потім рядка

6; sum (isnan (xt))

7; find (isnan (xt))

і отримаємо, що точки первісної xt з номерами 1, 21, 41 звертаються в inf (тут nan = inf / inf). Для радіуса 5 їх буде 21, для радіуса 5.6 їх буде 41, тобто все.

Кордон області збіжності зазвичай важко досліджувати аналітично, навіть у такому простому прикладі, як цей, і сама вона не визначається з умови | F '(x) | = 1 або | x | = 2.5. Умова | F '(X) | <1 гарантує стійкість нерухомої точки X перетворення F (x), умова | F' (X) |> 1 є достатнім для її нестійкості, а в разі | F '(X) | = 1 вона може бути як стійкої, так і нестійкою. Нестійкі нерухомі точки порівняно рідкісні в обчислювальних завданнях, але їх корисно мати на увазі при дослідженні чисельних алгоритмів. У нашому випадку ми встановили тільки, що безліч | x | £ 3 належить області збіжності ітерацій F, але зовсім не описали межу цієї області. Ми встановили також, що при | x | ³ 5.6 ітерації сходяться до inf, і тому inf є стійкою нерухомою точкою F. Щоб отримати уявлення про межі області збіжності, виконаємо рядки

8; r = 3: .1:5.6; z = []; n = 100; for kr = r, xt = kr * exp (i * fi); for k = 1: n, xt = eval (TF); end, z = [z; xt]; end

9; zn = isnan (z); Z = r '* exp (i * fi); plot (Z (zn)), axis equal

і побачимо наближений графік кордону - це зовсім не коло. Комагнда axis e qual вибирає по осях однаковий масштаб (це діє тільки на поточний plot; у команди axis багато інших опцій). Щоб побудувати кордон акуратніше, виконаємо рядок 4 з кроком по fi, рівним pi / 180, потім рядок 8 з r = 3: .02:5.6 і рядок 9. Отримаємо графік кордону, виконавши рядки

10; zn = isnan ([z; z (end ,:)]); zn = diff (zn) ~ = 0; plot (Z (zn)), hold on

11; y = 3 * exp (i * fi); plot (y, 'm'), axis equal, hold off

і побачимо відміну області збіжності від кола | z | £ 3.

Знайдемо ті точки, які перейдуть у z = 3 на перших 10 ітераціях. Для цього доведеться розглянути зворотне до F перетворення x = ± (5 y -6) ^ .5 і зробити з ним 10 ітерацій, задавши початкове значення y = 3. Рядок

12; n = 10; yt = 3; for k = 1: n, yt = (5 * yt-6). ^ .5; Yt = [yt, yt-]; end, plot (yt ,'.')

показує, що при цьому вийде практично та ж лінія, що і в рядку 10. Дивно, що це вірно і для інших точок z з області збіжності перетворення F, але з деякими вкрапленнями всередину області збіжності. Брати n великим не можна, тому що в кінці довжина вектора yt дорівнює 2 n.

2.Теперь представимо (2) у вигляді

x = F (x), де y = F (x) = 5-6 / x, так що F '(x) = 6 / x 2, F' (3) = 2 / 3 <1, | F '( x) | <1 при | x |> 2.45.

Отже, стійка нерухома точка - це x 0 = 3. Отримаємо резултате відразу для "всіх" речових x 0 = xt:

1; n = 100; xt =- 5: .5:5; x = xt; TF = '5 -6./xt '; for k = 1: n, xt = eval (TF); end, plot (x , xt ,'.')

Всі межі дорівнюють 3, випадає тільки нестійка точка xt = 2, для якої ітерації знову йдуть без помилок округлення. Виникає припущення, що вся комплексна площину з виколоти точкою x 1 = 2 стягується итерациями в точку x 2 = 3, хоча не скрізь | F '(x) | <1. Подивимося, як це можна побачити на графіку, виконавши програму

2; n = 4; fi =- pi: pi / 20: pi; xt = 4 * exp (i * fi); TF = '5 -6. / Xt '; z = xt;

3; for k = 1: n, xt = eval (TF); z = [z; xt]; end, plot (z','.'), axis equal

Початкове наближення (коло радіуса 4 з центром в нулі) ітеріруется 4 рази і всі результати видаються на графік. Кола (на графіку їх 5) досить швидко стягуються в точку x 2 = 3. Застосуємо zoom, щоб побачити малі кола.

Зробимо розрізи на початковій окружності:

4; n = 4; fi =- pi: pi/20: pi *. 75; fi (end-4 )=[]; xt = 4 * exp (i * fi); TF = '5 -6./xt '; z = xt;

і знову виконаємо рядок 3. Ми побачимо, що 1-а ітерація змінює напрямок обходу кола на протилежне, інші - ні. Виконаємо рядок 4 з n = 20 і потім рядок 3. За допомогою zoom дійдемо до найменшої кола (вона білого кольору) і переконаємося, що вона лежить правіше точки z = 3 приблизно в смузі від 3.0001 до 3.0004.

Втрат і придбань інших нерухомих точок тут немає.

3. Вирішимо нашу задачу методом Ньютона. Якщо x''задовольняє рівнянню f (x'') = 0, а x 'знаходиться поблизу x''і використовується для наближеної апроксимації f (x''), то

0 = f (x'') = f (x ') + f' (x ') (x''-x') або x''= x'-f (x ') / f' (x ') при умови, що f '(x'') ¹ 0 і тим самим f' (x ') ¹ 0.

Це і є ітерації по Ньютону для рівняння f (x) = 0. При переході до індексного формі запису для f (x) = x 2-5x +6 і f '(x) = 2x-5 з (2) отримаємо

x (k +1) = x (k)-f (x (k)) / f '(x (k)) або y = F (x), де F (x) = xf (x) / (f' (x), F '(x) = 2f (x) / (f' (x)) 2.

Так як f '(x 1,2) ¹ 0, можна використовувати ці ітерації по Ньютону (часто їх називають методом дотичній), а оскільки F' (x 1,2) = 0, слід чекати їх швидкої збіжності і того, що тепер обидва рішення x 1,2 будуть стійкими нерухомими точками ітераційного перетворення F. Нерухомою точкою буде і inf, але вона "не наша", і природа її, як ми побачимо нижче, набагато складніше.

Виділимо TF в окремий рядок і проведемо наші стандартні обчислення

1; TF = 'xt-(xt. ^ 2-5 * xt +6). / (2 * xt-5)';

2; xt = 0; n = 100; x = xt; for k = 1: n, xt = eval (TF); x = [x, xt]; end, plot (x), grid

3; w = 2: n; v = (x (w +1)-x (w ))./( x (w)-x (w-1)); plot (v), grid

Межа ітерацій при x 0 = 0 дорівнює 2, а збіжність має порядок вище першого, оскільки v (k) до втрати ними точності встигли подоойті до нуля. Щоб визначити порядок збіжності, відредагуємо рядок 2 на предмет оцінки квадратичної збіжності:

4; w = 2: n; v = (x (w +1)-x (w ))./( x (w)-x (w-1)). ^ 2; plot (v (1:6) ), grid

Тепер до втрати точності v (k) встигають підійти до 1 (біля v (7) точність вже втрачена), так що збіжність ітерацій тут квадратична. Нагадаємо, що точність втрачається у v (k), але не у x (k). Щоб побачити втрату точності у v (k), виконаємо рядок 4, видаючи в plot 7 точок.

При x 0 = 4 межа ітерацій дорівнює 3, а v (k) встигають підійти до -1 (у v (6) точність вже втрачена). Таким чином, в залежності від початкового наближення x 0 = xt ітерації можуть сходитися до обох рішенням завдання.

Так як тепер перетворення y = F (x) вже не є дробово-лінійним, розглянемо трансформацію комплексної прямокутної області в процесі ітерацій. Створимо таку область xt з 41 2 = тисячу шістсот вісімдесят одна комплексних точок і проітеріруем її (при цьому в командному вікні буде багато повідомлень про поділ на нуль):

5; x =- 10: .5:10; [X, Y] = meshgrid (x); xt = X + i * Y; n = 100; for k = 1: n, xt = eval (TF); end , plot (xt ,'.')

На графіку будуть обидва рішення задачі і деяку множину точок на прямій Re (z) = 2.5. Зробимо тільки 10 ітерацій і вставимо видачу графіка на кожній ітерації:

6; x =- 10: .5:10; [X, Y] = meshgrid (x); xt = X + i * Y; n = 10; for k = 1: n, xt = eval (TF); plot (xt ,'.'), pause (0), end

Кінцевий результат той же, що і при 100 ітераціях, але видно динаміка його формування.

Полуплоскость Re (z) <2.5 стягується итерациями в точку x 1 = 2, а полуплоскость Re (z)> 2.5 - в точку x 2 = 3: пряма Re (z) = 2.5 переходить сама в себе. Точка x 1 = 2 має червоний колір тому, що в неї перейшли перші 25 стовпців матриці xt, а залишок rem (25,7) = 4, але 4-й колір - червоний. Далі йдуть зелені точки (26-й стовпець) і синя x 2 = 3, оскільки rem (41,7) = 6 і 6-й колір - синій. Але завжди краще перевірити такі висновки чисельно: знайдемо

7; z = xt (:); [sum (abs (z -2) <.1), sum (abs (z -3) <.1), sum (abs (real (z) -2.5) <.1 )]

(= 1025 = 41 * 25), (= 615 = 41 * 15), (= 38 <41 на 3).

Так як втрат у матриці-таблиці xt MATLAB не допускає, з'ясуємо, у що перейшли ці 3 точки - в inf або NaN:

8; z = xt (:, 26); [sum (isinf (z)), sum (isnan (z))] = 0, = 3.

Індекси цих трьох точок з 26-го стовпця знаходяться командою

9; find (isnan (z)) '= 20 21 22,

і це є точки

z 1 = 2.5-0.5 i, z 2 = 2.5, z 3 = 2.5 +0.5 i.

Тепер розберемося, як перетвориться итерациями F пряма Re (z) = 2.5. Так як пряма Re (z) = 2.5 переходить в себе, її можна параметризовані за допомогою змінної y = Im (z), і нехай y k = Im (F k (Re (z) = 2.5), k = 1, 2, ..., тобто після k ітерацій y переходить в y k (y). Ясно, що на k-й ітерації в inf перейдуть тільки ті точки з усієї безлічі y k-1, що дорівнюють нулю. Оскільки-inf < y <inf, на 1-й ітерації це станеться тільки з однією точкою z 1 = 2.5 (для неї y = 0). Видамо на графік y 1 після 1-ї ітерації:

10; xt = 2.5 + (-10: .1:10) '* i; y = imag (xt); n = 1; for k = 1: n, xt = eval (TF); end, plot (y, imag (xt)), grid

і побачимо, що функція y 1 (y) перетинає вісь абсцис тільки два рази (при цьому y 1 двічі безперервно пробігає від-inf до inf), так що на 2-й ітерації в inf перейдуть тільки дві точки (це вже відомі нам z 1 = 2.5-0.5i, і z 3 = 2.5 +0.5 i). Виконаємо рядок 10 з n = 2 і знімемо з графіка рядком

11; q = ginput; q (:, 1) '

чотири точки перетину кривої y 2 (y) з віссю абсцис: наближено це -1.2079 -0.2056 0.2055 1.2079,

тоді як для y 1 (y) це були значення -0.5 і 0.5. Кожна переходить в inf на k-й ітерації точка породжує ліворуч і праворуч від себе дві такі точки, які перейдуть у inf на (k +1)-й ітерації. Тому із зростанням k точки y k (y) = 0, з одного боку, йдуть в обидва боки від нуля, а з іншого - все частіше з'являються в тих місцях, де вони одного разу вже з'явилися. Всього на k-й ітерації в inf перейде 2 k-1 точок. Насправді з зростанням k вони всюди щільно заповнюють пряму Re (z) = 2.5, а між ними y k (y) безперервно пробігає усі значення від-inf до inf. Щоб бути в цьому упевненими, виконаємо останню програму

12; hy = 1.0033 e -4; xt = 2.5 + (-1: hy: 1) '* i; w = 2: length (xt); n = 100; for k = 1: n, xt = eval (TF); end

13; pzn = sum (sign (xt (w-1)) .* sign (xt (w)) <0)

яка зафіксує 674 зміни знаку (pzn) у функції y 100 (y) для -1 £ y £ 1 з hy = 1.0033e-4 (з-за недостатньої малості hy далеко не всі вони враховані). Для -10 £ y £ -8 з hy = 1.0033e-4 pzn = 675. Для симетричного відносно точки y = 0 відрізка 8 £ y £ 10 з тим же hy отримаємо pzn = 702. Помилки при послідовному обчисленні F k (z) на прямій Re (z) = 2.5 будуть недійсними, тому що значення y k (y) обчислюються з високою точністю, але від кроку hy кількість знайдених змін знака у y k (y), звичайно, залежить: при k = 100 число всіх нулів функції y k (y) дорівнює 2 100-1 = 6.3383e29.

Зробимо коротке резюме щодо нерухомої точки inf перетворення F. Її слід розглядати як нестійку, оскільки будь-яка її околиця з розрізом Re (z) = 2.5 в комплексній площині з зростанням k "йде" від неї, а з цього розрізу все більше точок потрапляє в неї, і ці її дискретні прообрази все щільніше (в межі всюди щільно) заповнюють цей розріз, але всі їх безліч, оскільки воно лічильно, має міру нуль. Звідси випливає, що на пряме Re (z) = 2.5 майже всюди не існує теоретичної межі ітерацій. Через помилки округлення, хоча вони і малі, прообрази inf майже ніколи не потрапляють в inf при обчисленнях і тому ведуть себе так само (тобто хаотично), як і не прообрази.

Зауважимо тепер, що порядок дій в F можна переставити:

F (x) = x / 2 +1.25 +0.25 / (2x-5) і F (x) = (x 2 -6) / (2x-5) - це дві інші математично еквівалентні запису F. Виконаємо рядок

14; TF = 'xt / 2 +1.25 + .25. / (2 * xt -5)';

а потім рядки 5 і 7 і отримаємо той самий результат, що й вище. Але рядок

15; TF = '(xt. ^ 2-6). / (2 * xt-5)';

і рядки 5, 8 дадуть тільки стійкі нерухомі точки перетворення F і, як і раніше, три точки NaN, так що з такою TF ми не побачили б розглянутої вище картини. Виконаємо рядок 6 з n = 100 і побачимо, як все відбувається до кінця.

Приклад 3 показує, що, користуючись досить простими командами MATLAB'а та керуючись здоровим глуздом, можна проаналізувати досить тонкі математичні деталі завдання.

7. Системи лінійних алгебраїчних рівнянь

У MATLAB'e є кілька команд для вирішення таких завдань. Розглянемо тут найбільш вживану з них - оператор \ (зворотний слеш або ділення ліворуч). Для чисел a \ b = a -1 b на відміну від того, що a / b = ab -1. Для невиродженої квадратної матриці A і вектора-стовпця f вектор-стобец u = A \ f є рішення системи лінійних рівнянь Au = f, тобто в записі по аналогії з числами u = A -1 f, але матриця A -1 в дійсності не обчислюється, а система Au = f вирішується методом виключення Гауса (це значно дешевше), з яким кожен з нас якось знайомий ще зі школи . Права частина f може бути і матрицею, і тоді кожен стовпець матриці u = A \ f є рішення системи для відповідного стовпця з f, тобто командою A \ f можна відразу вирішити багато систем, і чинити так завжди вигідно, бо тоді матриця A буде розкладатися на множники за методом Гауса тільки один раз. Якщо матриця A не квадратна, вектор u = A \ f є рішення системи Au = f в сенсі методу найменших квадратів, але ми не будемо тут розбирати цей варіант оператора \, оскільки не повинні скільки-небудь помітно вдаватися в деталі чисельних методів. Таким оброазом, оператор \ - одна з найпотужніших команд системи. А оскільки рішення систем Au = f з квадратними матрицями - найбільш часто зустрічається в обчислювальній практиці завдання, необхідно мати уявлення хоча б про один з методів її реалізації.

1.Решім і проаналізуємо на точність систему Au = f, виконавши наступний рядок:

1; x = 1: .1:5; m = length (x); A = toeplitz (exp (x)); ut = sin (x) '; f = A * ut; u = A \ f;

Тут задається вектор-рядок x = 1: .1:5 з 41 елемента, потім по вектору exp (x) командою toeplitz обчислюється квадратна матриця A розмірів 41 * 41, далі у вигляді вектора-стовпця задається точне рішення ut, по ньому будується права частина f і, нарешті, знаходиться чисельне рішення системи Au = f. З'ясуємо вид A, побудувавши графіки

2; mesh (A) plot (A (1,:)) plot (A (20,:)) plot (A (41,:))

Тепер виконаємо команди

3; plot (ut) plot (f) plot ([ut, u])

і побачимо, що f сильно відрізняється від u, а точне і чисельне рішення графічно збігаються (зображає u фіолетова лінія накрила жовту лінію ut). Іноді потрібно порівнювати сильно різномасштабні криві (у нас це u і f). Для такого порівняння слід провести нормування кожної з них на свій абсолютний максимум. У нашому випадку графік

4; plot ([u / max (abs (u)), f / max (abs (f))])

дозволяє зобразити динаміку зміни кривих u і f відносно один одного.

2.Команда det обчислює визначник. Знайдемо

1; max (abs (u - ut)) (= 9.7167 e -13) det (A) (= 4.1058 e -9)

звідки видно, що u і ut збігаються приблизно в 12 знаків, хоча визначник det (A) системи на перший погляд дуже малий. Однак спробуємо вирішити нашу систему за правилом Крамера, позначивши таке рішення через uc:

2; d = det (A); uc = zeros (m, 1); for k = 1: m, uc (k) = det ([A (:, 1: k-1), f, A (:, k +1: m)]) / d; end, plot ([ut, uc])

Тут при k = 1 A (:, 1: k-1 )=[], а при k = m A (:, k +1: m), оскільки 1:0 = [] і m +1: m = [ ], так що всі матриці C k = [A (:, 1: k-1), f, A (:, k +1: m)], k = 1: m, сформовані правильно. З графіка видно, що так знайдене uc знову збігається з точним рішенням ut. Тепер

3; max (abs (uc - ut)) (= 9.6934 e -13) тобто похибка практично не змінилася. Тривалий час вважали, що прімалих d систему вирішувати чисельно не можна.

3.В дійсності труднощі чисельного рішення системи пов'язані з можливою близькістю до нуля власних значень матриці A, які визначаються як всі рішення l 1, ..., l m характеристичного рівняння det (A-l E) = 0, де E = eye ( m) - одинична (з нулями поза головної діагоналі) матриця порядку m. За останні 10 років були створені досить хороші програми для знаходження власних значень. У MATLAB'е це робить команда eig, яку ми і застосуємо сайті:

1; sp = eig (A); plot (sp), [y, yi] = min (abs (sp)), [z, zi] = max (abs (sp)), c = z / y

З графіка видно, що всі власні значення речовинні (по осі абсцис відкладені значення індексу) і ніяк не впорядковані. Найближче до нуля (= 0.136) третя з чисел l k, найдалі (= 898) - 41-е, а модуль їх відносини c (A) º c = 6.6040e +3 називається числом обумовленості (condition number) матриці A і відображає набагато більш глибокі її властивості, ніж величина її визначника (яка по суті взагалі нічого не відображає ні в практичному, ні в теоретичному плані). При зверненні [V, D] = eig (A) в діагональної матриці D розташовані l k (раніше вони були у векторі-стовпці sp), а в матриці V у вигляді стовпців - відповідні їм власні вектори v k з одиничною середньоквадратичної нормою

sum (v k 2)) 1 / 2 = 1, k = 1: m.

Команда eig (і цілий ряд заснованих на ній) має непогану точність при вирішенні спектральних задач середньої складності і є дуже інформативною. У нашому прикладі є відносно близькі один до одного l k. Щоб у цьому переконатися, виконаємо рядок

2; ssp = sort (sp); v = 1: m-1; di = diff (ssp); min (abs ([di. / ssp (v); di. / ssp (v +1)]))

(Нормировка різниць робиться на кожен член пари) і отримаємо 0.0044, тобто є така пара, у якої збігається більше двох знаків, так що наш приклад в спектральному плані не зовсім тривіальний. Тепер зробимо матрицю A виродженої шляхом дублювання її рядків і подивимося, як це відіб'ється на результатах eig:

3; r = zeros (1, m-1); for k = 1: m-1, v = [1: k, k, k +2: m]; C = A (v,:); r (k ) = min (abs (eig (C))); end, plot (r)

Результати слід визнати непоганими - max (r) = 1e-13 і є навіть декілька чистих нулів.

Щоб отримати уявлення про власні вектори перетворення A, виконаємо рядок

4; [V, D] = eig (A); D = V '* V; D (1: m +1: m ^ 2) = 0; mcv = max (abs (D (:)))

і отримаємо mcv = 7.5460e-16. Це означає, що власні вектори з високим ступенем точності ортогональні, як і повинно бути для симетричної матриці A.

4.Теоретіческая оцінка похибки при вирішенні лінійних систем. Пояснимо сенс числа обумовленості c (A), який має фундаментальне значення для обчислювальних методів лінійної алгебри. Нехай

A * u = f ¹ 0, A * du = df ¹ 0, re (f) = max (abs (df)) / max (abs (f)).

Тоді u ¹ 0 і du ¹ 0 зважаючи невирожденості A. Будемо розглядати df як помилку в f (зважаючи лінійності завдання це не призводить до втрати спільності - вважати вектор df малим не потрібно). Оцінимо величину

re (u) = max (abs (du)) / max (abs (u)) для всіх допустимих f і df,

тобто виникає при вирішенні нашої системи відносну помилку, яка є вельми загальною характеристикою матриці A. Якщо y, z і c - числа, визначені в рядку

1; sp = eig (A); plot (sp), [y, yi] = min (abs (sp)), [z, zi] = max (abs (sp)), c = z / y

(Це рядок 1 з попереднього прикладу), то

max (abs (du)) £ z -1 max (abs (df)),

оскільки du = A -1 u і z -1 - максимальне за модулем власне значення для A -1, причому рівність для max (abs (du)) теоретично можливо. Аналогічно

max (abs (u)) ³ y -1 max (abs (f)),

оскільки y -1 - мінімальне по модулю власне значення для A -1, і знак рівності тут також можливий. Об'єднуючи оцінки для чисельника і знаменника re (u), будемо мати

re (u) £ z -1 / y -1 (max (abs (df)) / max (abs (f))) або re (u) £ (y / z) re (f), тобто re (u) £ c * re (f),

причому рівність обов'язково досягається хоча б для однієї пари f і df. Іншими словами, число зумовленості c (A) є точна верхня межа зростання відносної помилки при вирішенні нашої системи: якщо відносна помилка правої частини дорівнює re (f), то відносна помилка рішення re (u) не перевершить її більш ніж у c (A) разів.

Зважаючи на важливість цього результату корисно розуміти його і на пальцях: df переходить в du зі зростанням компонент не більше ніж у z -1 раз, а f переходить в u із зменшенням компонент не більше ніж у y -1 раз, так що

re (u) £ (z -1 / y -1) re (f).

Ми встановили сенс числа c (A), виходячи із максимум-норми векторів f, коли norm (f) = max (abs (f)), але більш точне розгляд показує, що це вірно і для середньоквадратичної норми

norm (f) = (sum (abs (f. ^ 2)) ^ 0.5.

Команда cond (A) обчислює значення c (A) для квадратної невиродженої матриці A, але в cond спектр sp шукається не для A, а для B = A'A, бо тоді вся спектральну задачу для B стає свідомо речової і, більше того, всі її власні вектори взаємно ортогональні (остання обставина суттєво зменшує середньоквадратична помилка при цих складних обчисленнях). Для отримання наших y і z потрібно лише витягти квадратний корінь з y (B) і z (B) - це також робиться в команді cond (A). Для нашого прикладу c (A) =

2; cond (A) (= 6.6040 e +3)

і збігається зі значенням c з рядка 1, оскільки A - симетрична матриця.

Незважаючи на зовнішню простоту, поняття числа обумовленості було чітко сформульовано лише в середині 1960-х рр.., Тобто приблизно через 15 років після появи перших ЕОМ, а широко використовуватися стало ще пізніше. Для обчислювальних методів воно виявилося важливішим поняття лінійної залежності, яке тепер деяким чином виражається через нього, але ми вже не будемо тут цього уточнювати.

5.Практіческая оцінка похибки. Число обумовленості може бути не завжди підходить для оцінки похибки. Це так при великій n - тоді розв'язок спектральної задачі в команді cond буде занадто дорогим. Це так і в тому випадку, коли помилка df специфічно нерівномірна і має характерний профіль - тоді бажано мати і профіль для поточечно помилок du (k), чого, звичайно, ніяк не може дати cond (A). У таких випадках доводиться вдаватися до безпосереднього моделювання помилок, яке полягає в завданні деякого безлічі випадкових збурень правій частині і вирішенні всіх таких систем з подальшим поточечно описом меж одержані рішень. Проілюструємо цей спосіб оцінки помилки на нашому прикладі.

Відновимо цей приклад:

1; hx =. 1; x = 1: hx: 5; m = length (x); A = toeplitz (exp (x)); ut = sin (x) '; f = A * ut; u = A \ f;

(Тут введено крок hx). Нехай помилка

df (x) в f (x) є рівномірно розподілена випадкова величина,

за модулем не перевершує значення g (x) = 1e-4 * ò f (t) dt в межах від 1 до x.

У рядку

2; g = 1e-4 * hx * abs (cumsum (f)); rand ('state', 0); n = 30; V = (1: m) '* ones (1, n); plot ([ f / max (abs (f)), g / max (g)]), grid

задається профіль g = g (x) максимально допустимої помилки в точці x, приводиться в початковий стан лічильник випадкових чисел rand, задається число n = 30 збурень правій частині f і матриця V розмірів m * n з продубльованого n раз стовпця (1: m) '. Графік показує, що найбільші помилки там, де abs (f) мало. У рядку

3; F = f (V) + g (V) .* (2 * rand (m, n) -1); U = A \ F; plot (U), pause, plot ([g (V) .* (2 * rand (m, n) -1), g])

створюється матриця F обурених правих частин з m рядків і n стовпців шляхом додавання до розмноження n раз вектору f збурень у потрібних межах, вирішуються всі системи AU = F і виводяться на графік всі рішення U і всі обурення разом з їх кордоном. З останнього графіка видно, що обурення задано правильно. У рядку

4; ua = max (U, [], 2); ui = min (U, [], 2); plot ([ui, ut, ua])

шляхом обробки U вздовж рядків знаходяться поточечно кордону рішень (ua - верхня, ui - нижня) і будується графік точного рішення і цих кордонів.

Результат здається нам поганим тому, що на графіку

5; plot (F)

збурень просто не видно (масштаб f забив їх), а в U вони відбилися занадто сильно. Але формально він не суперечить оцінкою, пов'язаної з числом обумовленості c (A). Дійсно, без урахування профілю df максимальна помилка me обчислюється як

6; sp = eig (A); me = min (abs (sp ))^(- 1) * max (g)

і дорівнює 0.6064, тоді як

7; max ([ua - ut; ut - ui]) (= 0.5180)

і лише трохи підросте зі збільшенням n (при n = 100 це 0.5540), тобто не перевершить значення me, що і свідчить про формальну узгодженості обох методів оцінки похибки. Щоб остаточно позбутися відчуття якоїсь неузгодженості наших методів, застосуємо критерій cond (A) для середньоквадратичної норми. Для етогог нам потрібно пройтися по всіх стовпцях k = 1: n (зараз n = 30) матриць U і F і обчислити

max ((norm (du) / norm (u)) / (norm (df) / norm (f))).

Запишемо це у вигляді (максимально використовуйте кишеню при наборі рядків)

8; max ((sum ((U-ut (V)). ^ 2). ^ .5./sum (U. ^ ​​2). ^ .5). / (Sum ((Ff (V)). ^ 2). ^ .5./sum (F. ^ 2). ^ .5))

і отримаємо 3.6484e +3, що менше c (A) = 6.6040e +3, так що і тут ми не виявили протиріччя.

Якщо підійти до оцінки похибки спрощено, побудувавши графік

9; gm = max (g); plot (A \ [f-gm, f, f + gm])

то всі три лінії на ньому практично співпадуть. Таким способом можна моделювати помилку, якщо перетворення A -1 монотонно, тобто якщо при f 1 <f 2 обов'язково u 1 <u 2 ілі, наоборот, обов'язково u 1> u 2. Але у нас це не так: на графіку (цей рядок виходить з рядка 9)

10; gm = 100 * max (g); plot (A \ [f-gm, f, f + gm])

жовта лінія (вона відповідає рішенню для правої частини f-gm <f) то вище, то нижче фіолетовою. Саме через немонотонності перетворення A -1 виходить такий помітний розкид у U. За допомогою цього прийому можна швидко з'ясувати монотонність перетворення y = Q (x) (не обов'язково лінійного), що далеко не завжди вдається визначити теоретично: якщо всі відгуки Y = Q (X), xa <X <x+a, a> 0 , лежать між Q (xa) і Q (x + a), то перетворення Q монотонно, і потрібно лише взяти значення a таким, щоб результат був видний на графіку (для цього нам довелося збільшити max (g) в 100 разів).

6.Посмотрім, як зміняться результати нашого прикладу при збільшенні m - порядку матриці A. Виконаємо, не змінюючи змісту задачі, відредаговану рядок 1 попереднього прикладу

1; clear all, hx =. 01; x = 1: hx: 5; A = toeplitz (exp (x)); ut = sin (x) '; f = A * ut; u = A \ f; plot ( [ut, u]), c = cond (A)

Тут крок hx зменшено в 10 разів, так що тепер порядок m = 401 - досить високий; c (A) = 6.0804e5 зросла майже в 100 разів, тобто обумовленість A помітно погіршилася (приблизно в 10 2 разів), але

2; max (abs (u-ut)) (= 1.3153e-9)

ще досить малий, хоча і зріс приблизно в 10 3 разів, тобто більше, ніж c (A). Така розбіжність з теорією хіба що попереджає про те, що навіть при збереженні сенсу завдання збільшення її розмірності не дозволяє автоматично застосовувати критерій числа обумовленості до оцінки помилок округлення. До вибору числа m потрібно завжди ставитися з підвищеною увагою.

Щоб отримати уявлення про власні вектори перетворення A, виконаємо рядок

3; [V, D] = eig (A); D = V '* V; m = length (x); D (1: m +1: m ^ 2) = 0; mcv = max (abs (D ( :)))

і отримаємо mcv = 2.2985e-15, тобто ступінь ортогональності залишається напрочуд високою. Жорданова клітини порядку вище першого можуть бути тоді, коли mcv (A)> 0.99.

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

Висновок

MATLAB - головна система програмування, що дозволяє різко скоротити витрати праці при перевірці алгоритмів та проведенні приблизних розрахунків. Можливість проведення великих розрахунків на MATLAB'е визначається в основному тими витратами часу, на які може піти користувач: тут доводиться вибирати між легкістю і наочністю програмування і представлення результатів, з одного боку, і витратами часу на рахунок - з іншого. Система дуже зручна для освоєння й апробації чисельних методів, що ми і хочемо показати тут перш за все. Саме тому вона рекомендується як одна з основних для фізиків і багатьох інших природно-наукових спеціальностей у провідних американських університетах. Детальне освоєння будь-який великий програмної системи - це досить тривалий процес, основу якого складають індивідуальна робота, і наші заняття покликані дати лише початковий імпульс цьому процесу щодо MATLAB'а. Теми 2 - 4 являють порівняно елементарне введення, а в решті розглядаються більш складні приклади, що показують, як можна використовувати програмні і графічні можливості системи для дослідження чисельних алгоритмів.

Література

1. Using MATLAB. Version 5.2. The Mathworks, Inc., 1997. 531 p. MATLAB 5.2 Product Family New Features. Version 5.2. The Mathworks, Inc., 1998. 202 p.

2. Using MATLAB Graphics. Version 5.2. The Mathworks, Inc., 1997. 372 p.

3. MATLAB Functions Reference (Volumes 1 and 2). Version 5. The Mathworks, Inc., 1998. 819 p., 586 p.

4. Дьяконов В.П. Довідник по застосуванню системи PC MatLab. М., Фізматліт, 1993 112 с.

5. Потьомкін В.Г. Система MATLAB. Довідковий посібник. М., "Діалог-МІФІ", 1997. 350 с.

6. Гультяєв А. MATLAB 5.2. Імітаційне моделювання в середовищі Windows. СПб, "Коронс-принт", 1999, 288 с.

7. Дьяконов В.П., Абраменкова К. В. MATLAB 5. Система символьної математики. М., Нолидж, 1999, 633 с.

8. Лазарєв Ю. Ф. MATLAB 5.х. Київ, Вид. група BHV, 2000, 384 с. ("Б-ка студента").

9. Медведєв В. С., Потьомкін В. Г. Control System Toolbox. MATLAB 5 для студентів. М., "Діалог-МІФІ", 1997, 287 с.

10. Потьомкін, В.Г. Введення в MATLAB. М., "Діалог-МІФІ", 2000, 350 с.

11. Потьомкін, В.Г. Система інженерних розрахунків MATLAB 5.х. У 2-х томах. М., "Діалог-МІФІ", 1999, 366 с., 304 с.

12. Рудаков П.І., Сафонов В.І. Обробка сигналів та зображень. MATLAB 5 x. М., Діалог-МІФІ ", 2000, 413 с. (" Пакети прикладних програм ").

Додати в блог або на сайт

Цей текст може містити помилки.

Програмування, комп'ютери, інформатика і кібернетика | Реферат
188кб. | скачати


Схожі роботи:
MatLab
Програмування в MATLAB
Програма Matlab та її використання
Система математичних розрахунків MATLAB
Теорія кодування в середовищі MATLAB
Заняття з MATLAB в комп`ютерному класі
Рішення задачі за допомогою програм Mathcad та Matlab
Основи роботи в системі символьної математики MATLAB 5 2
Моделювання структурних схем в середовищі SIMULINK пакета MATLAB
© Усі права захищені
написати до нас