Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
Міністерство освіти і науки України
Вінницький національний технічний університет
МАТЕМАТИЧНІ МЕТОДИ
РОЗВ’ЯЗАННЯ ІНЖЕНЕРНИХ ЗАДАЧ
У БУДІВНИЦТВІ
В. І. Риндюк_____________
І.В. Коц _________________
Вінниця ВНТУ 2006
Міністерство освіти і науки України
Вінницький національний технічний університет
МАТЕМАТИЧНІ МЕТОДИ
РОЗВ’ЯЗАННЯ ІНЖЕНЕРНИХ ЗАДАЧ
У БУДІВНИЦТВІ
Вінниця ВНТУ 2006
Міністерство освіти і науки України
Вінницький національний технічний університет
МАТЕМАТИЧНІ МЕТОДИ
РОЗВ’ЯЗАННЯ ІНЖЕНЕРНИХ ЗАДАЧ
У БУДІВНИЦТВІ
Затверджено Вченою радою Вінницького національного технічного університету як практикум для магістрантів будівельних спеціальностей 8.092101- «Промислове та цивільне будівництво» і 8.092108 – «Теплогазопостачання і вентиляція». Протокол № 4 від 25 листопада 2004 р.
Вінниця ВНТУ 2006
Зміст
Вступ .....................................................................................................................5
І Елементи теорії похибок ..............................................................................5
точного врахування похибок .................................................................12
ІІ Математичне моделювання, модель і експеримент ..................................15
2.1 Інтерполяція в математичному моделюванні ......................................16
2.2 Основні характерні риси моделювання ................................................16
2.2.1 Основні типи математичних моделей .............................................18
2.3 Планування експерименту .....................................................................22
2.3.1 Підбір емпіричних формул та згладжування .................................23
2.3.2 Лінійне згладжування .......................................................................25
2.3.3 Нелінійне згладжування ...................................................................26
2.3.4 Лінійна кореляція і регресія .............................................................26
2.3.5 Рівняння прямих регресії ….............................................................29
2.4 Інтерполяція функцій однієї і декількох змінних .................................32
2.4.1 Одновимірна інтерполяція …............................................................33
2.4.2 Інтерполяція сплайнами ....................................................................34
2.4.3 Багатовимірна інтерполяція ..............................................................36
ІІІ Елементи чисельних методів ......................................................................38
3.1 Розв’язання системи лінійних алгебраїчних рівнянь методом
простої ітерації .........................................................................................38
3.1.1 Загальні положення .........................................................................38
3.1.2 Приклад ............................................................................................39
3.1.3 Текст програми розрахунку ............................................................40
3.2 Обчислення дійсного кореня рівняння з даною точністю із
застосуванням комбінованого методу .................................................41
3.2.1 Загальні положення .........................................................................41
3.2.2 Приклад ............................................................................................42
3.2.3 Текст програми розрахунку ….......................................................43
3.3 Наближене розв’язання систем нелінійних рівнянь методом
Ньютона .................................................................................................43
3.3.1 Загальні положення ........................................................................43
3.3.2 Приклад ...........................................................................................45
3.3.3 Текст програми розрахунку …......................................................46
3.4 Точкове апроксимування функцій за методом найменших
квадратів ................................................................................................47
3.4.1 Загальні положення ........................................................................47
3.4.2 Приклад ...........................................................................................48
3.4.3 Текст програми розрахунку ............................................................50
3.5 Одновимірна мінімізація функцій ........................................................51
3.5.1 Загальні положення ..........................................................................51
3.5.2 Приклад .............................................................................................52
3.5.3 Текст програми розрахунку ............................................................53
3.6 Багатовимірна безумовна мінімізація функцій. Градієнтні методи..54
3.6.1 Загальні положення .........................................................................54
3.6.1.1 Градієнтний метод з дрібненням кроку .....................................54
3.6.1.2 Метод покоординатного спуску ..................................................55
3.6.1.3 Приклад ..........................................................................................57
3.6.1.4 Текст програми розрахунку .........................................................61
3.7 Чисельне інтегрування та диференціювання функцій однієї
змінної ....................................................................................................62
IV Елементи варіаційного числення та метод скінчених елементів ............66
4.1 Метод скінчених елементів ...................................................................68
4.1.1 Загальні положення ..........................................................................68
4.1.2 Побудова елемента ...........................................................................69
4.1.3 Побудова інтерполяційних многочленів .......................................70
4.1.4 Одновимірний приклад варіаційного методу скінчених
елементів...........................................................................................71
4.1.5 Двовимірний приклад варіаційного методу скінчених
елементів ..........................................................................................78
V Постановка задач лінійного програмування ..............................................84
5.1 Симплекс-метод .....................................................................................88
5.2 Транспортна задача …………………...……........................................98
5.2.1 Метод північно-західного кута .....................................................101
5.2.2 Метод мінімального елемента ......................................................103
5.2.3 Метод апроксимації Фогеля ..........................................................104
5.2.4 Розв’язання транспортної задачі методом потенціалів...............106
VI Градієнтний метод розв’язування задач нелінійного програмування...110
VII Література....................................................................................................116
Вступ
Важливу роль у формуванні надійних і корисних для науково-технічної практики результатів відіграє етап математичного моделювання й оптимізації об’єктів управління, проектування та їх дослідження. Математичні моделі інженерних задач будівельного матеріалознавства і технології можуть бути розроблені тільки за участю інженерів-технологів з виробництва будівельних матеріалів та конструкцій.
Для цього спеціалістам необхідні знання, які дозволяють сформувати проектні, матеріалознавчі і технологічні задачі в математичних термінах, вибрати ті чи інші шляхи розв’язання математичної задачі і отримати з результатів математичного моделювання оптимальну (корисну) інженерну інформацію. Оптимізація неможлива без формулювання конкретної математичної задачі (наведеної в математичних виразах) з обов’язковими елементами оптимізації (тобто множини варіантів розв’язків і критерію оптимальності).
Даний посібник буде корисним при вивченні дисциплін:
1.1 Розв'язування задачі. Джерела і класифікація похибок
Процес розв'язування будь-якої реальної фізичної задачі можна розбити на такі етапи:
1. Побудова математичної моделі (математичне формулювання задачі), що охоплює найважливіші для даної задачі сторони явища.
Важливо вміти оцінити точність розв'язку задачі, який здебільшого ми дістаємо з похибками. Похибки результатів зумовлюються такими причинами:
1. Математична модель лише наближено відображає реальні явища;
2. Вхідні дані, як правило, — це числа неточні (дані для обчислень часто дістають з експерименту, а кожний експеримент може дати результат лише з обмеженою точністю);
3. Метод розв'язування задачі часто є наближеним. У багатьох задачах точний результат можна дістати лише після нескінченної або досить великої кількості арифметичних операцій, які практично здійснити неможливо. Тому замість точного розв'язку здебільшого доводиться відшукувати наближений (наприклад, замість суми ряду беруть суму скінченої кількості його членів, нескінченний ітераційний процес обривають після скінченого числа ітерацій, інтеграл замінюють скінченою сумою і т. п.);
4. Округлення при обчисленнях. Усі обчислення (вручну і на ЕОМ) можна виконувати лише з обмеженою кількістю значущих цифр. Тому при виконанні арифметичних дій потрібно вдаватися до округлень, які зумовлюють похибки, що нагромаджуються в процесі обчислень.
Похибки, що породжуються вищезгаданими причинами, відповідно називаються:
Повна похибка результату дорівнює сумі всіх перерахованих похибок.
Зауважимо, що похибок, пов'язаних з особливостями будови математичної моделі, не розглядатимемо. Похибки методів досліджуватимемо в кожному окремому випадку.
При виконанні обчислень потрібно дотримуватися такого правила: у проміжних результатах похибка округлення має бути значно меншою від інших похибок. При обчисленнях вручну це правило забезпечується збереженням запасних цифр. На сучасних ЕОМ числа записуються також з достатньою кількістю значущих цифр, тому похибкою окремого округлення здебільшого можна нехтувати на противагу похибці методу і неусувній похибці. У процесі виконання великої кількості операцій похибка округлень збільшується. До значного нагромадження таких похибок може призвести віднімання близьких за величиною чисел, знаходження коренів многочленів високих степенів тощо.
1.2 Похибки наближених чисел
Часто на практиці з тих чи інших причин доводиться замість точного числа брати його наближене значення а. При цьому виникає похибка:
.
Абсолютною похибкою наближеного числа а називають величину
.
Граничною абсолютною похибкою вважають будь-яке число , що задовольняє умову:
.
Це означає, що .
Надалі ці нерівності скорочено записуватимемо так:
.
Зауважимо, що найчастіше точне значення буває невідомим, а тому невідома й похибка наближеного числа а. Наприклад, неможливо точно визначити числа тощо, однак досить легко обчислити межі, у яких вони знаходяться. Можна знайти наближені значення їх та граничні абсолютні похибки. Зрозуміло, що чим менша гранична абсолютна похибка, тим точніше число а наближує число . Та гранична абсолютна похибка не завжди достатньо характеризує точність обчислень (чи вимірювань). Так, при вимірюванні величини, значення якої дорівнює 5107, граничну абсолютну похибку 0,01 можна вважати досить малою, однак при вимірюванні величини, значення якої дорівнює 510-5, похибка 0,01 завелика. Тому часто використовують відносну похибку наближеного числа а, яка означується рівністю:
.
Очевидно, похибку , як і саме число а, не завжди можна визначити. Тому на практиці здебільшого використовують граничну відносну похибку , що визначається умовою:
.
Отже,
.
Якщо , то вважатимемо, що для і справджуються співвідношення:
, . (1.2.1)
Запишемо додатне число а у вигляді скінченного десяткового дробу:
.
Або
(де всі коефіцієнти > 0; а > 0 і менші за число 10).
K-та цифра наближеного числа а називається правильною, якщо абсолютна похибка цього числа не перевищує половини одиниці k-го розряду, тобто коли
.
В іншому випадку цифру k-го розряду називають сумнівною.
Зауважимо, що інколи цифру k-го розряду називають правильною, якщо абсолютна похибка числа не перевищує одиниці цього розряду. Цифру за останнім означенням називають правильною у широкому розумінні, а за першим – у вузькому. Надалі ми вживатимемо поняття правильних цифр в розумінні першого означення.
Значущими цифрами наближеного числа а називатимемо всі його правильні цифри, починаючи з першої зліва, що не дорівнює нулю (її десятковий розряд позначується m ), до першої сумнівної цифри включно. Усі інші цифри називатимемо незначущими.
Записуючи остаточні результати наближених обчислень, незначущі цифри числа відкидатимемо. При цьому кожне число записуватимемо як добуток деякого степеня числа 10 на число, всі цифри якого значущі. Якщо з наближеними числами виконуватимуться обчислення, то в них, крім значущих, треба зберігати ще одну або дві сумнівні цифри.
На практиці при виконанні обчислень часто виникає потреба округляти числа.
Якщо округляється наближене чи точне число а до n значущих цифр, то за правилом доповнення число а замінюють числом а1 з n значущими цифрами так, щоб похибка округлення не перевищувала половини одиниці розряду, що зберігається, тобто, щоб справджувалась умова:
.
Отже, якщо при округленні числа а до n значущих цифр відкидається менше ніж – , то в числі а1 зберігаються без зміни всі перші n цифр числа а. А коли відкидається не менше як – , то в числі а1 остання цифра береться на одиницю більшою, ніж у числі а.
Зауважимо, що внаслідок округлення наближеного числа a до n значущих цифр похибка округлення 0 може досягти величини 0,5·10m-n+1. Якщо врахувати ще й похибку a числа а, то n-а цифра округленого числа може вже бути неправильною. Та коли похибка a досить мала порівняно з максимально можливим значенням похибки округлення 0 (наприклад, a= 0,5·10m-n ), то похибкою a можна знехтувати і n-у цифру числа а після округлення вважати правильною, навіть якщо похибка округлення досягає 0,5·10m-n+1.
Отже, щоб знайти результат з точністю до 0,5·10m-n+1, треба виконувати обчислення так, щоб дістати результат з точністю принаймні до 0,5·10m-n, тобто з n+1 правильними цифрами. Тоді після округлення можна буде вважати, що результат має n правильних цифр.
Легко встановити зв'язок відносної похибки наближеного числа а з кількістю правильних значущих цифр. Якщо додатне наближене число а має п правильних значущих цифр, то його відносна похибка δ задовольняє співвідношення:
,
де – перша значуща цифра наближеного числа а.
Справді,
,
тому
.
Таким чином, за граничну відносну похибку можна взяти
.
Якщо число а має більше ніж дві правильні значущі цифри, то на практиці можна користуватись оцінкою
.
Навпаки, за відносною похибкою δ можна визначити кількість правильних цифр числа а. Для цього знаходимо граничну абсолютну похибку a із співвідношення (1.2.1), після чого кількість правильних цифр числа а легко встановлюється.
Приклади
1. Якою буде гранична відносна похибка, якщо замість числа взяти 1,41?
Маємо: n = 3, = 1. Отже .
2. Знайти наближено так, щоб відносна похибка не перевищувала 0,001 (або 0,1 %).
Маємо: = 7. Отже, повинно справджуватись співвідношення . Звідки і .
Отже, щоб відносна похибка не перевищувала 0,001, треба взяти не менше як з чотирма правильними цифрами.
3. Число 7432 має граничну відносну похибку =0,01 (1%). Скільки в ньому правильних цифр?
Знайдемо граничну абсолютну похибку:
.
Отже, число 7432 має лише одну правильну цифру (друга цифра сумнівна). Його треба записати так: 0,74·104 , або 74·102, або 7,4·103 і т. д.
Вправи
1. Знайти граничну абсолютну та граничну відносну похибки, якщо замість точного числа взято наближене число а: а) = e, а = 2,7; б) = lg 2, а=0,30; в) =, а = 2,236; г) =, а =9,95.
2. Знайти граничну відносну похибку числа а:
а) =273,4 ±0,1; б) = 82 900 ± 500; в) = 0,00126 ±0,00003.
3. Знайти межі точного числа , якщо відома гранична відносна похибка числа а: а = 1427,394; = 0,00005; а = 97,45; = 0,015; а = 389·105; = 3%.
4. Записати число:
а) π з п'ятьма правильними значущими цифрами;
б) з трьома правильними значущими цифрами;
в) 99,995 з чотирма правильними значущими цифрами.
5. Скільки правильних цифр треба взяти в наближенні числа , щоб відносна похибка не перевищувала δ:
=, δ = 0,0001; = , δ = 0,01%; = lg 127, δ = 1%?
6. Число а має відносну похибку δ. Скільки в ньому правильних цифр, якщо:
а) а = 97,432, δ =0,05%;
б) а = 11,4814, δ =0,01;
в) а = 547,92, δ =0,001?
1.3 Загальна формула для обчислення похибки
Нехай потрібно обчислити значення функції при заданих значеннях незалежних змінних.
Якщо при цьому замість точних значень , які нам не відомі, підставляються їх наближені значення (причому , ), то при обчисленні значення функції виникає похибка
.
Припустивши, що функція диференційовна (достатню кількість разів), ∆y можна подати у вигляді:
(1.3.1)
плюс члени другого і вищих порядків відносно ∆xi.
Якщо похибки ∆xi за абсолютною величиною досить малі, то членами другого і вищих порядків у виразі (1.3.1) можна знехтувати. Тоді для абсолютної похибки наближеного значення матимемо:
, (1.3.2)
де через ∆xi позначені граничні абсолютні похибки аргументів. Таким чином,
.
Для граничної відносної похибки δy маємо:
,
або
. (1.3.3)
З формули (1.3.2) випливає, що гранична абсолютна похибка суми наближених чисел a1+a2+…+an… дорівнює сумі їх граничних абсолютних похибок. Для різниці двох чисел a = a1-a2 маємо: .
При цьому відносна похибка різниці може виявитись досить великою, якщо числа а1 та а2 між собою мало відрізняються і їх різниця близька до нуля. У цьому випадку для забезпечення потрібної точності треба мати в зменшуваному і від'ємнику достатню кількість цифр, щоб гранична абсолютна похибка різниці а була меншою від самої різниці. Тому, по можливості, слід уникати віднімання близьких чисел.
Для похибки добутку додатних чисел а = а1а2... аn з (1.3.3) маємо:
,
тобто гранична відносна похибка добутку кількох наближених додатних чисел, що відрізняються від нуля, дорівнює сумі граничних відносних похибок цих чисел.
Легко переконатись, що гранична відносна похибка частки дорівнює сумі граничних відносних похибок діленого і дільника.
Справді, нехай (а1,а2 >0). Тоді з (1.3.3) маємо:
.
Можна переконатися також, що відносна похибка k-то степеня числа а у k раз більша, ніж відносна похибка числа а. Аналогічно, відносна похибка кореня k-ro порядку з числа a у k раз менша, ніж відносна похибка числа а.
Вправи
1. Знайти суму наближених чисел так, щоб усі цифри результату були значущими (усі цифри доданків значущі):
а) а1 = 28,85; а2 = 491,000285; а3 = 1,18294; а4 = 0,18498932; а5 = 9544,2;
б) а1 = 25,4; а2 = 2900,3241; а3 = 9,7·102; а4 = 0,00097; а5 = 469,99999.
з чотирма правильними значущими цифрами.
1.4 Поняття про ймовірну оцінку похибки.
Обчислення без точного врахування похибок
Наближене число а може відхилятися від точного числа в той чи інший бік на певну величину. Тому наближене число можна розглядати як випадкову величину з математичним сподіванням . Сама похибка також є випадковою величиною з математичним сподіванням, що дорівнює нулю, .
Розглянемо суму наближених чисел a=a1+a2+…+an. Гранична абсолютна похибка суми .
Вважатимемо, що випадкові величини мають один і той самий нормальний розподіл ймовірностей з математичним сподіванням, що дорівнює нулю, і середнім квадратичним відхиленням σ. Тоді, як відомо з курсу теорії ймовірностей, математичне сподівання суми випадкових величин дорівнює сумі їх математичних сподівань, тобто . Дисперсія суми незалежних випадкових величин дорівнює сумі дисперсій цих величин, тобто (нагадаємо, що за означенням середнє квадратичне відхилення випадкової величини X дорівнює , де - дисперсія випадкової величини X, а тому .
Сума незалежних випадкових величин з нормальним розподілом ймовірностей також має нормальний розподіл ймовірностей. Як відомо, випадкова величина з нормальним розподілом ймовірностей з імовірністю 0,997 відхиляється в той чи інший бік від свого математичного сподівання не більш як на три середніх квадратичних відхилення. Отже, якщо практично , то .
Таким чином, абсолютна похибка суми наближених чисел a=a1+a2+…+an пропорційна числу , а не числу n, як це випливало з формули (1.3.2). Це пояснюється тим, що практично похибки різних чисел можуть частково «гасити» одна одну, в той час як формула (1.3.2) передбачає, так би мовити, крайній випадок.
Аналогічно приходимо до висновку, що й відносна похибка добутку наближених чисел a1, a2, …, an зростає пропорційно числу , а не числу n.
Точний підрахунок похибок результатів обчислень наближених чисел досить громіздкий. Тому здебільшого на практиці точно не враховують похибок, а користуються правилами підрахунку цифр за В. М. Брадісом:
Зауважимо, що в основу цих правил покладено принцип академіка О.М. Крилова – основний принцип звичайних обчислень, тобто обчислень без строгого врахування похибок: наближене число слід писати так, щоб у ньому всі цифри, крім останньої, були правильними і лише остання була сумнівною, притому не більш як на одну одиницю.
1.5 Обернена задача теорії похибок
Досить часто доводиться розв'язувати таку задачу: визначити допустимі похибки аргументів так, щоб гранична абсолютна похибка функції не перевищувала заданої величини Δу.
Користуючись формулою (1.3.2), помічаємо, що ця задача не визначена, оскільки, по-різному вибравши значення , можна дістати одне й те саме значення Δу.
Отже, задачу можна розв'язувати кількома способами:
1. Доберемо так, щоб в (1.3.2) всі величини були рівними між собою (принцип рівних впливів).
Тоді, очевидно, і .
2. Поклавши, що всі рівні між собою, дістанемо:
3. Вважаючи всі відносні похибки аргументів рівними між собою
та підставляючи в (1.3.2) вирази для , дістанемо:
Визначивши δ, для матимемо:
При розв'язуванні задачі можуть бути використані і деякі інші варіанти.
Аналогічно можна розв'язувати і другу обернену задачу теорії похибок, коли потрібно підібрати граничні абсолютні чи відносні похибки аргументів так, щоб відносна похибка функції не перевищувала наперед заданої величини.
Вправи
1. Визначити точність, з якою треба виконувати обчислення, щоб знайти результат з точністю до 0,5·10-3:
а) , ;
б) , .
2. Корені рівняння x2-2x+lg3=0 треба знайти з чотирма правильними значущими цифрами. З якою точністю треба взяти вільний член рівняння?
ІІ Математичне моделювання, модель і експеримент
Математичне моделювання – це теоретично-експериментальний метод пізнавально-складової діяльності, це метод визначення і пояснення явищ, процесів і систем (об’єктів-оригіналів) на основі створення нових об’єктів математичних моделей. Під математичною моделлю прийнято розуміти сукупність відношень (рівнянь, нерівностей, логічних умов, операторів тощо), які визначають характеристики перебування об’єкта моделювання, а через них і вихідні реакції в залежності від параметрів об’єкта-оригінала W, вхідних впливів , початкових і граничних умов, а також часу.
Математична модель, як правило, враховує лише ті властивості (атрибути) об’єкта-оригінала W, які відображають, визначають і є цікавими з точки зору цілей і задач конкретного дослідництва.
Відповідно, в залежності від цілей моделювання, при розгляді одного і того ж об’єкта-оригінала W з різних точок зору і в різних аспектах останній може мати різні математичні описи і, як наслідок, бути представленим різними математичними моделями.
Визначення математичної моделі найбільш загально і конструктивно сформулював П.Дж. Косн. Суть його така.
Математична модель – це формальна система, що є кінцевою сукупністю символів і досить суворих правил, визначених цими символами для певного об’єкта із деякими власними відношеннями, символами або константами.
Отже, кінцева сукупність символів (алфавіт) і досить сурові правила оперування цими символами (“граматика” і “синтаксис” математичних виразів) приводять до формування абстрактних математичних об’єктів. Тільки інтерпретація робить абстрактний об’єкт математичною моделлю.
2.1. Інтерпретації в математичному моделюванні
Інтерпретація (лат. “interpretatio” – пояснення, тлумачення) визначається як сукупність значень (змістів), що належать яким-небудь способом елементам деякої системи (теорії), наприклад, формулам і окремим символам. В математичному аспекті інтерпретація – це екстраполяція вихідних положень якої-небудь утримуваної системи, вихідне положення якої визначається незалежно від формальної системи.
Відповідно, можна запевняти, що інтерпретація – це встановлення відповідності між деякою формальною і утримуваною системою. В тих випадках, коли формальна система є застосованою (інтерпретованою) до утримуваної системи, оскільки встановлено, що між елементами формальної системи і елементами утримуваної системи існує взаємно однозначна відповідність, всі вихідні положення формальної системи отримують підтвердження в утримуваній системі. Інтерпретація вважається повною, якщо кожному елементу формальної системи відповідає деякий елемент (інтерпретант) утримуваної системи. Якщо вказана умова порушується, має місце частинна інтерпретація.
При математичному моделюванні в результаті інтерпретації задаються значення елементів математичних виразів (символів, операцій, формул) і цілих конструкцій.
Отже, інтерпретація в математичному моделюванні – інформаційний процес перетворення абстрактного математичного об’єкта (АМО) у відповідну математичну модель (ММ) конкретного об’єкта на основі відображення не пустої інформаційної кількості даних і знань, визначеної АМО і названої областю інтерпретації, в конкретну область – інформаційну кількість даних і знань, визначену предметною областю і об’єктом моделювання і названу областю значень інтерпретації.
Таким чином, інтерпретацію слід розглядати як один із основоположних механізмів (інструментів) технології математичного (наукового) моделювання.
Повний спектр елементів інтерпретації, що відображає перехід від АМО-опису до конкретної ММ, включає чотири види інтерпретації: синтаксичну (структурну), семантичну (змістовну), якісну і кількісну.
2.2 Основні характерні риси моделювання
Метою прикладної математики є математичне сприйняття дійсності [1]. З іншої сторони, інженеру-практику, наприклад, важливіше знати, чи витримає його міст заплановане навантаження. Іншими словами, отримати конкретну відповідь на запитання, а не пориватись до більш вищих цілей.
На практиці вихідним пунктом часто є деяка емпірична ситуація, яка ставить перед дослідником “задачу”, на яку необхідно знайти відповідь. Але вживання таких слів як “задача” і “відповідь” може привести до блуду. Перш за все необхідно встановити, в чому саме суть “задачі”. Це пов’язано з тим, що реальні ситуації рідко бувають чітко окресленими, а складна взаємодія з навколишнім середовищем часто робить точний опис ситуації утрудненим. Часто (але не завжди) паралельно з цією стадією постановки задачі йде процес виявлення основних або суттєвих особливостей явища (рис. 2.1), зокрема для фізичних явищ цей процес схематизації або ідеалізації грає вирішальну роль, оскільки в реальному явищі бере участь множина процесів.
Деякі риси явища є важливими, багато інших – несуттєвими. Розглянемо, наприклад, рух математичного маятника. В даній ситуації суттєвим є реальний характер коливань маятника, а не те, що нитка біла, а тягарець – чорний.
Рисунок 2.1
Після того, як суттєві факти виявлені, наступний крок полягає в переведенні цих факторів на мову математичних понять і величин та постулюванні співвідношень між цими величинами. Як правило, це найважча стадія процесу моделювання.
Після побудови моделі її необхідно перевірити на адекватність. По-перше, сама математична основа моделі повинна мати непротиріччя і підкорятись всім звичайним законам математичної логіки. Потім справедливість моделі залежить від її здатності адекватно описати відповідну ситуацію.
Математична модель являє собою спрощення реальної ситуації. Спрощення наступають тоді, коли несуттєві особливості ситуації відкидаються і складна задача зводиться до ідеалізованої задачі, яка піддається математичному аналізу. Саме при такому підході в класичній прикладній механіці виникли блоки без тертя, ідеальні стержні, невагомі нерозтяжні нитки, нев’язкі рідини та багато інших понять. Ці поняття не існують в реальній дійсності, вони є абстракціями. Але їх часто можна з успіхом вважати задовільним наближенням до реальної ситуації.
Це лише одна сторона спрощення. Друга сторона пов’язана в порівнянні порядку різноманітних величин, які фігурують в моделі. Наприклад, зміна деякої величини х протягом часу можна описати рівнянням виду:
.
Нехай в результаті спостережень було зазначено, що член набагато більший за величиною, ніж . В такому випадку можна зекономити багато часу, спростивши рівняння (відкинувши член ), в результаті чого розв’язок знайдеться швидше, і він буде правильно відображати ситуацію. Зауважимо, що більш точне розв’язання даного рівняння в математичному відношенні в дійсності може привести до неправильних висновків.
Після отримання формул або результатів необхідно здійснити зворотний переклад з математичної мови на мову, якою спочатку формулювалась задача. Тобто, необхідно чітко зрозуміти математичний зміст отриманих розв’язків і те, що вони означають на мові реального світу.
Критерієм достовірності моделі служить інженерна або технічно-економічна корисність нової інформації, отриманої за математичною моделлю при наступній перевірці її в системі об’єкта-прототипу. Щоб модель виявилась корисною, вона повинна мати дві основні властивості.
По-перше, модель повинна бути економічнішою (дешевшою, швидкіснішою, безпечнішою тощо), ніж об’єкт-прототип.
По-друге, модель повинна забезпечувати можливість розповсюдження без суттєвих витрат.
При математичному моделюванні використовуються практично всі розділи математики. В основі вибору того чи іншого математичного методу повинні лежати: ціль моделювання, теоретичні уявлення про об’єкт, накопичення інформації про нього, обмеження на матеріальні і часові ресурси тощо.
2.2.1 Основні типи математичних моделей
З врахуванням існуючої практики розв’язання інженерних задач в області будівельного матеріалознавства і технології виділяють п’ять основних типів математичних моделей. В подальшому будемо їх позначати, відповідно, М1, М2, ..., М5 [2].
Моделі М1. В інженерних розрахунках використовується перехід від нескінченно малих до кінцевих змін в реальній системі. Але такий перехід, як правило, настільки аналітично складний, що інтегрування замінюється наближеними методами. В кінці XIX – на початку XX століття сформувався спеціальний науковий напрямок – теорія подібності, яка поєднала диференціальні рівняння з експериментальною інформацією про систему. В зв’язку з розвитком ЕОМ головну роль почали грати чисельні методи розв’язання задач математичної фізики.
Зауважимо, що при складанні таких моделей – рівнянь математичної фізики – використовуються фундаментальні закони природи. Наприклад, рівняння Нав’є-Стокса для нестисливої рідини, що відображає закон збереження енергії:
,
де – оператор Лапласа;
– вектор швидкості руху поля потоку рідини;
– масова сила, яка віднесена до одиниці об’єму ( – вектор при- скорення вільного падіння в полі земного тяжіння);
– похідна субстанції (для стаціонарного руху в напрямку осі х за- мість можна записати , де – є зміна швидкості при переході від однієї точки до іншої);
Р – тиск рідини.
Моделі М2. Мають в своїй основі перш за все деяку інженерну концептуальну модель, яка виражається на початку словами в термінах даної науки (матеріалознавство, технологія). Ця концептуальна модель описується далі абстрактно-знаково за допомогою диференціальних або алгебраїчних рівнянь, геометричних співвідношень, логічних операцій тощо. На рис. 2.2 розглянуто реологічне тіло Максвелла з періодом релаксації, який дорівнює відношенню між ньютонівською в’язкістю 2 і модулем пружності Гука G1 [3]: .
Рисунок 2.2
В основі математичної моделі лежить уява про те, що поведінка реального матеріалу під навантаженням аналогічна поведінці послідовно з’єднаної пружини і поршня, який рухається у в’язкій рідині. Перевага моделі М2 полягає у “фізичній” інтерпретації технологом, або ця інтерпретація вже закладена в концептуальній основі моделі. Недолік моделі М2 полягає в тому, що описати математичними символами можна концепцію, яка далека від реальності. Для використання моделей М2 необхідно визначати практично різні параметри об’єкта (наприклад, в’язкість 2 і пружність G1 конкретних рідин тощо).
Моделі М3. Описують з відомою точністю зв’язок між входами X і виходом Y системи на основі експериментальних даних без аналізу внутрішньої структури системи, яку можна розглядати як “чорний ящик” [4]. Модель М3 відрізняється універсальністю методології збирання експериментальної інформації і будови моделей. До недоліків варто віднести складність різносторонньої інтерпретації параметрів моделей М3 і обмеженість області їх дії з відповідною точністю рамок тих конкретних об’єктів, для яких вони побудовані. Моделі М3 і надалі будуть займати основне місце в інженерній практиці в зв’язку з розширенням випуску та використанням мікроЕОМ і персональних комп’ютерів.
Моделі М4. Досліджують операцію як сукупність дій, направлених на досягнення деякої мети, характерні, головним чином, для техніко-економічних задач (управління запасами, розподілення ресурсів між підприємствами, задачі заміни обладнання і контролю якості тощо).
На рис. 2.3 показано використання моделі М4: фільтрація рідини через бетон розглядається як рух через сітку випадкової конфігурації з запірними вузлами [5].
Рисунок 2.3
Імітаційна модель М5 включає в єдину програму для ЕОМ моделі різних типів, які об’єднані правилами взаємопереходу, а також генератор випадкових чисел, який генерує збурювальні дії середовища на систему (рис. 2.4). На сьогодні в різноманітних галузях науки створюються імітаційні моделі, до складу яких входять програмні блоки, які забезпечують активний діалог людини з ЕОМ [6].
Математичні моделі є ефективним засобом розв’язання задач у всіх основних областях діяльності інженерів-матеріалознавців; технологів – при дослідженні і проектуванні систем, при управлінні їх функціонуванням. В задачах аналізу системи, який направлений на пізнання внутрішнього механізму явищ, які в ній відбулися, особливе значення мають моделі М1 і М2 , природа яких забезпечує поглиблену інтерпретацію результатів моделювання. В той же час в задачах управління системою, коли головною задачею є досягнення конкретної мети функціонування системи в конкретних умовах, пріоритетними стають моделі М3-М5. При проектуванні системи на етапі пошукового синтезу перевага надається моделям М1 і М2, які дають більш загальну інженерну інформацію про систему деякого класу, а на стадії підготовки синтезованої системи до реалізації підвищується значення моделей М3, М4 і особливо М5.
Рисунок 2.4 Принципова схема управління системи і її зв’язки із зовнішнім середовищем
Робота з моделями всіх класів суттєво спрощується, якщо можливе використання їх графічного відображення.
Так однофакторна модель відображається лінією, яка може бути гладкою або кусково-гладкою. Двофакторна може бути відображена сімейством однофакторних функцій Y = f1(x1) при деяких фіксованих значеннях іншої змінної X2 = C = const (і навпаки). Але найбільш зручне її відображення – площинна діаграма, яка складається з ліній рівного рівня відгуку системи (ізоліній). Ізолінії створюються як лінії перерізу поверхні відгуку з горизонтальними площинами, які проведені через відповідні мітки вертикальної шкали Y. Ізолінії з поверхні відгуку проектуються на площину факторного простору і створюють діаграму.
2.3 Планування експерименту
Будь-який експеримент, яким би він не був складним, закінчується поданням результатів, формулюванням висновків і, можливо, дачею рекомендацій іншим особам. Ця інформація може надаватися у вигляді графіків або кривих, математичних формул або номограм, таблиць, статичних даних або словесних описів. Із графіків отримується залежність результату R від змінної х або залежність R від х і y у випадку параметричних кривих. Для отримання залежності R від х, y і z необхідно побудувати ізометричні координати. Графічне зображення більш складних функцій неможливе, оскільки людина не в змозі наочно уявити більш складні співвідношення. Записуючи результати у вигляді формул, можна виразити залежність R від більшої кількості змінних, але лише в небагатьох експериментах одночасно досліджується більше трьох незалежних змінних.
Статичний показник може давати інформацію про всю сукупність даних і про змінність окремих елементів сукупності. Він може давати інформацію про значимість причинного співвідношення або вказувати ймовірність появи конкретної події в майбутньому на основі попереднього досліду.
Подання результатів експериментів в словесній формі завжди було проблемою при проведенні наукових дослідів. Це найнеефективніший спосіб подання результатів, але його не можна ігнорувати повністю. Безумовно, результати деяких експериментів, які проводяться в сучасних фізичних лабораторіях, просто неможливо подати в словесній формі. В техніці такі експерименти зустрічаються порівняно рідко, і значна частина технічного звіту звичайно містить словесні описи і пояснення.
Більшість інженерних експериментів веде до конкретної дії – прийняття рішення, продовження дослідів або визнання невдач.
Найбільш складною проблемою є правильне формулювання питань, пов’язаних з побудовою плану експерименту.
План експерименту – це загальний термін. Він визначає набір інструкцій до проведення експерименту, в яких вказується послідовність робіт, характер і величина змін змінних і даються вказівки щодо проведення повторних експериментів.
Неможливо уявити собі експеримент, в якому основним продуктом не були б дані, які записані в тій або іншій формі. В простих експериментах кінцеві дані можуть бути досить елементарними і не містити взагалі будь-якої помилки. З іншого боку, у виключно складних експериментах може бути низька точність, а кінцевий результат буде містити великі помилки.
При дослідженні теоретичних питань, пов’язаних з проведенням експериментів, важливу роль грає вивчення помилок і невизначеностей. Ні один експеримент неможливо правильно спланувати без детального вивчення цього важливого фактора.
Експериментатор, в якій би області він не працював, майже завжди притримується такої послідовності: на початку відбувається планування, потім додають обладнання, після цього відбуваються досліди і накінець виконується аналіз та складається звіт. При плануванні і придбанні обладнання аналіз помилок повинен бути на першому плані.
Якщо технічні обмеження приводять до досить великої невизначеності або фінансові можливості не дозволяють придбати прилади високої точності, то проведення експерименту неможливе.
План експерименту є компактним і ефективним, якщо завчасно встановлюються інтервали між значеннями змінних. Якщо нас цікавить функціональне співвідношення між незалежною змінною х і залежною змінною у, то цю функцію можна зобразити у вигляді деякої кривої (або прямої) в системі координат (х, у). Така крива складається з нескінченної кількості окремих точок, з яких вибирається деяка кінцева необхідна кількість точок, що характеризують дану функцію. Якщо розглядається функція двох незалежних змінних, то повній сукупності даних відповідає деяка площина.
Таким чином, вибір кінцевої сукупності експериментальних точок є необхідним і важливим етапом планування.
Наступним етапом після проведення експерименту є обробка результатів.
2.3.1 Підбір емпіричних формул та згладжування
При експериментальному (емпіричному) вивченні функціональної залежності однієї величини y від другої величини x проводять ряд вимірювань величини у при різних значеннях величини х. Задача полягає в аналітичному поданні шуканої функціональної залежності, тобто в підборі формули, що описує результати експерименту. Особливість задачі полягає в тому, що наявність випадкових помилок вимірювань роблять нерозумним підбір такої формули, яка точно описувала б всі експериментальні значення. Іншими словами, графік шуканої функції не повинен проходити через всі точки, а повинен, по можливості, згладжувати “шум” (випадкові помилки вимірювань). Зрозуміло, що згладжування “шуму” буде тим більш точним і надійним, чим більша кількість проведених експериментів, тобто чим більше ми маємо надлишкової інформації.
Наприклад, для проведення прямої досить двох точок – (x1, y1) і (x2, y2), якщо ці точки відомі точно. Але при наявності більш або менш значного „шуму”, як на рис.2.5, для тієї ж мети може знадобитися декілька десятків точок.
Рисунок 2.5
Емпіричну формулу звичайно вибирають з формул конкретного типу, наприклад:
(2.3.1)
Іншими словами, задача зводиться до визначення параметрів a, b, c формули (2.3.1), в той же час як вигляд формули відомий наперед з яких-небудь теоретичних міркувань або з міркувань простоти аналітичного подання емпіричного матеріалу.
Позначимо вибрану функціональну залежність через
(2.3.2)
з явним показом всіх параметрів, які необхідно визначити. Ці параметри не можна визначити точно за емпіричними значеннями функції , тому що останні містять випадкові помилки. Мова йде тільки про отримання “достатньо добрих” оцінок шуканих параметрів. Метод найменших квадратів дозволяє отримати задовільні оцінки всіх параметрів . Див. підрозділ 3.4.
Але не завжди, наприклад, з фізичного аналізу задачі відомий вигляд формули, щоб потім визначити невідомі параметри. Іноді доводиться мати діло з задачами, для яких попередній фізичний аналіз не дає точного вигляду формули або дозволяє вибрати формули з досить широкого кола. Крім того, для простоти подальших розрахунків часто віддають перевагу залежності у вигляді многочлена, але при цьому степінь многочлена не задається наперед.
Розглянемо деякі рекомендації щодо вибору формул серії задач, які часто зустрічаються [10].
Якщо кількість експериментальних точок велика, то підбір емпіричної формули може виявитися досить важким: формули з малим числом параметрів можуть давати більші спотворення, а велике число параметрів не зручне для аналізу. Для аналізу важливо ліквідувати “шум” експерименту, зберігаючи інформацію про істинність функції. Для цієї мети використовується згладжування емпіричних даних, тобто заміна даної таблиці експериментальних точок іншою таблицею близьких до них точок, які лежать на достатньо гладкій кривій.
Згладжування відбувається за допомогою многочленів (бажано оптимального степеня), які наближають за методом найменших квадратів вибрані групи експериментальних точок. Оскільки найкраще згладжування отримується для середніх точок (коли враховується інформація про поведінку функції по обидва боки від згладжуваної точки), то кількість точок для згладжування вибирають непарну, а групи точок – ковзними вздовж всієї таблиці: беруть, наприклад, перші п’ять точок у1, у2, у3, у4, у5 і згладжують за їх допомогою середню точку у3, потім беруть наступну групу у2, у3, у4, у5, у6 і згладжують точку у4, і так далі до кінця таблиці. При цьому залишаються декілька кратних точок (наприклад, дві на початку і дві в кінці таблиці), які згладжуються з меншою точністю.
Нижче наводяться найбільш вживані з простих формул згладжування для таблиць з постійним кроком. В цих формулах прийняті такі позначення. Середній точці групи приписується індекс 0, симетричні точки отримують при цьому індекси ±1, ±2,...
Згладжувані значення позначаються хвилькою зверху. Основною формулою слугує формула згладжування середньої точки, тобто формула для , інші формули використовуються тільки на краях таблиці.
2.3.2 Лінійне згладжування
Лінійним згладжуванням називається згладжування многочленом першого степеня.
Формули лінійного згладжування за трьома точками:
Формули лінійного згладжування за п’ятьма точками:
2.3.3 Нелінійне згладжування
З множини формул нелінійного згладжування наведемо тільки формули згладжування многочленів третього степеня за сімома точками:
Формули згладжування многочленами більш високих степенів майже не використовуються, а формули згладжування за більшим числом точок використовуються досить рідко, тому що вони залишають погано згладженими досить велику кількість точок по краях таблиці.
Зауважимо, що результат згладжування сильно залежить від використаних формул, вибір яких підказується швидше інтуїцією, ніж якими-небудь правилами.
2.3.4 Лінійна кореляція і регресія
При вивченні залежності між двома величинами, кожна з яких піддається випадковому розсіюванню (безконтрольному розподілу), використовуються методи кореляційного аналізу. Кореляційний аналіз вивчає усереднений закон поведінки кожної з величин в залежності від значень іншої величини, а також міру залежності між розглянутими величинами. Підставляючи кожному значенню однієї величини, наприклад, х, середнє з відповідних значень другої величини, наприклад, у, ми отримуємо функцію регресії або просто регресію у на х. Аналогічно визначається регресія величини х на величину у. Функція регресії зображається графічно лінією регресії. Міра залежності між величинами характеризується коефіцієнтом кореляції або кореляційними відношеннями.
Кореляція між величинами х та у називається лінійною, якщо обидві функції регресії лінійні. В цьому випадку лінії регресії перетворюються в прямі регресії. Кутові коефіцієнти цих прямих виражаються через коефіцієнт кореляції, який служить також мірою лінійної залежності між величинами.
Коефіцієнтом кореляції між випадковими величинами х та у називається математичне сподівання добутку їх нормованих відхилень:
,
де і – центри розподілення величин х і у, і – їх дисперсії.
Коефіцієнт кореляції може бути також записаний в одній з таких форм:
Величина називається кореляційним моментом.
Коефіцієнт кореляції – це безрозмірна величина, його абсолютна величина не перевищує одиницю: .
Для незалежних величин х і у коефіцієнт кореляції дорівнює нулеві. Якщо для залежних величин коефіцієнт кореляції дорівнює нулю, то в цьому випадку величини х і у називаються некорельованими.
Коефіцієнт кореляції характеризує міру лінійної залежності між величинами х і у. Пояснимо це. Нехай – лінійна функція середньоквадратичного наближення до величини у, тобто така лінійна функція, для якої математичне сподівання досягає найменшого значення.
Для експериментального вивчення залежності між двома величинами х і у проводять деяку кількість п незалежних випробувань (дослідів, спостережень). Результат і-го дослідження дає пару значень хі, уі . За цими значеннями визначаються точкові оцінки як середніх значень, так і коефіцієнта кореляції.
Незмінними і самостійними оцінками теоретичних середніх значень a і b служать емпіричні середні значення:
(2.3.3)
Незміщеними і самостійними оцінками дисперсії і служать емпіричні дисперсії:
(2.3.4)
Для оцінки кореляційного моменту служить емпіричний кореляційний момент:
За цими оцінками будують емпіричний коефіцієнт кореляції:
(2.3.6)
Зауважимо, що емпіричний коефіцієнт кореляції, як і теоретичний, за абсолютною величиною не перевищує одиниці.
Коефіцієнт r не змінюється при зміні початку відліку і масштабу вимірювання величин х і у. Ця властивість дозволяє суттєво спростити обчислення за допомогою вибору зручного початку відліку (хо, уо) і одиниць масштабу.
Після заміни:
,
тобто , емпіричний коефіцієнт кореляції обчислюється за формулою:
, (2.3.7)
де .
Приклад розрахунку коефіцієнта кореляції
Проведено вимірювань значень х і у. В перших трьох стовпцях таблиці знаходяться різні значення вимірювань (хі, уі) з кількістю вимірювань ті , в яких зустрічалась кожна точка (хі, уі). Для величини х вибрано початок відліку , масштабний коефіцієнт . Для величини у, відповідно, і . Таким чином, величини наближають тільки цілі значення. При підрахунку сум, які входять в формулу (2.3.4), кожний доданок враховується стільки раз, скільки він зустрічається в таблиці, тому для розрахунку введені стовпці і . В другому рядку умовно вказаний порядок дій. В останньому рядку підраховані суми, які необхідні для розрахунків:
Зауважимо, що тут суми результатів беруться не по всіх п = 26 вимірюваннях, як в формулі (2.3.7), а лише по різних точках (хі, уі), що і викликає необхідність врахування множників ті. Наприклад, повна сума тут записується у вигляді , а середнє значення – .
Для контролю розрахунків рекомендується повторити розрахунок з іншим початком відліку.
За отриманими даними знаходимо середнє значення і емпіричний коефіцієнт кореляції:
.
Таблиця 2.1
|
х |
у |
т |
u |
um |
U2m |
||||
|
(1) |
(2) |
(3) |
(4) |
(5)=(3)(4) |
(6)=(4)(5) |
(7) |
(8)=(3)(7) |
(9)=(7)(8) |
(10)=(5)(7) |
|
23,0 |
0,48 |
2 |
-6 |
-12 |
72 |
-2 |
-4 |
8 |
24 |
|
24,0 |
0,50 |
4 |
-4 |
-16 |
64 |
0 |
0 |
0 |
0 |
|
24,5 |
0,49 |
3 |
-3 |
-9 |
27 |
-1 |
-3 |
3 |
9 |
|
24,5 |
0,50 |
2 |
-3 |
-6 |
18 |
0 |
0 |
0 |
0 |
|
25,0 |
0,51 |
1 |
-2 |
-2 |
4 |
1 |
1 |
1 |
-2 |
|
25,5 |
0,52 |
1 |
-1 |
-1 |
1 |
2 |
2 |
4 |
-2 |
|
26,0 |
0,49 |
2 |
0 |
0 |
0 |
-1 |
-2 |
2 |
0 |
|
26,0 |
0,51 |
1 |
0 |
0 |
0 |
1 |
1 |
1 |
0 |
|
26,0 |
0,53 |
2 |
0 |
0 |
0 |
3 |
6 |
18 |
0 |
|
26,5 |
0,50 |
1 |
1 |
1 |
1 |
0 |
0 |
0 |
0 |
|
26,5 |
0,52 |
1 |
1 |
1 |
1 |
2 |
2 |
4 |
2 |
|
27,0 |
0,54 |
2 |
2 |
4 |
8 |
4 |
8 |
32 |
16 |
|
27,0 |
0,52 |
1 |
2 |
2 |
4 |
2 |
2 |
4 |
4 |
|
28,0 |
0,53 |
3 |
4 |
12 |
48 |
3 |
9 |
27 |
36 |
|
Суми |
26 |
- |
-26 |
248 |
- |
22 |
104 |
87 |
2.3.5 Рівняння прямих регресії
Пряма регресії у на х має рівняння:
, (2.3.8)
де – параметри.
Пряма регресії х на у має рівняння:
(2.3.9)
Коефіцієнти
називаються коефіцієнтами регресії у на х та х на у, відповідно. Обидва коефіцієнти регресії і мають той же знак, що і коефіцієнт кореляції . Обидві прямі регресії проходять через точку (a, b) – центр розподілення. Прямі регресії у на х та х на у збігаються тільки тоді, коли , тобто у випадку лінійної функціональної залежності між величинами х та у.
За результатами експерименту параметри рівнянь (2.3.8) і (2.3.9) оцінюються за допомогою формул (2.3.3, 2.3.4, 2.3.6).
Ці оцінки визначають емпіричні прямі регресії. Емпірична пряма регресії у на х має рівняння:
, (2.3.10)
причому
є емпіричним коефіцієнтом регресії у на х.
Параметри лінійної функції (2.3.10) задовольняють принцип найменших квадратів за у: сума квадратів відхилень спостережуваних значень уі від розрахункових за рівнянням прямої регресії (2.3.10) менше, ніж сума квадратів відхилень їх від будь-якої іншої прямої. Таким чином, має місце нерівність:
,
або, з врахуванням (2.3.7), маємо:
Емпірична пряма регресії х на у має рівняння:
,
де
є емпіричним коефіцієнтом регресії х на у. Для цієї прямої принцип найменших квадратів виконується відносно величини х: сума квадратів відхилень спостережуваних значень хі від розрахункових за рівнянням прямої регресії менше, ніж сума квадратів відхилень їх від будь-якої іншої прямої. Таким чином, має місце нерівність:
.
При цьому .
На рис. 2.6 i рис. 2.7 показано, які відхилення маються на увазі при розгляді прямих регресії у на х і х на у.
Рисунок 2.6 Рисунок 2.7
Обидві емпіричні прямі регресії проходять через центр емпіричного розподілення – точку . Пряма регресії х на у завжди проходить крутіше (з більшим нахилом до осі х), ніж пряма регресії у на х.
Емпіричні коефіцієнти регресії і мають той же знак, що і емпіричний коефіцієнт кореляції , і пов’язані співвідношенням .
Приклад Визначимо рівняння емпіричних прямих регресії попереднього прикладу. Маємо:
Емпіричний коефіцієнт регресії у на х:
.
Емпіричний коефіцієнт регресії х на у:
;
Рисунок 2.8
2.4 Інтерполяція функцій однієї і декількох змінних
Повним описом функціональної залежності відповідних результатів випробувань (вимірів) вважають задання формули або нескінченний набір чисел. Розв’язання задач на ЕОМ відбувається тільки з кінцевою сукупністю чисел. Виникає необхідність охарактеризувати функцію кінцевим набором чисел.
За методом сіток функції описують їх значеннями в кінцевому числі k–x точок. Це означає, що на деякому інтервалі [a,b] у вузлах сітки маємо відповідні значення функції
2.4.1 Одновимірна інтерполяція
На практиці найбільш часто необхідно відновити функцію за її сітковим поданням. Це роблять за допомогою інтерполяції.
Задача одновимірної інтерполяції полягає в побудові такої функції, що f(xk)=yk для всіх k і при цьому f(x) повинна приймати розумні значення для X, які містяться між заданими точками. Критерій розумності для кожної задачі свій і, напевне, ніколи точного визначення не буде мати.
Найбільш важливим класом інтерполяційних функцій є множина алгебраїчних поліномів. Задача поліноміальної інтерполяції полягає в побудові полінома
степені n-1, який приймає у вузлах x1, x2,…,xn значення f1, f2,…, fn. З курсу математичного аналізу відомо, що будь-яка неперервна функція може бути представлена на кінцевому інтервалі [a,b] поліном Рn степені n, а через n точок можна провести єдиний поліном степені n-1.
Відновлення коефіцієнтів полінома Рn-1(x;f) можна провести, розв’язуючи систему n рівнянь відносно невідомих лінійних коефіцієнтів аk (k=0, n-1). Але в багатьох випадках рівняння недостатньо обумовлені. Найбільш задовільний спосіб визначення полінома, який інтерполює точки (xk, fk), полягає у використанні базису лагранжевих поліномів на множині (xk). Це такі поліноми lj(x) (j=0, n-1) степені n-1, що
Ці умови задовольняє поліном
. (2.4.1)
Кожний множник чисельника (2.4.1) обертається в нуль при . Відповідні множники знаменника нормують поліномом так, що . Поліном приймає значення в точці і дорівнює нулеві у всіх точках . Таким чином, інтерполяційний поліном степені , який проходить через n точок (хk, fk), виражається формулою
(2.4.2)
і називається інтерполяційним поліномом Лагранжа.
Похибка поліноміальної інтерполяції на інтервалі [а, в] оцінюється за формулою:
де - похідна порядку n.
Якщо похідні високих порядків великі за абсолютним значенням, то точність апроксимації, згідно з (2.4.3), буде малою, тому немає потреби вибирати велику кількість вузлів. Крім того, при підвищенні степені полінома (2.4.2) результат через погану обумовленість стає досить чутливим до помилок округлення. При недостатній функції відрізок а, в розбивають на часткові інтервали і на кожному з них використовують інтерполяцію невисокого порядку. Наприклад, якщо відомо, що має обмежену тільки другу похідну , то бажано обмежитись лінійною інтерполяцією. Тоді похибка на рівномірній сітці з кроком h буде не перевищувати М2h2/8. Таким чином, процес кусково-лінійної інтерполяції збігається при зі швидкістю h2, тоді як процес інтерполяції в цілому може розбігатись.
Для наближеного відновлення функції за її значеннями на сітці використовують, крім поліномів, інші функції. Наприклад, для подання періодичних функцій використовують тригонометричну інтерполяцію.
Приклад
За дослідними даними залежність водопоглинання w від температури опалення T становить
Т, ºС 1000 1200 1300
w, % 15 7 2 ,
визначити водопоглинання для проміжного значення Т=1250ºС.
Розв’язання: Інтерполяційний поліном Лагранжа для триточкової інтерполяції такий:
. (2.4.4)
Якщо підставити значення , х1 = 1000, х2 = 1200, х3 = 1300, , і , отримаємо:
.
Отже, при Т=1250ºС водопоглинання w=4,6%.
2.4.2 Інтерполяція сплайнами
. (2.4.5)
Будемо шукати кубічний поліном у вигляді:
(2.4.6)
З умови маємо:
(2.4.7)
Визначимо похідні:
і запишемо умову їх неперервності при :
(2.4.8)
Загальна кількість невідомих коефіцієнтів дорівнює 4n, але кількість рівнянь (2.4.7) і (2.4.8) дорівнює 4n-2. Інші два рівняння отримаємо з умови (2.4.5) при і :
Виразивши з (2.4.8) та підставляючи його в (2.4.7) з врахуванням , отримаємо:
(2.4.9)
З першого рівняння (2.4.8) і (2.4.9) та виразу di отримуємо для визначення сі різницеве рівняння другого порядку:
(2.4.10)
з граничними умовами:
(2.4.11)
Різницеве рівняння (2.4.10) з умовами (2.4.11) розв’язується методом прогонки.
Використаємо кубічну сплайн-інтерполяцію для розв’язання попереднього прикладу.
Оскільки необхідно визначити значення , то кубічний поліном будемо шукати у вигляді
(2.4.12)
на рівномірній сітці з кроком на інтервалі .
Коефіцієнт .
Різницеве рівняння (2.4.10) запишеться у вигляді:
.
З врахуванням (2.4.11) , тоді . Звідки . Тоді .
Визначимо за формулою (2.4.9) значення в2:
.
В результаті кубічний поліном набуде вигляду:
Звідки
Інтерполяція використовується також в задачі оберненого інтерполювання: задана таблиця ; знайти хі як функцію від уі. Прикладом оберненого інтерполювання може служити задача про знаходження коренів рівняння.
Інтерполяційні формули використовуються також при обчисленні інтегралів, при написанні різницевих апроксимацій для диференціальних рівнянь на основі інтегральних тотожностей.
В різних додатках широко використовують двовимірні, тривимірні таблиці. Наприклад, теплофізичні властивості різноманітних речовин залежать від температури і тиску, а оптичні характеристики – ще і від довжини хвилі випромінювання.
При багатовимірній інтерполяції через громіздкість таблиць необхідно брати достатньо великі кроки по аргументах, тобто сітка вузлів, на якій будують таблицю, отримується достатньо грубою. Тому необхідно вводити перетворення змінних , , , підбираючи відповідні формули. При вдалому виборі таких формул можна використовувати в нових змінних інтерполяційний поліном невисокого степеня.
Здійснюючи багатовимірну інтерполяцією, необхідно пам’ятати, що розміщення вузлів не може бути довільним. Наприклад, при інтерполяції полінома першого степеня вузли не повинні лежати на одній прямій в площині. Дійсно, визначник системи трьох рівнянь
(2.4.13)
записується у вигляді:
. (2.4.14)
Умови розміщення трьох точок на одній прямій виглядають таким чином:
.
Після простих перетворень маємо:
, (2.4.15)
тобто, якщо вузли лежать на одній прямій, то визначник дорівнює нулю і побудувати поліном вигляду (2.4.15) неможливо. Перевіряти умови потрібного типу достатньо складно, тому на практиці намагаються будувати регулярні сітки, як правило, прямокутні і рівномірні, коли вузли є точками перетину двох взаємно перпендикулярних систем паралельних прямих. На цій сітці проводять просту послідовну інтерполяцію спочатку за рядками, а потім за стовпцями.
При послідовній інтерполяції завищується степінь інтерполяційного полінома. При трикутній конфігурації розміщення вузлів степінь многочлена буде мінімальною. Многочлен n-ої степеня у формі Ньютона в даному випадку можна уявити як узагальнення одновимірного варіанта запису:
. (2.4.16)
Приклад
Записати многочлен Ньютона першого і другого степеня для двовимірної інтерполяції функції .
З (2.4.16) отримаємо:
В деяких випадках використовують нерегулярні сітки. Тоді обмежуються інтерполяційним поліномом першого степеня і його коефіцієнти знаходять за трьома вузлами:
Коефіцієнти а, в, с не обов’язково визначати, оскільки перший стовпець є лінійною комбінацією трьох інших стовпців. В результаті з цих стовпців можна скласти визначник, який буде дорівнювати нулю. Розглядаючи цей визначник за першим стовпцем, отримаємо залежність .
ІІІ Елементи чисельних методів
3.1 Розв’язання системи лінійних алгебраїчних рівнянь
методом простої ітерації
3.1.1 загальні положення
Розглядається система n лінійних рівнянь з n невідомими
(3.1.1)
де
для застосування методу простої ітерації використовується зведена система
(3.1.2)
де
Перехід від системи (3.1.1) до зведеної системи (3.1.2) здійснюється таким чином:
(3.1.3)
(3.1.4)
Довільне (k-1)-e наближення обчислюється за формулою:
(3.1.5)
або в розгорнутому вигляді
(3.1.6)
(3.1.7)
Процес ітерації (3.1.5) стосовно зведеної системи (3.1.2) збігається незалежно від початкового наближення, якщо будь-яка канонічна норма матриці a менше 1,тобто ½a½<1.
Оцінка похибки методу простої ітерації здійснюється за формулою:
. (3.1.8)
Якщо необхідно обчислити значення коренів з точністю e ,то, як тільки D £ e, виконується умова ½x - x (k+1)½< e, тобто x (k+1) – шукане рішення.
Алгоритм методу:
3.1.2 Приклад
Знайти розв’язок системи алгебраїчних рівнянь за методом простої ітерації з точністю e=10-5, де
Переходимо до зведеної системи, яка має вигляд (3.1.2):
Звідки:
Перевіряємо умову збіжності ½a½<1. Знайдемо, наприклад, l-норму матриці a:
Маємо:
Отже, процес ітерації для даної системи є збіжним.
За початкове наближення візьмемо
Далі процес обчислення проводимо у відповідності з наведеним алгоритмом до досягнення заданої точності.
3.1.3 Текст програми розрахунку
10 REM МЕТОД ІТЕРАЦІЙ
11 DIM X(4) , A(4, 4) , XI(4) , B(4) , XP(4) , AL(4, 4) , BT(4)
12 INPUT E
20 PRINT " ВВЕДІТЬ СТОВПЕЦЬ B(I) "
25 FOR I=1 TO 4
30 INPUT B(I)
40 NEXT I
50 FOR I=1 TO 4 : FOR J=1 TO 4
60 PRINT "A("I" ,"J")="
70 INPUT A(I ,J) : NEXT J : NEXT I
71 FOR I=1 TO 4 : BT(I)=B(I)/A(1, 1) : XI(I)=BT(I)
72 RI=RI+XI(I)
80 NEXT I
90 FOR I=1 TO 4 : FOR J=1 TO 4
92 AL(I, J)=-A(I, J)/A(I, I)
93 IF I=J THEN AL (I, J)=0
95 XP(I)=XP(I)+AL(I, J)*XI(J)
96NEXT J : NEXT I
97FOR I=1 TO 4
100 X(I)=XP(I)+BT(I)
101 RR=RR+X(I)
104 XI(I)=X(I) : XP(I)=0 : B(I)=0
105 NEXT I
110 IF ABS (RI-RR) < E GOTO 133
111 RI=RR
120 RR=0
130 GOTO 90
133 FOR I=1 TO 4 : PRINT "X("I")=" ; X(I) , "XI("I")=" ; XI(I) : NEXT I
136 REM ПЕРЕВІРКА РОЗРАХУНКУ
140 I=1 TO 4 : FOR J=1 TO 4 : B(I)=B(I)+A(I, J)*X(J) : NEXT J : LPRINT "B("I")=" ; B(I) : NEXT I
150 END
В програмі: стовпець В(І)→`; елементи матриці А(I, J)→ аij; BT(I)→ ві/ij; AL(I, J)→ ij; XI(I)→xi(0); XP(I)→xi(K+1); RI→X; RR→X(K+1).
3.2 Обчислення дійсного кореня рівняння з даною точністю із
застосуванням комбінованого методу
3.2.1 Загальні положення
Розглядається рівняння виду:
(3.2.1)
де функція визначена та неперервна на якомусь кінцевому або нескінченному інтервалі.
Процес знаходження ізольованих дійсних коренів рівняння (3.2.1) складається з двох етапів: відділення коренів та обчислення кореня рівняння з даною точністю.
Комбінований метод хорд та дотичних для обчислення дійсного кореня рівняння (3.2.1) полягає в одночасному застосуванні методів хорд та дотичних, при цьому на кожному кроці здобуваємо значення кореня за відповідними формами.
Нехай корінь є відділеним на відрізку [a, b] та виконані умови, які гарантують єдиність кореня на [a, b] та збіжність методу: ; і – неперервні, збігаються і мають постійні знаки на відрізку [a, b].
Введемо такі позначення:
– значення наближеного кореня, здобутого за методом дотичних;
– значення кореня знаходять за формулами:
Початкове наближення для методу дотичних вибирають у відповідності з умовою:
(3.2.4)
За початкове наближення для методу хорд тоді приймається , або , відповідно.
Процес обчислення кореня закінчується, коли виконується умова:
(3.2.5)
де – дана точність.
За наближене значення кореня рівняння приймається:
(3.2.6)
Алгоритм методу містить такі етапи:
3.2.2 Приклад
Обчислити корінь рівняння з точністю до комбінованим методом.
Графічно відділяємо корені. Для цього задане рівняння записуємо у вигляді . Будуємо графік функцій та (рис.3.1).
Точний корінь рівняння . Другий точний корінь рівняння , відрізок – інтервал ізоляції кореня.
Перевіряємо умови, що гарантують єдиність кореня на і збіжність методу:
Рисунок 3.1
Тоді маємо:
неперервна на і не змінює знак
( неперервна на ).
неперервна на і не змінює знак
( неперервна на ).
У відповідності з умовою (3.2.4) за початкове наближення для методу дотичних беремо , для методу хорд ().
Далі процес обчислення кореня здійснюється відповідно з наведеним алгоритмом до досягнення заданої точності .
3.2.3 Текст програми розрахунку
10 REM МЕТОД ДОТИЧНИХ ТА ХОРД
20 PRINT " ВВЕДІТЬ ПОЧАТКОВІ ЗНАЧЕННЯ XP , XO , E "
25 INPUT XP , XO , E
30 DEF FNA (X)=X^2+4*SIN(X)
40 DEF FNP (X)=X*2+4*COS(X)
50 D=XP- XO
60 XN=XO-(FNA (XO) / (FNA (XP)-FNA (XO)))*D
70 XR=XP-FNA (XP) / FNP (XP)
80 IF ABS (XN-XR) < E GOTO 100
90 XP=XR : XO=XN : GOTO 50
100 X=(XN+XR) / 2
110 LPRINT " НАБЛ. КОРІНЬ= " ; X , " ЗНАЧ. F("X")=" ; FNA (X)
120 END
НАБЛ. КОРІНЬ = -1.933832 ЗНАЧ. F(-1.933832)=4.155636E-04
Для запуску програми необхідно уточнити:
FNA(X) - вираз функції f(x);
FNP(X) - вираз похідної від f(x);
XP - початкове наближене значення кореня за методом дотичних;
ХО - початкове наближене значення кореня за методом хорд;
3.3 Наближене розв’язання систем нелінійних
рівнянь методом Ньютона
3.3.1 Загальні положення
Розглядається система n лінійних рівнянь з n невідомими:
(3.3.1)
У матричному вигляді система (3.3.1) запишеться як:
(3.3.2)
Для розв’язання системи нелінійних рівнянь (3.3.2) застосовується метод послідовних наближень, в якому значення кореня на (k+1)-y кроці визначається таким чином:
(3.3.3)
де – вектор поправок.
Суть методу Ньютона полягає у тому, що розв’язання системи нелінійних рівнянь зводиться до знаходження послідовності розв’язків на кожному кроці системи n лінійних рівнянь відносно поправок`D(k), в якій матриця коефіцієнтів є матрицею Якобі:
(3.3.4)
де
(3.3.5)
. (3.3.6)
За нульове наближення беремо грубе значення кореня.
Процес обчислення кореня системи (3.3.2) продовжується до виконання однієї з умов:
У разі виконання (3.3.7) або (3.3.8) за корінь системи (3.3.2) приймаємо =(k+1).
Алгоритм методу Ньютона містить такі етапи:
1. Беремо k=0. Вводимо початкове наближення (k) та точність e.
2. Обчислюємо вектор правих частин системи лінійних рівнянь` f((k)) у відповідності до (3.3.5).
3. Обчислюємо матрицю Якобі W((k)) у відповідності до (3.3.6).
4. Розв’язуємо систему лінійних рівнянь (3.3.4).
5. Обчислюємо (k+1) у відповідності до (3.3.3).
6. Перевіряємо умову (3.3.7) або (3.3.8). Якщо умова виконана, беремо k=k+1 і переходимо до п.2, інакше – до п.7.
7. Виводимо (k+1), кількість ітерацій.
3.3.2 Приклад
Знайти розв’язок системи нелінійних рівнянь з точністю до ε =10-3.
Початкове наближення:
Складаємо матрицю Якобі:
Складаємо систему лінійних рівнянь вигляду (3.3.4) для початкового наближення.
Маємо:
Система лінійних рівнянь для k=0 має вигляд:
Розв’язуючи систему, маємо:.
Звідки
Перевіряємо виконання умови (3.3.7):
Звідки одержуємо:
Процес обчислення продовжується у відповідності з наведеним алгоритмом до досягнення точності ε.
3.3.3 Текст програми розрахунку
10 REM PROGRАM METON
20 PRINT " МЕТОД НЬЮТОНА"
30 PRINT " УКАЖІТЬ ПОРЯДОК НЕЛІНІЙНОЇ СИСТЕМИ – N, ТОЧНІСТЬ - E"
40 INPUT N, E
50 DIM A(N, N), B(N), X(N), Y(N)
60 PRINT " ПОЧАТКОВІ НАБЛИЖЕННЯ X(I) "
70 FOR I=1 TO N: PRINT " X("I")= " : INPUT X(I) : NEXT I
80 REM УТОЧНІТЬ ЕЛЕМЕНТИ МАТРИЦІ ЯКОБІ –А(I, J)
90 A(1, 1)=1+2*X(1)
100 A(1, 2)=-2*X(3)
110 A(1, 3)=-2*X(2)
120 A(2, 1)=3*X(3)
130 A(2, 2)=1-2*X(2)
140 A(2, 3)=3*X(1)
150 A(3, 1)=2*X(2)
160 A(3, 2)=2*X(1)
170 A(3, 3)=1+2*X(3)
180 REM УТОЧНІТЬ ВИРАЗИ ФУНКЦІЙ F(I) ДЛЯ B(I)
190 B(1)=-(X(1)+X(1)^2-2*X(2)*X(3)-.1)
200 B(2)=-(X(2)+X(2)^2-3*X(1)*X(3)+.1)
210 B(3)=-(X(3)+X(3)^2-2*X(1)*X(2)-.3)
220 REM РОЗВ’ЯЗОК СИСТЕМИ W(X(K)) Y(K) =-F(X(K))
230 FOR I=1 TO N=1 : FOR J=I+1 TO N
240 A(J, I)=-A(J, I) / A(I, I) : FOR K=I+1 TO N
250 A(J, K)=A(J, K)+A(J, I)*A(I, K) : NEXT J
260 B(J)=B(J)+A(J, I)*B(I) : NEXT J
270 Y(N)=B(N) / A(N, N)
280 NEXT I
290 FOR I=N-1 TO 1 STEP –1 : H=B(I)
300 FOR J=I+1 TO N : H=H-Y(J)*A(I, J) : NEXT J
310 Y(I)=H/A(I, I)
320 NEXT I
330 LPRINT " КОРЕНІ ДЛЯ ДЕЛЬТА - Y " : PRINT
340 FOR I=1 TO N : LPRINT " Y(";I;")=" Y(I) : NEXT I
350 FOR I=1 TO N : X(I)=X(I)+Y(I) : NEXT I
360 KI=KI+1
370 F= ABS (B(1)+B(2)+B(3))/3
380 IF F>E GOTO 90
390 LPRINT " КОРЕНІ СИСТЕМИ РІВНЯНЬ " : PRINT
400 FOR I=1 TO N : LPRINT "X("I")=" ; X(I) , "F("I")=" ; B(I)
410 NEXT I
420 LPRINT " ІТЕРАЦІЙ – КІ= " ; KI , " СЕРЕДНЄ АРИФМ. ЗНАЧЕНЬ Ф. F(I) " ; F
430 STOP
3.4 Точкове апроксимування функцій
за методом найменших квадратів
3.4.1 Загальні положення
Нехай функція y=f(x) задається таблично на множині точок x0,x1,…,xn, тобто, y0=f(x0), y1=f(x1),…, yn=f(xn).
Потрібно апроксимувати функцію f(x) поліномом
(3.4.1)
таким чином, щоб відхилення функції f(x) від полінома Q(x) на даній множині точок було мінімальним.
Мірою відхилення будемо вважати квадратичне відхилення:
(3.4.2)
Для побудови апроксимованого полінома потрібно вибрати його коефіцієнти a0,a1,…,an таким чином, щоб величина S була мінімальною.
В точці мінімуму функції її перші похідні обов’язково повинні дорівнювати нулю. Диференціюючи функцію (3.4.2) за змінними та прирівнюючи похідні dS/dai до нуля, отримуємо систему лінійних рівнянь.
Коефіцієнти полінома a0,a1,…,am є розв’язком системи лінійних рівнянь:
, (3.4.3)
де
(3.4.4)
Поліном (3.4.1) з коефіцієнтами a0,a1,…,am, що визначаються з системи (3.4.3), матиме мінімальне квадратичне відхилення від даної функції f(x).
3.4.2 Приклад
Побудувати апроксимувальний поліном для функції f(x), заданої таблично:
Апроксимувальний поліном матиме вигляд:
Для знаходження коефіцієнтів a1, a2, a3 складаємо систему вигляду:
(3.4.5)
Знайдемо коефіцієнти системи S0,…,S4 та вільні члени t0,t1,t2.
Таким чином, система має вигляд:
Розв’язуючи цю систему, маємо:
Апроксимувальний поліном має вигляд
Обчислюємо значення полінома Q(x) у точках xi:
Підрахуємо мінімальне квадратичне відхилення за формулою:
Побудуємо графіки функції f(x) та полінома Q(x) (рис.3.2).
3.4.3 Текст програми розрахунку
5 REM АПРОКСИМАЦІЯ ПОЛІНОМОМ 2-ГО ПОРЯДКУ
10 DIM S(6) , X(6) , Y(6) , T(6) , Q(6)
20 FOR I=0 TO 3
30 REM ВВЕДЕННЯ ПОЧАТКОВОГО НАБЛИЖЕННЯ
40 INPUT X(I) , Y(I)
50 NEXT I
60 FOR J=0 TO 4
65 FOR I=0 TO 3
70 S(J)=S(J)+X(I)J
80 NEXT I
90 NEXT J
100 FOR J=0 TO 2
110 FOR I=0 TO 3
120 T(J)=T(J)+X(I)J*Y(I)
130 NEXT I
131 LPRINT " T("J")= "; T(J)
140 NEXT J
150 A(1,1)=S(0)
160 A(2,1)=S(1)
170 A(1,2)=A(2,1)
180 A(3,1)=S(2)
190 A(2,2)=S(2)
200 A(1,3)=S(2)
210 A(3,2)=S(3)
220 A(2,3)=S(3)
230 A(3,3)=S(4)
240 FOR I=1 TO 3
250 FOR J=1 TO 3
260 LPRINT " A("I","J")= "; A(I,J)
270 NEXT J : NEXT I
271 PRINT "A0,A1,A2"
275 INPUT A0,A1,A2
276 FOR I=0 TO 3
277 Q(I)=Q(I)+A0+A1*X(I)+A2*X(I)2
278 D=(Y(I)-Q(I))2
279 LPRINT " Q("I")= "; Q(I) , " D= "; D
280 NEXT I
290 STOP
3.5 Одновимірна мінімізація функцій
3.5.1 Загальні положення
Розв’язується задача
(3.5.1)
Нехай функція y(x) унімодальна на відрізку [a,b], точка мінімуму x· належить [a,b]. Треба визначити точку мінімуму функції y(x) x· з заданою точністю ε.
Метод дихотомії. Нехай на k-му кроці точка мінімуму x· належить відрізку [ak,bk]. Відрізок [ak,bk] ділиться точкою (k) навпіл, потім обчислюються значення функції y(x) у точках
Порівнюючи значення y(x1(k)) та y(x2(k)), робимо висновок, до якого з відрізків [ak,(k)] чи [(k),br] належить точка мінімуму, тобто звужуємо відрізок [ak,bk]. Процес пошуку точки мінімуму продовжується до досягнення заданої точності.
Алгоритм методу дихотомії містить такі етапи:
5. Порівняти y(x1(k)) та y(x2(k)). Якщо y(x1(k)) ³ y(x2(k)), то ak+1=`(k), bk+1=bk, якщо y(x1(k)) < y(x2(k)), то ak+1=ak, bk+1= x(k).
6. Перевірити умови завершення пошуку: ½bk+1-ak+1½< ε. Якщо умова виконується, то перейти до наступного етапу, в іншому випадку – до етапу 3, припустивши k=k+1.
7. Взяти за точку мінімуму
Метод золотого перерізу. На кожному кроці методу точки x1(k),x2(k) знаходяться на відрізку [ak,bk] симетрично відносно його кінців і обчислюються за формулами, що витікають з пропорції золотого перерізу:
(3.5.2)
де
(3.5.3)
Порівняння значень функції y(x1(k)) та y(x2(k)) дає змогу визначити відрізок [ak,bk], який містить точку мінімуму x·.
Алгоритм методу золотого перерізу складається з таких етапів:
8. Вивести на друк або монітор x·, y(x·).
Знайти мінімум функції
на [0;0,5] з точністю до ε=0,05.
За допомогою методу дихотомії маємо:
1. a=0, b=0,5;
2. k=0; a0=a=0; b0=b=0,5;
4. y(x1(0))=-6,11; y(x2(0))=-5,453.
5. y(x1(0))<y(x2(0)) тому a1=a0=0; b1= (0)=0,25.
6.½b1-a1½>ε ®k=1 переходимо до пункту 3 і продовжуємо процес пошуку до досягнення потрібної точності.
За допомогою методу золотого перерізу маємо:
1. a=0, b=0,5;
2. k=0; a0=0; b0=0,5;
3. x1(0)=a0+0,382(b0-a0)=0,191; x2(0)=a0+0,618(b0-a0)=0,309;
5. y(x1(0))<y(x2(0)) ®
6.½b1-a1½>ε ® переходимо до пункту 5, припустивши k=1.
Друга ітерація
6.½b1-a1½>ε ® беремо k=2 і переходимо до пункту 5.
Процес пошуку продовжується до досягнення заданої точності.
3.5.3 Текст програми розрахунку
10 REM METOД ДИХОТОМІЇ
20 INPUT A , B , E
30 K=0
40 AK=A
50 BK=B
60 XG=(AK+BK) / 2
65 DEF FNF (X1)=-(803*X1^3)/3+275* X1^2-74*X1-1/3
66 DEF FNY (X2)=-(803*X2^3)/3+275* X2^2-74*X2-1/3
70 X1=XG+E / 2
80 X2=XG+E / 2
103 F1=FNF (X1)
104 F2=FNY (X2)
110 IF F1>=F2 GOTO 180
120 AK=AK
130 BK=XG
140 IF ABS (BK-AK) < G GOTO 160
150 GOTO 60
160 XX=(AK+DR) / 2
165 X1=XX
170 PRINT " ТОЧКА МІНІМ.-XX= " ; XX , "F("XX")=" ; F1
175 STOP
176 GOTO 210
180 AK=XG
190 BK=BK
200 GOTO 140
210 REM МЕТОД ЗОЛОТОГО ПЕРЕРІЗУ
220 K1=.382
230 K2=.618
240 AK=A
250 BK=B
260 X1=AK+K1*(BK-AK)
270 X2=AK+K2*(BK-AK)
280 F1=FNF (X1)
290 F2=FNY (X2)
300 IF F1>=F2 GOTO 400
310 AK=AK
320 BK=X2
330 X2=X1
340 F2=F1
350 X1=AK=K1*(BK-AK)
360 F1=FNF (X1)
370 IF ABS (BK-AK) / 2>E GOTO 300
380 XX=(AK+BK) / 2
390 GOTO 165
400 AK=X1
410 BK=BK
420 X1=X2
430 F1=F2
440 X2=AK+K2*(BK-AK)
450 F2=FNY (X2)
460GOTO 370
3.6 багатовимірна безумовна мінімізація функцій.
Градієнтні методи
3.6.1 Загальні положення
Треба розв’язати задачу
тобто знайти точку мінімуму функції y(x) з точністю до ε.
За допомогою градієнтних методів послідовність наближень до точки мінімуму визначається у відповідності з виразом:
де l(k) – параметр, що регулює величину кроку, 0<l(k) <1.
Градієнтні методи відрізняються засобами вибору параметра l(k). Процес пошуку точки мінімуму припиняється, коли необхідні умови існування точки мінімуму виконуються із заданою точністю ε.
3.6.1.1 Градієнтний метод з дрібненням кроку
У цьому методі параметр l(k) вибирається таким чином, щоб виконувалася умова
процес пошуку точки мінімуму продовжується до виконання умови:
Алгоритм методу містить такі етапи:
6. Порівняти y((k)) та y((k+1)). Якщо y((k))<y((k+1)), перейти до п.7, якщо y((k))<y((k+1)), то зменшити l(k) вдвічі і перейти до п.4
7. Перевірити умову:
Якщо умова не виконується, то, взявши k=k+1, перейти до п.3, у іншому випадку перейти до наступного пункту.
8. За точку мінімуму взяти ·= (k+1).
9. Вивести результати на друк або монітор: ·, y(·).
3.6.1.2 Метод координатного спуску
Суть цього методу полягає в тому, що на кожному r-ому кроці k-ої ітерації змінюється лише одна змінна xr, для якої не виконується необхідна умова існування точки мінімуму:
Параметр lr(k) обчислюється за формулою:
Алгоритм методу містить такі етапи:
8. Перевірити виконання умови:
Якщо умова не виконується, то припустити k=k+1 та перейти до п.2. У іншому випадку перейти до п.9.
9. Взяти за точку мінімуму ·= (k+1).
10. Вивести на друк або монітор`·, y(· ).
3.6.1.3 Приклад
Знайти мінімум функції
з точністю до e = 0,01. Початкове наближення:
Градієнтний метод з дрібненням кроку
Обчислюємо y((0)) = 4,25.
Знаходимо похідні:
Обчислюємо:
Таким чином,
Нульова ітерація: k=0.
Беремо l(0)=1.
Знаходимо:
Обчислюємо y((1)) = 16,5.
Адже y((1))>y ((1)), дрібнимо параметр l(0).
Знаходимо:
Умова y((1)) < y ((1)) не виконується, знову дрібнимо l(0):
Тепер умова y((1))<y ((1)) виконана, тобто
Обчислюємо:
тобто, задану точність не досягнуто, переходимо до наступної ітерації.
Перша ітерація: k=1.
Беремо l(0)=1.
Знаходимо:
Оскільки y((1)) < y ((1)), беремо l(0)=0,5.
Обчислюємо знову:
Умова y((1)) < y ((1)) виконана, тобто
Знаходимо:
тобто, задану точність розв’язку не досягнуто, треба взяти k=2 i перейти до наступної ітерації.
Метод координатного спуску
Обчислюємо y((0))=4,25.
Знаходимо вирази для перших та других похідних:
Нульова ітерація: k=0.
І-й крок: x=1.
Обчислюємо:
Задану точність за змінною x1 не досягнуто, адже
Знаходимо значення x1(1):
Перевіряємо виконання умови:
Враховуючи, що змінну x1 вже обчислено (дорівнює –0,705),
задану точність за змінною x2 не досягнуто.
ІІ-й крок. Знаходимо x2(1):
Перевіримо, що після ітерації ми наблизимося до точки мінімуму:
Перевіримо умову:
Маємо:
Необхідна умова існування точки мінімуму не виконується, тому переходимо до першої ітерації.
Перша ітерація: k=k+1=1.
І-й крок
Знаходимо x1(2):
ІІ-й крок
Знаходимо x2(2):
Оскільки y((2))=3,022<y((1))=3,316, то в результаті ітерації ми наближаємося до точки мінімуму (значення функції змінюється).
Перевіряємо виконання необхідної умови існування точки мінімуму:
Умова
не виконується, переходимо до наступної ітерації.
Процес пошуку мінімуму функції y() продовжується до досягнення заданої точності ε.
3.6.1.4 Текст програми розрахунку
10 REM ГРАДІЄНТНИЙ МЕТОД З ДРІБНЕННЯМ КРОКУ
20 DEF FNXX (X1, X2)=X1*X2-1/X1-1/X2
30 DEF FNXA (X1, X2)=X2-1/X1^2
40 DEF FNXB (X1, X2)=X1-1/X2^2
45 E=.001 : L=1
50 X1=-.8
60 X2=-.8
65 FS=FNXX (X1, X2)
70 F1=FNXA (X1, X2)
75 IF R=1 GOTO 230
80 F2=FNXB (X1, X2)
90 IF N=2 GOTO 310
95 GR= SQR(F1^2+F2^2)
100 IF GR<E GOTO 180
101 K=K+1
110 X1=X1-L*F2
120 X2=X2-L*F2
130 FH=FNXX (X1, X2)
140 IF FH>FS GOTO 160
150 GOTO 65
160 L=L/2
170 GOTO 65
180 PRINT " X1= " ; X1, " X2=" ; X2, " FH="; FH , " ЧИСЛО ІТЕР.= " ; K
185 STOP
190 REM МЕТОД КООРДИНАТНОГО СПУСКУ
200 DIM F(2)
210 DEF FNDD (X1)=-2/X1^3
220 DEF FNDE (X2)=-2/X2^3
221 N=2 : R=1 : K=0 : GOTO 50
230 IF ABS (F1)<E GOTO 290
240 F(1)=FNDD(X1)
250 L=1/F(1)
260 X1=X1-F1*L
270 FF=FNXX (X1, X2)
271 K=K+1
280 GOTO 65
290 R=N
300 GOTO 80
310 IF ABS (F2)<E GOTO 350
315 F(2)=FNDE (X2)
320 L=1/F(2)
330 X2=X2-F2*L
331 K=K+1
340 GOTO 80
350GOTO 180
диференціювання функцій однієї змінної
Задача чисельного інтегрування полягає в знаходженні наближеного інтеграла
,
де – задана функція.
На відрізку вводиться сітка: , і за наближене значення інтеграла розглядається число
, (3.7.1)
де – значення функції у вузлах , – деякі числа (вагові множники), які залежать тільки від вузлів, але не залежать від вибору . Формула (3.7.1) називається квадратурною формулою.
Найчастіше використовуються в чисельному інтегруванні квадратурні формули прямокутників, трапецій і Сімпсона.
Розглянемо найпростіші квадратурні формули.
Формула прямокутників:
.
Формула трапецій:
.
Формула Сімпсона є наслідком інтерполяції на інтервалі . Нехай числова область складається з двох інтервалів довжиною h: , . Тоді, замінюючи криву АВС інтерполяційною параболою, яка проведена через точки , , , отримують формулу Сімпсона (рис.3.3):
.
Рисунок 3.3
На відрізку а, в при формулу Сімпсона можна подати в такому вигляді:
де
При грубому оцінюванні інтеграла від остаточного члену квадратичної інтерполяції похибка формули складає .
Приклад
За дослідними даними про залежність ентропії від температури процесу Т
|
T(S), K |
278 |
303 |
323 |
340 |
350 |
354 |
347 |
330 |
325 |
|
S, кДж/(кгК) |
0,1 |
0,2 |
0,3 |
0,4 |
0,5 |
0,6 |
0,7 |
0,8 |
0,9 |
визначити середньоінтегральну температуру газу Тс.і., використовуючи для цього термодинамічне співвідношення:
.
Інтеграл в даній формулі обчислюється за методом Сімпсона:
Тоді:
Чисельне диференціювання або визначення похідних є теж не менш важливою задачею в прикладному аналізі.
Нехай відомі значення функції в точках , близьких до точки х (точка х може збігатися з однією з точок ). Необхідно знати наближене значення похідної та .
Один із способів розв’язання цієї задачі полягає в заміні функції її інтерполяційним поліномом:
В іншому способі – невизначених коефіцієнтів – використовується формула розкладання в ряд Тейлора функції , неперервної в околі точки :
Звідки
Оскільки , то похідна
(3.7.2)
називається, відповідно, правою і лівою похідною. Похибка обчислення похідної буде порядку О(h).
Формула (3.7.2) може бути покращена за рахунок триточкової апроксимації, яка має другий порядок точності. На підставі формули Тейлора:
;
.
Помножимо перший рядок на 4 і віднімемо його від другого, отримаємо:
де - похибка обчислення похідної.
Отримана формула може бути використана для визначення похідної в крайній лівій точці, а для останньої (кінцевої) точки використовується формула:
.
Найбільш практичною формулою для внутрішніх точок є центральна різницева похідна:
Визначимо симетричну апроксимацію другої похідної функції .
Нехай На підставі формули Тейлора маємо:
;
.
Складемо ці рівняння і визначимо :
.
Отримані формули мають здатність примножувати ефекти завад як при надмірно малому кроці сітки h, так і при великих його значеннях. Отже, в кожному випадку необхідний підбір оптимального кроку для максимально можливого уникнення завад.
Так, при розрахунку похідних на ЕОМ величину h (практичні результати) вибирають в межах 0,001 або 0,005. Покажемо це на наступному прикладі.
Визначимо першу та другу похідні функції в точці х = 0 при h = 0,001; 0,005 та порівняємо з аналітичними розрахунками:
;; ; .
При h = 0,001:
При h = 0,005:
Отже, правомірність наших передбачень очевидна.
Приклад
Використовуючи дані про залежність границі міцності R від водоцементного відношення В/Ц
|
В/Ц |
0,3 |
0,35 |
0,4 |
0,45 |
0,5 |
0,55 |
0,6 |
0,65 |
0,7 |
0,75 |
|
R, МПа |
32 |
36 |
41 |
41 |
38 |
36 |
34 |
32 |
30 |
28 |
оцінити за значенням чутливість R до зміни В/Ц при значеннях останнього: 0,3; 0,45 і 0,6.
Підстановкою відповідних значень в триточкову формулу для отримаємо:
Отже, в великій частині розглянутого діапазону В/Ц з його збільшенням міцність бетону зменшується, причому при великих В/Ц такого спаду менше (рис.3.4).
Рисунок 3.4
IV ЕЛЕМЕНТИ ВАРІАЦІЙНОГО ЧИСЛЕННЯ
ТА МЕТОД КІНЦЕВИХ ЕЛЕМЕНТІВ
Розглянемо деякі поняття варіаційного числення [13].
Основна теорема математичного аналізу встановлює, що будь-яка функція y=f(x), неперервна в області , досягає своїх мінімальних і максимальних значень в цій області. Тобто виконується рівність:
, (4.1)
з якої можна визначити значення х0, де f(x) досягає свого оптимального значення.
Твердження (4.1) може бути сформульоване в такій формі. Якщо функція у=f(х), неперервна в області , досягає оптимального значення в точці х=х0, то похідна df/dx в цій точці повинна дорівнювати нулеві. Отже, перша варіація функції - δf, яка іноді означається як диференціал df, дорівнює нулеві для будь-якої змінної δх в х, тобто
δf=(df/dx)δx=0. (4.2)
подальше диференціювання функції y=f(x) визначає відносний мінімум, максимум або мінімаксну умову функції f(x) в точці х=х0 :
мінімум, якщо d2f/dx2>0 при х=х0 , .
максимум, якщо d2f/dx2<0 при х=х0 , (4.3)
мінімакс, якщо d2f/dx2=0 при х=х0 . .
Рівняння (4.1), (4.2) дають необхідні умови існування екстремуму функції f(x) в точці х0, але вирази (4.3) вказують на недостатність умови. Рівняння (4.2) є необхідною і достатньою умовою стаціонарності функції f(x) в точці х=х0, тобто для f(x) виконуються умови (4.3).
Розглянемо умову стаціонарності функції двох незалежних змінних z=f(x,y). Якщо функція z стаціонарна в точці (х0,у0), то х0 і у0 будуть розв’язками системи рівнянь:
(4.4)
(4.5)
З іншого боку, якщо z= f(x,у) стаціонарна в точці (х0,у0), то варіація цієї функції повинна бути нульовою для будь-яких варіацій δх і δу, тобто
(4.6)
в точці (х,у)=( х0,у0). Рівняння (4.1) і (4.2) справедливі через довільність δх і δу.
Формулу (4.6) можна узагальнити і для функції n незалежних змінних: f(x1, x2, …, xn). Якщо існує точка (x1, x2,..., xn) = (x10, x20,..., xn0)= Х0, в якій функція f(xi) стаціонарна, то її необхідною і достатньою умовою стаціонарності є рівність нулеві в цій точці першої варіації для довільних варіацій δхі:
і=, тобто , (4.7)
де – частинна похідна від f(xі) на хі, яка вирахувана в точці Х0.
Рівняння (4.7) справедливе для будь-яких варіацій Нехай для незалежної змінної хк вибрана нульова варіація (). Тоді (4.7) запишеться у вигляді:
. (4.8)
Звідси маємо, що:
. (4.9)
Оскільки незалежна змінна хк була вибрана довільно, то вираз (4.7) справедливий для всіх . Тому необхідна і достатня умова стаціонарності функції f(xi) в точці Х0 зводиться до системи сумісних рівнянь:
(4.10)
4.1 Метод скінченних елементів
4.1.1 Загальні положення
МСЕ містить основні концепції методу сіток, пов’язані з дискредитацією областей неперервної зміни методу Гальоркіна, в яких використовуються варіаційні принципи для знаходження функцій у вузлах. В результаті розв’язання краєвої задачі зводиться до розв’язання системи лінійних алгебраїчних рівнянь (порядок яких визначається числом вузлів).
Основна ідея МСЕ полягає в тому, що неперервну величину (температуру, тиск і переміщення) можна апроксимувати дискретною моделлю, яка будується на множині кусково-неперервних підобластей. Кусково-неперервні функції визначаються за допомогою значень неперервної величини в кінцевому числі точок розглядуваної області.
Для побудови дискретної моделі неперервної величини чинять так.
Розбити область на елементи можна різними способами. Якщо обмежити кожний елемент двома сусідніми вузловими точками, то будуть створені чотири елементи (рис.4.1,а), або розбити область на два елементи, кожний з яких містить три вузли (рис.4.1,б). Відповідний поліном для кожного елемента визначають за значеннями Т(х) у вузлових точках елемента. Так, для елементів (рис.4.1,а) функції будуть лінійні за x, а для елементів з трьома вузловими точками функції будуть у вигляді полінома другого степеня (рис.4.1,б).
Рисунок 4.1,а Рисунок 4.1,б
В цьому випадку апроксимацією Т(х) є сукупність двох кусково-неперервних квадратичних функцій.
Якщо значення функції Т(х) у вузлових точках невідомі, то визначають множину вузлів і значення температури у вузлах Т1, Т2, Т3,..., які тепер є змінними, оскільки явно невідомі. Область розбивають на елементи, на кожну з яких визначають відповідну функцію елемента. Вузлові значення Т(х) повинні забезпечувати найкраще наближення до справжнього розподілення функцій Т(х). Для цього потрібно мінімізувати деяку величину, яка пов’язана з фізичною суттю задачі. Якщо розглядається задача розповсюдження, наприклад, теплоти, то мінімізується функціонал, пов’язаний з відповідними диференціальними функціями. Процес мінімізації зводиться до розв’язування системи лінійних алгебраїчних рівнянь відносно вузлових значень Т(х).
При побудові дискретної моделі неперервної величини, яка визначена в дво- або тривимірній області, основну ідею методу скінченних елементів використовують аналогічно.
В двовимірному випадку елементи описуються функціями х і у (трикутники або чотирикутники, тобто функція елемента подається площиною).
Якщо число вузлів, що використовується, більше мінімального, то функції елемента відповідає криволінійна поверхня. Крім того, надлишкове число вузлів дозволяє розглянути елементи з криволінійними границями. Кінцевою апроксимацією двовимірної неперервної величини φ(х,у) служить сукупність кусково-неперервних поверхностей, кожна з яких визначається на окремому елементі за допомогою значень φ(х,у) у відповідних вузлових точках.
4.1.3 Побудова інтерполяційних многочленів
а) Одновимірний лінійний елемент.
Елемент містить два вузли, невідому функцію записують у вигляді:
φ=α1+α2х. (4.1.1)
Коефіцієнти α виражаються через значення функцій у вузлах:
φ=NiФі+NjФj ,
де Ni , Nj – функції форми;
Фі ,Фj значення φ у вузлах i , j.
Ni=(xi-x)/l , Nj=(x-xj)/l ,
де xi , xj – координати вузлів елементів;
l – довжина елемента.
б) Двовимірний лінійний елемент, містить три вузли, шукана функція має вигляд:
(4.1.2)
Коефіцієнти виражаються через значення шуканої функції у вузлах (рис.4.2):
(4.1.3)
(оскільки площа трикутника ніколи не дорівнює нулеві).
Рисунок 4.2
(4.1.4)
Розв’язуючи систему (4.1.4) відносно 1, 2, 3, отримуємо такі вирази:
(4.1.5)
де
(4.1.6)
4.1.4 Одновимірний приклад варіаційного методу
скінченних елементів
Розглянемо одиничну масу, яка рухається під дією сили, що пропорційна пройденому шляху у (коефіцієнт пропорційності k = 1). Тоді диференціальне рівняння руху запишеться у вигляді:
. (4.1.7)
Необхідно знайти пройдений шлях за проміжок часу від t = 0 до t = 2, якщо дано:
(4.1.8)
Відомо, що розв’язок рівняння (4.1.7) можна отримати подвійним інтегруванням (4.1.7) з врахуванням початкових і кінцевих умов:
При t = 0, ,
t = 2, .
Звідки ,
Отже, у = еt.
Розв’язок рівняння (4.1.7) збігається з функцією, що мінімізує функціонал
(4.1.9)
за умови, що пробні функції неперервні, мають кусково-неперервні перші похідні і задовольняють умови (4.1.8).
Область [0≤t≤2] розіб’ємо на n кінцевих елементів е1, е2, ..., еn і у будемо апроксимувати пробною функцією в середині кожного елемента (рис.4.3).
Функціонал Х запишемо у вигляді суми елементів інтегралів:
. (4.1.10)
Визначаючи елементний внесок в формулу (4.1.10) через Хеі, можна записати:
. (4.1.11)
e1 e2 e3 en-1 en
t1=0 t2 t3 tn tn+1 t
Рисунок 4.3
Виберемо пробну функцію у вигляді:
(4.1.12)
де , - постійні, які можуть бути визначені із умов у вузлах:
(4.1.13)
Отримаємо:
(4.1.14)
З врахуванням (4.1.14) і (4.1.12) отримуємо елементну апроксимацію у вигляді:
(4.1.15)
Визначаючи коефіцієнти при і як базисні функції і , рівність (4.1.15) запишемо у вигляді апроксимації базисними функціями:
. (4.1.16)
Зауважимо, що і залежать тільки від t і зовні елемента пробна функція дорівнює нулеві.
Рівність (4.1.16) запишемо у матричній формі:
, (4.1.17)
де , . (4.1.18)
Величину називають матрицею базисних функцій, а - вузловим вектором.
Якщо елементи вибрані однакової довжини, тоді
. (4.1.19)
З врахуванням значень t з умови задачі рівність (4.1.19) запишеться у вигляді:
(4.1.20)
Підставимо (4.1.20) в (4.1.15), отримаємо елементарну пробну функцію:
(4.1.21)
Підстановка (4.1.21) в дає такий вираз елементного вкладу:
(4.1.22)
Проінтегруємо (4.1.22) з врахуванням (4.1.20), маємо:
. (4.1.23)
Підставляючи (4.1.23) в (4.1.11), приходимо до висновку, що Х є функцією вузлових значень , тобто
. (4.1.24)
З варіаційного числення відомо, що мінімумом для функціонала Х є:
(4.1.25)
Вузлові значення і замінимо постійними з врахуванням (4.1.8):
. Отже, в (4.1.25) р буде приймати значення 2, 3, ..., n. З рівності (4.1.23) визначимо :
. (4.1.26)
Аналогічно для елемента еі-1 маємо:
. (4.1.27)
Підставляючи (4.1.26) і (4.1.27) в (4.1.25), приходимо до висновку, що в рівняння (4.1.25) дають внесок тільки елементи еі, еі-1 , :
(4.1.28)
Для і = 1 та і = n +1 маємо:
(4.1.29)
(4.1.30)
Внесок в формули (4.1.29), (4.1.30) дають, відповідно, тільки елементи e1 і en.
Підставляючи рівняння (4.1.28) – (4.1.30) в (4.1.25), отримуємо систему n+1 лінійних алгебраїчних рівнянь, яку запишемо в матричній формі:
, (4.1.31)
де а = (4/(3n)2 + n), в = (1/(3n) – n/2).
З урахуванням граничних умов (4.1.8) запишемо одиниці на головній діагоналі в першому і (n+1)-й рядках, нулі в інших позиціях цих рядків і замінимо значення першого і (n+1)-го елементів в правому стовпці значеннями у в першій і (n+1)-й вузлових точках. Отже, маємо таку систему алгебраїчних рівнянь:
. (4.1.32)
Розв’язуючи (4.1.32), знаходимо значення 1 - n+1.
Як конкретний приклад розглянемо розв’язок задачі (4.1.7)-(4.1.8) з врахуванням чотирьох елементів однакової довжини. Для цього вирази і на підставі формули (4.1.23), (4.1.26) і (4.1.27) запишемо у вигляді:
.
, (4.1.33)
, (4.1.34)
, (4.1.35)
. (4.1.36)
(4.1.37)
Матричну систему рівнянь (4.1.31) п’ятого порядку отримуємо за рахунок асамблювання матриць (4.1.37) та з врахуванням (4.1.25).
Асамблювання матриць для першого та другого елементів таке:
Складемо отримані матриці. В результаті маємо систему рівнянь:
.
Після елементарних перетворень системи отримаємо:
.
Використовуючи перше правило для вузлових точок з умовами Діріхле (якщо р - вузол, в якому вузлове значення задається точно, тобто , то процедура введення умов Діріхле в матричне рівняння системи полягає в наступному: в р-ий рядок матриці системи вносяться нулі, крім діагональної позиції, на яку поміщається 1, в р-ий рядок стовпця вноситься ур), маємо:
Використавши метод Гаусса в програмному забезпеченні, маємо:
=1, = 1,643897, = 2,716637,
= 4,498064, = 7,4529,
де , , , , .
Аналітичний розв’язок:
у1 = 1, у2 = е0,5 , у3 = е , у4 = е1,5 , у5 = е2.
4.1.5 Двовимірний приклад варіаційного
методу скінченних елементів
Як приклад розглянемо рівняння Лапласа:
. (4.1.38)
З граничними умовами Діріхле на одній ділянці границі:
Т=50 , у=0, (4.1.39)
Т=100 , у=L. (4.1.40)
З умовами Неймана на іншій ділянці границі:
(4.1.41)
(4.1.42)
Задача визначена на множині R, яка складається з області D і межі S, тобто R=D+S.
За допомогою варіаційного числення розв’язок Т(х,у), що задовольняє рівняння (4.1.38)-(4.1.42), збігається з функцією, яка мінімізує функціонал
(4.1.43)
де – функція з допустимої множини пробних функцій, заданих в D.
Граничні умови Неймана ( 4.1.41-4.1.42) виконуються автоматично для функції, яка мінімізує функціонал, і, як наслідок, варіаційні формулювання будуть називатися натуральними граничними умовами.
Розіб’ємо область на і скінченних елементів (трикутників). Загальне число вузлів позначимо через n.
В двовимірному і тривимірному випадках немає чіткого зв’язку між загальним числом елементів і загальним числом вузлів. Розбиття області і умови неперервності, що накладаються на пробні функції, дозволяють записати функціонал у вигляді:
, (4.1.44)
де – елементний вклад, що визначається рівністю:
. (4.1.45)
Для довільного елемента еі пробну функцію виберемо у вигляді (4.1.2).
Враховуючи, що в даному випадку Фі=, Фj=, Фk=, згідно з формулами (4.1.5-4.1.6) для кожного елемента отримуємо:
.
(4.1.46)
Тоді маємо:
, (4.1.47)
. (4.1.48)
Підставимо (4.1.47)-(4.1.48) в (4.1.45):
(4.1.49)
Оскільки підінтегральний вираз в (4.1.49) не залежить від х і у, і, крім того, , то рівність (4.1.49) перепишемо в такому вигляді:
. (4.1.50)
Вираз вигляду (4.1.50) може бути отриманий для кожного елемента. Підставляючи (4.1.50) в (4.1.44), перетворимо функціонал (4.1.42) в функцію всіх вузлових значень Т1 –Тn , тобто
(4.1.51)
Умови мінімуму Х можна записати у вигляді:
. (4.1.52)
Підставляючи (4.1.44) в (4.1.52), маємо:
(4.1.53)
Диференціювання (4.1.50) за приводить до виразу:
, (4.1.54)
де значення р, q, і r – відповідно, номери вузлів i, j, k елемента еі в множині номерів вузлів системи.
Об’єднання компонентів елементарних рівнянь, що задаються рівністю (4.1.54), називається об’єднанням за вузлами, тому що процес об’єднання повинен бути виконаний окремо для кожного вузла системи. Формування вкладів та їх об’єднання проілюструємо нижче.
Спочатку область D розіб’ємо на 16 елементів з загальним числом вузлів 15 (рис. 4.4) та запишемо номери вузлів i, j, k для кожного елемента системи в таблицю 4.1. Обхід кожного елемента будемо робити проти годинникової стрілки. Розміри всіх елементів однакові. Координати вузлів запишемо в таблицю 4.2.
Рисунок 4.4
Таблиця 4.1 – Співвідношення між глобальними і локальними номерами вузлів
|
Елемент |
Глобальний номер вузла |
Елемент |
Глобальний номер вузла |
||||
|
I |
J |
k |
I |
J |
K |
||
|
1 |
4 |
1 |
5 |
9 |
10 |
7 |
11 |
|
2 |
2 |
5 |
1 |
10 |
8 |
11 |
7 |
|
3 |
5 |
2 |
6 |
11 |
11 |
8 |
12 |
|
4 |
3 |
6 |
2 |
12 |
9 |
12 |
8 |
|
5 |
7 |
4 |
8 |
13 |
13 |
10 |
14 |
|
6 |
5 |
8 |
4 |
14 |
11 |
14 |
10 |
|
7 |
8 |
5 |
9 |
15 |
14 |
11 |
15 |
|
8 |
6 |
9 |
5 |
16 |
12 |
15 |
11 |
В таблиці 4.3 згідно з формулою (4.1.6) наведені вi, вj, вk, сi, сj і сk. Для ілюстрації розглянемо елемент 9. З таблиці 4.1 видно, що його вузли i, j, k мають номери 10, 7 і 11, відповідно.
Таблиця 4.2 – Координати вузлів
|
Вузол |
Координати |
Вузол |
Координати |
Вузол |
Координати |
|||
|
х |
У |
х |
У |
Х |
У |
|||
|
1 |
0 |
0 |
6 |
2 |
0,5 |
11 |
1 |
1,5 |
|
2 |
1 |
0 |
7 |
0 |
1 |
12 |
2 |
1,5 |
|
3 |
2 |
0 |
8 |
1 |
1 |
13 |
0 |
2 |
|
4 |
0 |
0,5 |
9 |
2 |
1 |
14 |
1 |
2 |
|
5 |
1 |
0,5 |
10 |
0 |
1,5 |
15 |
2 |
2 |
Таблиця 4.3 – Характерні параметри елементів
|
Елементи |
Параметри |
|||||
|
вi |
вj |
вk |
сi |
сj |
сk |
|
|
1 |
-0,5 |
0 |
0,5 |
1 |
-1 |
0 |
|
2 |
0,5 |
0 |
-0,5 |
-1 |
1 |
0 |
|
3 |
-0,5 |
0 |
0,5 |
1 |
-1 |
0 |
|
4 |
0,5 |
0 |
-0,5 |
-1 |
1 |
0 |
|
5 |
-0,5 |
0 |
0,5 |
1 |
-1 |
0 |
|
6 |
0,5 |
0 |
-0,5 |
-1 |
1 |
0 |
|
7 |
-0,5 |
0 |
0,5 |
1 |
-1 |
0 |
|
8 |
0,5 |
0 |
-0,5 |
-1 |
1 |
0 |
|
9 |
-0,5 |
0 |
0,5 |
1 |
-1 |
0 |
|
10 |
0,5 |
0 |
-0,5 |
-1 |
1 |
0 |
|
11 |
-0,5 |
0 |
0,5 |
1 |
-1 |
0 |
|
12 |
0,5 |
0 |
-0,5 |
-1 |
1 |
0 |
|
13 |
-0,5 |
0 |
0,5 |
1 |
-1 |
0 |
|
14 |
0,5 |
0 |
-0,5 |
-1 |
1 |
0 |
|
15 |
-0,5 |
0 |
0,5 |
1 |
-1 |
0 |
|
16 |
0,5 |
0 |
-0,5 |
-1 |
1 |
0 |
Підставляючи ці дані і відповідні значення параметрів для елемента 9 з таблиці 4.3 в рівність (4.1.50), отримуємо вираз:
. (4.1.55)
Верхній індекс 9 у площі Δ вказує на номер елемента, для якого вона визначена. Оскільки площі всіх елементів одинакові , то в подальшому верхній не будемо зображувати.
Згідно з (4.1.55) визначимо похідні :
, (4.1.56)
, (4.1.57)
. (4.1.58)
Рівняння (4.1.56-4.1.58) можна записати у вигляді елементарного матричного рівняння:
, (4.1.59)
або , (4.1.60)
де (4.1.61)
елементарна матриця жорсткості елемента е і Те – елементарний вузловий вектор. Для спрощення записів в подальшому верхній індекс еі в рівняннях (4.1.60) і (4.1.61) будемо опускати.
Оскільки елементи 1, 3, 5, 7, 9, 11, 13, 15 мають однакові розміри і орієнтацію відносно системи Оху , то елементарна матриця жорсткості k однакова з врахуванням (4.1.60) у відповідних номерах вузлів. Так, для елемента 5 матричне рівняння має вигляд:
(4.1.62)
Інші елементи 2, 4, 6, 8, 10, 12, 14, 16 також мають загальну елементну матрицю жорсткості k, яка в даному випадку тотожна матриці жорсткості з непарними номерами. Так, для елемента 10, матричне рівняння має вигляд:
. (4.1.63)
Після визначення похідних для всіх елементів покажемо їх об’єднання по вузлах згідно з формулою (4.1.53) для отримання матричного рівняння системи. Наприклад, для вузла 7 дають внесок лише елементи 5, 9 і 10, отже, маємо:
. (4.1.64)
Підставляючи в (4.1.64) внески елементів 5, 9 і 10 з врахуванням (4.1.59), (4.1.62) і (4.1.63), отримаємо:
(4.1.65)
Рівняння (4.1.65) запишемо в розширеній формі:
,
(4.1.66)
де вузловий вектор Т системи має вигляд:
. (4.1.67)
Для кожного вузла системи рівняння аналогічні (4.1.65).
В результаті матричне рівняння системи має вигляд:
. (4.1.68)
Для корегування отриманої системи необхідно враховувати граничні умови Діріхле. Оскільки у вузлах 1, 2, 3 і 13, 14, 15 значення задані явно ((4.1.39),(4.1.40)), то, згідно з методом Пейна-Айронса, в рядках матриці k 1, 2, 3, 13, 14 і 15 внесемо нулі, крім діагональної позиції, на яку запишемо одиницю, а в відповідні рядки стовпця в правій частині вводимо значення вузлів. В результаті отримаємо таку матричну систему рівнянь:
(4.1.69)
Використовуючи чисельні методи (програми) розрахунку системи алгебраїчних рівнянь, визначаємо розв’язок (4.1.69):
Т1 = 50, Т2 = 50, Т3 =50
Т4 =62,5, Т5 = 62,5, Т6 = 62,5
Т7 = 75, Т8 = 75, Т9 = 75
Т10 = 87,5, Т11 = 87,5, Т12 = 87,5
Т13 = 100, Т14 = 100, Т15 = 100.
Отже, наприклад, для досліджуваної області в точці А з координатами (l/2,l/4) температура Т/А≈ Т11=87,5˚С.
Одним з методів розв’язання задач оптимізації є метод програмування, який полягає в побудові, аналізі і чисельному розв’язанні лінійних оптимізаційних моделей. Задача лінійного програмування має місце лише тоді, коли обмеження і цільова функція об’єкта управління лінійні відносно факторів Х. Метод використовується для розв’язання виробничих і управлінських задач в різноманітних галузях народного господарства, якщо параметри об’єкта моделювання відповідають вимогам аксіом пропорційності (ділення) і адитивності.
Принцип ділення (пропорційності), наприклад, у виробничому процесі означає, що загальний (сумарний) об’єм застосування кожного виду ресурсів і, відповідно, отриманий при цьому прибуток підприємства прямо пропорційні об’єму виробленої продукції. Так, якщо збільшити вдвічі об’єм продукції, то вживання усіх ресурсів при цьому також виростуть вдвічі.
Адитивність в цих випадках полягає в тому, що прибуток підприємства від реалізації продукції є сума вкладів від реалізації різноманітних її видів. Одночасно повна кількість кожного з використовуваних ресурсів дорівнює сумі однойменних ресурсів, які затрачуються на кожному виробництві всіх видів продукції (наприклад, загальне використання на заводі цементу за день складається з його витрат на виробництво всієї денної номенклатури виробів).
Приклад
Завод ЗБВ випускає декоративні плити двох типів – П1 і П2, витрачаючи кольорові мінеральні наповнювачі (КМН) і кольорові мінеральні заповнювачі (КМЗ), одночасні ресурси яких обмежені; прибуток від реалізації продукції П1 і П2 різний (табл.5.1).
Показник |
Тип плит |
Ресурси матеріалу, кг |
|
|
П1 |
П2 |
||
|
Витрата кольорових матеріалів, кг наповнювач КМН заповнювач КМЗ |
70 75 |
40 100 |
168000 240000 |
|
Прибуток від випуску однієї плити, грн. |
3 |
2 |
— |
Змінні Х, тобто величини, які належить визначити в результаті розв’я-зання задачі ЛП, в даному випадку – об’єми виробництва плит Х1, шт., типу П1 і Х2, шт., типу П2.
Прибуток від реалізації всіх плит П1 буде 3Х1 грн. (виконується аксіома пропорційності), а від всіх плит П2 складе 2Х2 грн., тобто загальний прибуток у відповідності з принципом адитивності:
. (5.1)
Ця функція називається цільовою. Тобто, необхідно знайти такі значення Х1 і Х2, щоб прибуток С був максимальним. При цьому не можна перевищувати ресурси ні по одному з матеріалів:
, (5.2)
. (5.3)
До отриманих рівнянь необхідно додати умову невід’ємності змінних як натуральних результатів виробництва: Х1≥0, Х2≥0.
Розв’яжемо задачу геометрично в координатах Х1 і Х2. Область допустимих значень Ω{Х1,Х2} визначається на першому етапі розв’язання задачі. Внаслідок невід’ємності змінних вона знаходиться в першому квадранті і обмежена координатними осями (Х1=Х2=0) і двома прямими:
75Х1 + 100Х2=240000 (проходить через точки М{3200;0} і М´{0;240000}) і 70Х1 + 40Х2=168000 (проходить через точки К{2400;0} і К´{0;1680}), які створюють опуклий багатокутник. В кожній точці всередині нього (незаштрихована область Ω) виконуються всі обмеження задачі.
Х2, шт.
k´
4000
3000 7500
M
2000
1000 4500 L
. 6000
Ω
Х1, шт.
1000 2000 k 3000 M´
Рисунок 5.1
Наступний етап в розв’язанні задачі оснований на аналізі цільової функції С – на діаграму наносяться її ізолінії при послідовно зростаючих на ΔС=const значеннях загального прибутку. З рис.5.1 видно, що максимально можливе значення Сmax в допустимій області Ω{Х1,Х2} досягається у вершині L багатокутника – на перетині прямих. Координати її визначаються з рівнянь (5.2) і (5.3), прирівнюючи їх відносно однієї із змінних (наприклад, Х1):
Звідки Х2{L}=1050 і Х1{L}=1800. Тоді значення цільової функції в точці максимуму С{L}=3·1800+2·1050=7500 грн.
Отримані результати можна перевірити за рахунок аналітичного розв’язання задачі.
Запишемо задачу в канонічній формі. Для цього введемо в ліву частину (5.2), (5.3) додаткові невід’ємні змінні і :
(5.4)
(5.5)
Цільова функція для системи (5.4)-(5.5) аналогічна (5.1):
(5.6)
Отже, для визначення невідомих маємо меншу кількість рівнянь, тобто система є невизначеною.
Одне з можливих рішень невизначеної системи можна знайти, прирівнюючи які-небудь дві змінних до нуля. Наприклад, Х1=Х2=0, тоді Х3=168000 і Х4=240000, причому цільова функція (5.6) в цій точці С=0.
Розглянемо можливі варіанти розв’язання системи (5.4), (5.5):
.
Отже, в даному випадку, отримуємо недопустимий розв’язок (Х4<0).
Значення цільової функції С=4800.
Знову отримали недопустимий розв’язок (Х3<0).
Отже, отримали шість розв’язків задачі, з яких один оптимальний. Кількість розв’язків в загальному випадку можна отримати за формулою:
, (5.7)
де n – кількість рівнянь;
m – кількість невідомих.
Ми впевнились, що навіть при n=2 число варіантів розв’язання дорівнює 6. Тобто, з точки зору геометричного розв’язання задачі прямий перебір всіх вершин багатокутника є занадто громіздким.
Систематичний обчислювальний метод, який дозволяє знаходити чисельні рішення для лінійних оптимізаційних моделей за менше число кроків, розроблений Дж. Данцігом і називається симплекс-методом [11]. Він дозволяє аналізувати і правильно інтерпретувати отримані розв’язки на кожному кроці виконання алгоритму.
Процес розв’язання починається з визначення початкової вершини багатогранника розв’язків, координати якої створюють початковий допустимий базисний розв’язок, або початковий опорний план. Оптимальний розв’язок знаходиться як результат послідовного покращення проміжного базисного розв’язку в напрямку градієнта цільової функції. Перехід від одного допустимого розв’язку до другого відбувається по тому ребрі багатокутника розв’язків, яке має найбільшу за абсолютним значенням проекцію градієнта С.
У виразі цільової функції С ненулеві коефіцієнти відповідають небазисним змінним. Для отримання оптимального розв’язку С виражають через інші (базисні) змінні, які, таким чином, здобувають від’ємні коефіцієнти. Оскільки всі змінні невід’ємні, то С досягає максимуму, коли ці базисні змінні дорівнюють нулю. Ця процедура – сутність операції заміни базису.
Координати кожної наступної вершини багатокутника визначаються заміною однієї з нульових змінних проміжного розв’язання проміжною базисною змінною. Оскільки кількість вершин багатокутника кінцева, то оптимальні розв’язки знаходяться за кінцеву кількість кроків, причому меншу, ніж кількість вершин.
Ці положення узагальнюються ітеративним алгоритмом, який називається симплекс-метод.
Симплексний метод задачі лінійного програмування оснований на переході від одного опорного плану до іншого, при якому значення цільової функції збільшується. Вказаний перехід можливий, якщо відомий який-небудь початковий опорний план. Розглянемо задачу, для якої цей план можна записати.
Нехай необхідно знайти максимальне значення функції
за умов:
Тут , і () – задані постійні числа (mn і ві0).
Векторна форма даної задачі має такий вигляд: знайти максимум функції
(5.1.1)
за умов:
(5.1.2)
(5.1.3)
де
Оскільки
(5.1.4)
то за означенням опорного плану x=(в1;в2;…;0;…;0) є опорним планом даної задачі. План x=(в1;в2;…;вn) називається опорним планом основної задачі лінійного програмування, якщо система векторів Pj в (5.1.4) з додатними коефіцієнтами вj, лінійно незалежна. Цей план визначається системою одиничних векторів Р1, Р2, ..., Рm, які створюють базис m-вимірного простору. Тому кожний з векторів Р1, Р2, ..., Рn а також вектор Р0 можуть бути подані у вигляді лінійної комбінації векторів даного базису. Нехай
Покладемо: .
Оскільки вектори Р1, Р2, ..., Рm – одиничні, то:
і , а .
Ознаку оптимальності опорного плану визначають за теоремами:
Теорема 1.
Опорний план Х*= задачі (5.1.1) – (5.1.3) є оптимальним, якщо для будь-якого .
Теорема 2.
Якщо для деякого і серед чисел немає додатних , то цільова функція (5.1.1) задачі (5.1.1) – (5.1.3) не обмежена на множині її планів.
Теорема 3.
Якщо опорний план Х задачі (5.1.1) – (5.1.3) невироджений і , але серед чисел є додатні (не всі ), то існує опорний план такий, що
Дослідження опорного плану на оптимальність, а також подальший обчислювальний процес зручно вести, якщо умову задачі і початкові дані отримані після визначення відповідного опорного плану в таблиці 5.2.
Таблиця 5.2
|
і |
Базис |
Сб |
Р0 |
С1 |
C2 |
... |
Сr |
… |
Cm |
Cm+1 |
… |
Ck |
… |
Cn |
|
P1 |
P2 |
… |
Pr |
... |
Pm |
Pm+1 |
… |
Pk |
… |
Pn |
||||
|
m+1 |
F0 |
0 |
0 |
0 |
В стовпці Сб цієї таблиці записують коефіцієнти при невідомих цільової функції, які мають такі ж індекси, що і вектори даного базису. В стовпці Р0 записують додатні компоненти опорного плану, в ньому ж в результаті розрахунків отримують додатні компоненти опорного плану. Стовпці векторів являють собою коефіцієнти розкладання цих векторів за векторами даного базису.
В таблиці перші m рядків визначаються даними задачі, а показники (m+1)-ого рядка визначають. В цьому рядку в стовпці вектора P0 записують значення цільової функції, які вона приймає при даному опорному плані, а в стовпці Pj – значення
Значення zj знаходяться як скалярний добуток вектора Pj () на вектор Сб=(с1;с2;...;сm):
.
Значення F0 дорівнює скалярному добутку вектора Р0 на вектор Сб:
Після заповнення таблиці опорний план перевіряють на оптимальність. Для цього переглядають елементи (m+1)-ого рядка таблиці. В результаті можуть мати місце такі випадки:
Перший випадок є оптимальним. В другому випадку цільова функція не обмежена зверху на множині планів, а в третьому випадку можна перейти до нового опорного плану, при якому значення цільової функції збільшиться. Цей перехід від одного опорного плану до другого відбувається вилученням з базису якого-небудь з векторів, введенням в нього нового вектора. Як новий можна взяти будь-який з векторів , для якого . Нехай, наприклад, <0, і введемо в базис вектор .
Для визначення вектора, який належить вилученню з базису, знаходять для всіх >0. Нехай цей мінімум досягається при . Тоді з базису вилучають вектор Pr, а число називають елементом розв’язку.
Стовпці, на перетині яких знаходиться елемент розв’язку, називають напрямними.
Після виділення напрямного рядка і напрямного стовпця знаходять новий опорний план, коефіцієнти розкладу векторів через вектори нового базису, відповідного новому опорному плану. Це реалізується методом Жордана-Гаусса. При цьому додатні компоненти нового опорного плану обчислюються за формулами:
(5.1.5)
а коефіцієнти розкладу векторів – через вектори нового базису, відповідного новому опорному плану, за формулами:
(5.1.6)
Метод Жордана-Гаусса полягає в тому, що після приведення відповідної матриці до трикутного вигляду продовжують елементарні перетворення рядків останньої до діагональної форми. В результаті за методом Жордана-Гаусса розрахунок відбувається аналогічно способу Гаусса, але розв’язок системи отримується одразу без зворотного ходу (в лівій частині буде одинична матриця).
Після знаходження і згідно з формулами (5.1.5)-(5.1.6) їх значення заносять в таблицю 5.3. Елементи (m+1) - го рядка цієї таблиці можуть бути визначені за формулами:
Наявність двох стовпців знаходження елементів (m+1)-ого рядка дозволяє провести контроль правильності проведення розрахунків.
Таблиця 5.3
|
і |
Базис |
Сб |
Р0 |
С1 |
C2 |
... |
Сr |
… |
Cm |
Cm+1 |
… |
Ck |
… |
Cn |
|
P1 |
P2 |
… |
Pr |
... |
Pm |
Pm+1 |
… |
Pk |
… |
Pn |
||||
|
m+1 |
F0 |
0 |
0 |
0 |
З формули (5.1.7) маємо, що при переході від одного опорного плану до іншого найбільш раціонально ввести в базис вектор , при якому максимальним за абсолютною величиною є число <0>0). В подальшому вектор, який вводиться в базис, визначати, виходячи з максимальної абсолютної величини від’ємних чисел . Якщо таких чисел декілька, то в базис будемо вводити вектор, який має такий же індекс, як і максимальне з чисел Cj, що визначаються даними числами (<0).
Перехід від одного опорного плану до другого зводиться до переходу від однієї симплекс-таблиці до іншої. Елементи її вираховуються за допомогою рекурентних формул (5.1.5)-(5.1.8) за такими правилами:
Правило 1.
В стовпцях векторів, які входять в базис, на перетині рядків і стовпців однойменних векторів проставляються одиниці, а всі інші елементи даних стовпців вважають рівними нулеві.
Правило 2.
Елементи векторів Р0 і Рj в рядку нової таблиці, в якому записаний вектор, введений в базис, отримують з елементів цього ж рядка попередньої таблиці діленням їх на величину вирішуваного елемента. В стовпці Сб в рядку вектора, що вводиться, проставляють величину Сk, де k – індекс вирішуваного вектора.
Правило 3.
Інші елементи стовпців вектора P0 і Pj нової таблиці обчислюють за правилом трикутника. Для знаходження будь-якого з цих елементів знаходять три числа:
Ці три числа створюють трикутник, дві вершини якого відповідають числам, що знаходяться в даній таблиці, а третя – числу, що знаходиться в новій симплекс-таблиці. Для визначення шуканого елемента нової таблиці від першого числа віднімають добуток другого і третього.
Після заповнення нової симплекс-таблиці переглядають елементи (m+1)- го рядка. Якщо всі , то новий опорний план є оптимальним. Якщо серед вказаних чисел є від’ємні, то, використовуючи описану вище послідовність дій, знаходять новий опорний план. Цей процес продовжують до тих пір, поки не отримають оптимальний план задачі або встановлять неможливість її розв’язання.
Як приклад розглянемо розв’язання ще однієї задачі.
Для виготовлення різноманітних виробів А, В і С підприємство використовує три види сировини. Норма витрат сировини на виготовлення одного виробу кожного виду, ціна одного виробу А, В і С, а також загальна кількість сировини кожного виду, що може бути використана підприємством, наведені в таблиці 5.4.
Таблиця 5.4
|
Вид сировини |
Норма витрат сировини (кг) на один виріб |
Загальна кількість сировини (кг) |
||
|
А |
В |
С |
||
|
І ІІ ІІІ |
18 6 5 |
15 4 3 |
12 8 3 |
360 192 180 |
|
Ціна одного виробу (грн.) |
9 |
10 |
16 |
Скласти план виготовлення виробів, при якому загальна вартість всієї виробленої підприємством продукції є максимальною.
Розв’язання.
Складемо математичну модель задачі. Позначимо кількість виробів А через х1, виробів В – через х2, виробів С – через х3. Згідно з умовою задачі маємо таку систему нерівностей:
(5.1.9)
Загальна вартість виготовленої підприємством продукції запишеться у вигляді:
. (5.1.10)
За своїм економічним змістом змінні х1, х2 і х3 можуть приймати тільки невід’ємні значення:
х1,х2,х3≥0. (5.1.11)
Отже, необхідно визначити серед розв’язків (5.1.9) такі, при яких функція (5.1.10) приймає максимальне значення.
Введемо три додаткові змінні: х4, х5, х6 – невикористана кількість сировини, відповідно. Тоді систему нерівностей (5.1.9) запишемо у вигляді системи рівнянь:
(5.1.12)
Систему (5.1.12) запишемо у векторній формі:
де
Оскільки серед векторів Р1–Р6 є три одиничних вектори, то для даної задачі можна записати опорний план, який вважається системою тривимірних одиничних векторів Р4, Р5, Р6, що створюють базис тривимірного векторного простору, тобто
х=(0;0;0;360;192;180).
Складемо симплексну таблицю для першої ітерації, підрахуємо значення і перевіримо початковий опорний план на оптимальність:
Для векторів базису
Таблиця 5.5
|
і |
Базис |
Сб |
Р0 |
С1=9 |
С2=10 |
С3=16 |
С4=0 |
С5=0 |
С6=0 |
|
Р1 |
Р2 |
Р3 |
Р4 |
Р5 |
Р6 |
||||
|
1 2 3 |
Р4 Р5 Р6 |
0 0 0 |
360 192 180 |
18 6 5 |
15 4 3 |
12 8 3 |
1 0 0 |
0 1 0 |
|
|
F0=0 |
-9 |
-10 |
-16 |
0 |
0 |
0 |
З таблиці видно, що значення цільової функції дорівнює нулеві, тобто ніщо не виробляється, сировина не використовується. Отже, план не є оптимальним.
В 4-ому рядку таблиці 5.5 від’ємні числа свідчать про можливість збільшення загальної вартості виробленої продукції і показують, на скільки збільшиться ця сума при введенні в план одиниці того чи іншого виду продукції. Так, число 9 означає, що при включенні в план виробництва одного виробу А забезпечується збільшення вартості випуску продукції на 9 грн. Якщо включити в план виробництва по одному виробу В і С, то загальна вартість виготовленої продукції виросте, відповідно, на 10 і 16 грн. Тому з економічної точки зору найбільш раціональним є включення в план виробництва вибір С. Оскільки число 16 стоїть в 4-ому рядку стовпця Р3, то в базис вводимо вектор Р3. Визначимо вектор, який належить виключенню з базису. Для цього знаходимо:
Отже, вектор Р5 належить вилученню з базису. Стовпець вектора Р3 і другий рядок є напрямними. Складемо таблицю для 2-ї ітерації (табл.5.6).
Таблиця 5.6
|
і |
Базис |
Сб |
Р0 |
9 |
10 |
16 |
0 |
0 |
0 |
|
Р1 |
Р2 |
Р3 |
Р4 |
Р5 |
Р6 |
||||
|
9 |
0 0 1 |
||||||||
|
4 |
384 |
3 |
-2 |
0 |
0 |
2 |
0 |
Спочатку заповнюємо рядок вектора Р3. Елементи цього рядка отримуються з відповідних елементів рядка 2 таблиці 5.3 діленням на вирішуваний елемент 8. При цьому в стовпці Сб записуємо коефіцієнт С3=16. Потім заповнюємо елементи стовпців для векторів, які входять в повний базис. В цих стовпцях на перетині рядків і стовпців однойменних векторів проставляємо одиниці, а всі інші елементи вважають рівними нулеві.
Для визначення інших елементів таблиці 5.6 використовуємо правило трикутника.
Визначаємо елементи таблиці 5.6, які стоять в стовпці вектора Р0.
Для першого елемента маємо:
Для третього елемента - .
Значення F0 в 4-ому рядку стовпця цього ж вектора визначимо за формулою:
.
За правилом трикутника визначимо елементи стовпця вектора Р1. Перші два числа з таблиці 5.5, а третє число – з таблиці 5.6. Це число стоїть на перетині 2-ого рядка і стовпця вектора Р1 останньої таблиці. В результаті маємо значення невідомих елементів:
Число в 4-ому рядку стовпця вектора Р1 таблиці 5.6 визначимо за формулою:
Аналогічно знаходимо елементи стовпця вектора Р2. Елементи стовпця вектора Р5 визначаємо за правилом трикутника.
Для першого елемента маємо:
Для третього елемента маємо:
В результаті з таблиці 5.6 отримуємо новий опорний план задачі При даному плані виробництва виготовляються 24 вироби С і залишаються невикористаними 72 кг сировини 1-го виду і 108 кг сировини 3-го виду. Вартість всієї виробленої продукції при цьому плані дорівнює 384 грн. Але план не є оптимальним. Це видно з 4-ого рядка таблиці 5.6, оскільки у стовпці вектора Р2 цього рядка стоїть від’ємне число 2. Отже, в базис необхідно ввести вектор Р2, тобто в новому плані необхідно передбачити випуск виробів В. Можливий випуск виробів В визначається для >0, тобто знаходимо:
Отже, виключенню з базису підпадає вектор Р4. Число 9 є вирішуваним елементом, а стовпець вектора Р2 і перший рядок таблиці 5.6 є напрямними. Складаємо таблицю для третьої ітерації (таблиця 5.7).
В таблиці 5.7 спочатку заповнюємо елементи першого рядка, який являє собою рядок з вектором Р2 в базисі. Елементи цього рядка отримуємо з елементів першого рядка таблиці 5.6 діленням останніх на вирішуваний елемент (тобто на 9). При цьому в стовпці Сб даного рядка записуємо С2=10.
Потім заповнюємо елементи стовпців векторів базису і за правилом трикутника вираховуємо елементи останніх стовпців.
Так, для стовпця вектора Р0 маємо:
другий елемент:
третій елемент:
Значення F0 в 4-ому рядку стовпця цього ж вектора визначимо за формулою:
Таблиця 5.7
|
і |
Базис |
Сб |
Р0 |
9 |
10 |
16 |
0 |
0 |
0 |
|
Р1 |
Р2 |
Р3 |
Р4 |
Р5 |
Р6 |
||||
|
1 |
0 0 1 |
||||||||
|
4 |
400 |
5 |
0 |
0 |
0 |
Для стовпця вектора Р1 визначаємо елементи 2-4 рядків.
Другий елемент:
Третій елемент:
Аналогічно вираховуємо елементи для векторів Р4 і Р5 (елементи другого і третього рядків).
Наприклад, для вектора Р5:
елемент другого рядка:
;
елемент третього рядка:
Розглянемо четвертий рядок таблиці 5.7. В цьому рядку серед чисел немає від’ємних. Це означає, що знайдений опорний план є оптимальним і
Отже, план випуску продукції, включаючи виготовлення 8 виробів В і 20 виробів С, є оптимальним. Сировина першого і другого видів використовується повністю. При цьому залишається 96 кг сировини третього виду, а вартість виробленої продукції рівна 400грн.
5.2 Транспортна задача
Математична постановка задачі.
Загальна постановка транспортної задачі полягає у визначенні оптимального плану перевезення деякого однорідного вантажу з т пунктів відправлення А1, А2,..., Ат у п пунктів призначення В1, В2,..., Вп. При цьому як критерій оптимальності звичайно береться або мінімальна вартість перевезень усього вантажу, або мінімальний час його доставки. Розглянемо транспортну задачу, за критерій оптимальності якої взята мінімальна вартість перевезень усього вантажу. Позначимо через тарифи перевезення одиниці вантажу з і-го пункту відправлення в j-й пункт призначення, через аj – запаси вантажу в j-му пункті відправлення, через bj – потреби у вантажі в j-ому пункті призначення, а через хі j – кількість одиниць вантажу, перевезеного з і-ого пункту відправлення в j-ий пункт призначення. Тоді математична постановка задачі полягає у визначенні мінімального значення функції
(5.2.1)
за умов
; (5.2.2)
; (5.2.3)
. (5.2.4)
Оскільки змінні задовольняють системи лінійних рівнянь (5.2.2) і (5.2.3) і умову незаперечності (5.2.4), забезпечується доставка необхідної кількості вантажу в кожний з пунктів призначення, вивезення наявного вантажу з усіх пунктів відправлення, а також включаються зворотні перевезення.
Означення 1. Будь-який негативний розв’язок систем лінійних рівнянь (5.2.2) і (5.2.3), обумовлений матрицею , називається планом транспортної задачі.
Означення 2. План , при якому функція (5.2.1) приймає своє мінімальне значення, називається оптимальним планом транспортної задачі.
Звичайно вихідні дані транспортної задачі записують у вигляді табл.5.8.
Очевидно, загальна наявність вантажу в постачальників дорівнює, а загальна потреба у вантажі в пунктах призначення дорівнює одиниць. Якщо загальна потреба у вантажі в пунктах призначення дорівнює
запасу вантажу в пунктах відправлення, тобто
, (5.2.5)
то модель такої транспортної задачі називається закритою. Якщо ж зазначена умова не виконується, то модель транспортної задачі називається відкритою.
Таблиця 5.8.
|
Пункти відправлення |
Пункти призначення |
Запаси |
||||
|
В1 |
.... |
Вj |
… |
Bn |
||
|
А1 |
c11 X11 |
... |
c1j X1j |
... |
c1n X1n |
а1 |
|
... |
… |
... |
... |
... |
... |
... |
|
Аі |
ci1 Xi1 |
... |
cij Xij |
... |
cin Xin |
аі |
|
... |
... |
... |
... |
... |
... |
|
|
Ат |
cm1 Xm1 |
... |
cmj Xmj |
... |
cmn Xmn |
ат |
|
Потреби |
b1 |
... |
bj |
… |
bn |
Теорема 1. Для розв’язання транспортної задачі необхідно і достатньо, щоб запаси вантажу в пунктах відправлення були рівні потребам у вантажі в пунктах призначення, тобто, щоб виконувалась рівність (5.2.5).
У випадку перевищення запасу над потребою, тобто , вводиться фіктивний (п+1)-ий пункт призначення з потребою і відповідні тарифи вважаються рівними нулю: . Отримана задача є транспортною задачею, для якої виконується рівність (5.2.5).
Аналогічно, при вводиться фіктивний (т+1)-ий пункт відправлення з запасом вантажу і тарифи покладаються рівними нулю: . Цим задача зводиться до звичайної транспортної задачі, з оптимального плану якої виходить оптимальний план вихідної задачі. Надалі будемо розглядати закриту модель транспортної задачі. Якщо ж модель конкретної задачі є відкритою, то, виходячи із сказаного вище, перепишемо таблицю умов задачі так, щоб виконувалась рівність (5.2.5).
Число змінних у транспортній задачі з т пунктами відправлення і n пунктами призначення дорівнює nm, а число рівнянь у системах (5.2.2) і (5.2.3) дорівнює п+т. Оскільки ми припускаємо, що виконується умова (5.2.5), то число лінійно незалежних рівнянь дорівнює п+т-1. Отже, опорний план транспортної задачі може мати не більше п+т-1 відмінних від нуля невідомих.
Якщо в опорному плані число відмінних від нуля компонентів дорівнює в точності п+т-1, то план є невиродженим, а якщо менше – то виродженим.
Для визначення опорного плану існує кілька методів. Три з них – метод північно-західного кута, метод мінімального елемента і метод апроксимації Фогеля – розглядаються нижче.
Як і для всякої задачі лінійного програмування, оптимальний план транспортної задачі є й опорним планом.
Для визначення оптимального плану транспортної задачі можна використовувати викладені вище методи. Однак через виняткову практичну важливість цієї задачі і специфіку її обмежень (кожна невідома входить лише в два рівняння систем (5.2.2) і (5.2.3) і коефіцієнти при невідомих дорівнюють одиниці) для визначення оптимального плану транспортної задачі розроблені спеціальні методи. Один з них – метод потенціалів – розглянутий нижче.
Задача 1.
Чотири підприємства даного економічного району для виробництва продукції використовують три види сировини. Потреби в сировині кожного з підприємств відповідно рівні 120, 50, 190 і 110 од. Сировина зосереджена в трьох місцях, а запаси відповідно рівні 160, 140, 170 од. На кожне з підприємств сировина може завозитися з будь-якого пункту її отримання. Тарифи перевезень є відомими величинами і задаються матрицею:
.
Скласти такий план перевезень, при якому загальна вартість перевезень є мінімальною.
Розв’язання.
Позначимо через кількість одиниць сировини, перевезеної з 1-ого пункту її отримання на j-е підприємство. Тоді умови доставки і вивезення необхідної і наявної сировини забезпечуються за рахунок виконання таких рівностей:
(5.2.6)
При даному плані Х=(і=1,3; j=1,4) перевезень загальна її вартість складе:
F=7x11+8x12+x13+2x14+4x21+5x22+9x23+8x24+9x31+2x32+3x33+6x34. (5.2.7)
Таким чином, математична постановка даної транспортної задачі полягає в знаходженні такого невід’ємного розв’язку системи лінійних рівнянь (5.2.6), при якому цільова функція (5.2.7) приймає мінімальне значення.
5.2.1 Метод північно-західного кута
При знаходженні опорного плану транспортної задачі методом північно-західного кута на кожному кроці розглядаємо перший з пунктів відправлення, що залишилися, та перший з пунктів призначення, що залишилися. Заповнення клітинок таблиці умови починається з лівої верхньої клітинки для невідомого x11 ("північно-західний кут") і закінчується клітинкою для невідомого хmn, тобто наче йде по діагоналі таблиці.
На три бази А1, А2, А3 надійшов однорідний вантаж в кількостях, відповідно рівних 140, 180 і 160 одиниць. Цей вантаж треба перевезти в п'ять пунктів призначення В1, В2, В3, В4, В5 відповідно в кількостях 60, 70, 120, 130 і 100 одиниць. Тарифи перевезення одиниць вантажу з кожного з пунктів відправлення у відповідні пункти призначення вказані в таблиці:
Таблиця 5.9
|
Пункти відправлення |
Пункти призначення |
Запас |
||||
|
B1 |
В2 |
В3 |
В4 |
В5 |
||
|
А1 |
2 |
3 |
4 |
2 |
4 |
140 |
|
А2 |
8 |
4 |
1 |
4 |
1 |
180 |
|
Аз |
9 |
7 |
3 |
7 |
2 |
160 |
|
Потреби |
60 |
70 |
120 |
130 |
100 |
480 |
Знайти план перевезень даної транспортної задачі методом північно-західного кута.
Розв'язання.
Тут кількість пунктів відправлення m = 3, а кількість пунктів призначення n = 5. Отже, опорний план задачі визначається числами, які знаходяться в 5+3-1=7 заповнених клітинках.
Заповнення таблиці почнемо з клітинки для невідомого x11, тобто будемо намагатися задовольнити потреби першого пункту призначення за рахунок запасів першого пункту відправлення. Оскільки запаси пункту А1 більші, ніж потреби пункту В1 то вважаємо х11 = 60, записуємо це значення у відповідній клітинці таблиці 5.11 і тимчасово виключаємо з розгляду стовпець В1, вважаючи при цьому запаси пункту А1 рівними 80.
|
Пункти відправлення |
Пункти призначення |
Запас |
||||
|
B1 |
В2 |
В3 |
В4 |
В5 |
||
|
А1 |
2 60 |
3 70 |
4 10 |
2 |
4 |
140 |
|
А2 |
8 |
4 |
1 110 |
4 70 |
1 |
180 |
|
Аз |
9 |
7 |
3 |
7 60 |
2 100 |
160 |
|
Потреби |
60 |
70 |
120 |
130 |
100 |
480 |
Розглянемо перший з пунктів відправлення А1 та пункт призначення В2. Запаси пункту А1 більше потреб пункту В2. Покладемо х12 = 70, запишемо це значення у відповідній клітинці таблиці 5.10 і тимчасово виключимо з розгляду стовпець В2. В пункті А1 запаси вважаємо рівними 10 одиницям. Знову розглянемо перші з пунктів відправлення А1 та призначення В3, що залишилися. Потреби пункту В3 більші запасів, що залишилися в пункті А1. Покладемо x13 = 10 та виключимо з розгляду стовпець А1. Значення х13 = 10 запишемо у відповідну клітинку таблиці 5.10 і будемо вважати потреби пункту В3 рівними 110 одиницям.
Тепер перейдемо до заповнення клітинок для невідомого х23 і т.д. Через шість кроків залишається один пункт відправлення А3 з запасом вантажу 100 одиниць і один пункт призначення В5 з потребою в 100 одиниць. Відповідно є одна вільна клітинка, яку і заповнюємо, вважаючи x35 =100 (таблиця 5.10). В результаті отримуємо опорний план:
.
Згідно з даним планом перевезень, загальна вартість перевезень всього вантажу складає:
F = 2·60 + 3·70 +4·10 +1·110 +4·70 +7·60 +2·100 = 1380.
Суть методу мінімального елемента полягає у виборі клітинки з мінімальним тарифом. Слід відмітити, що цей метод, як правило, дозволяє знайти опорний план транспортної задачі, при якому загальна вартість перевезень вантажу менша, ніж загальна вартість перевезень при план, знайденому для даної задачі за допомогою методу північно-західного кута. Тому найбільш доцільно опорний план транспортної задачі знаходити методом мінімального елемента.
Приклад.
Вхідні дані задачі запишемо у вигляді таблиці 5.11.
Таблиця 5.11
|
Пункти відправлення |
Пункти призначення |
Запас |
||||
|
B1 |
В2 |
В3 |
В4 |
В5 |
||
|
А1 |
2 |
3 |
4 |
2 |
4 |
140 |
|
А2 |
8 |
4 |
1 |
4 |
1 |
180 |
|
Аз |
9 |
7 |
3 |
7 |
2 |
160 |
|
Потреби |
60 |
70 |
120 |
130 |
100 |
480 |
Покладемо х23 = 120 і запишемо це значення у відповідну клітинку таблиці 5.12. Виключимо тимчасово з розгляду стовпець В3.
В частині таблиці, що залишилася з трьома рядками і з чотирма стовпцями В1, В2, Вз і В5, клітинка з найменшим значенням тарифу С25 знаходиться на перетині рядка А2 і стовпця В5. Покладемо х25 = 60 і внесемо це значення у відповідну клітинку таблиці 5.12.
|
Пункти відправлення |
Пункти призначення |
Запас |
||||
|
B1 |
В2 |
В3 |
В4 |
В5 |
||
|
А1 |
2 10 |
3 |
4 |
2 130 |
4 |
140 |
|
А2 |
8 |
4 |
1 120 |
4 |
1 60 |
180 |
|
Аз |
9 20 |
7 90 |
3 |
7 |
2 40 |
160 |
|
Потреби |
60 |
70 |
120 |
130 |
100 |
480 |
Тимчасово виключимо з розгляду рядок А2. Оскільки тариф 2 зустрічається у клітинках тричі, то доцільно вибрати тариф С14. Покладемо х14 = 130, стовпець В4 виключаємо з розгляду. Після цього з врахуванням мінімальних тарифів заповнюємо стовпець В5: х35 = 40 і стовпець В1, x11 = 10.
Знову розглядаємо частину таблиці, що залишилася. В ній мінімальний тариф Сij знаходиться в клітинці на перетині рядка А3 і стовпця В2, що дорівнює 7.
Заповнимо наведеним вище способом цю клітинку і аналогічно заповнюємо (у відповідній послідовності) клітинку, що знаходиться на перетині рядка А3 і стовпця В1. В результаті отримаємо опорний план:
.
При даному плані перевезень загальна вартість транспортування складає:
F = 2·10 + 9·20 + 7·90 + 1·120 + 2·130 + 1·60 + 2·40 = 1170.
5.2.3 Метод апроксимації Фогеля
При визначенні оптимального плану транспортної задачі методом апроксимації Фогеля на кожній ітерації по всіх стовпцях і по всіх рядках знаходять різницю між двома записаними в них мінімальними тарифами. Ці різниці записують в спеціально відведених для цього рядку та стовпці в таблиці умови задачі. В рядку (чи в стовпці), якому дана різниця відповідає, визначають мінімальний тариф. Клітинку, в якій він записаний, заповнюють на даній ітерації.
Якщо мінімальний тариф однаковий для декількох клітинок даного рядка (стовпця), то для заповнення вибирають ту клітинку, яка знаходиться в стовпці (рядку), що відповідає найбільшій різниці між двома мінімальними тарифами, які знаходяться в даному стовпці (рядку).
Приклад.
Використовуючи метод апроксимації Фогеля, знайдемо опорний план попередньої задачі.
Таблиця 5.13
|
Пункти відправлення |
Пункти призначення |
Запас |
||||
|
B1 |
В2 |
В3 |
В4 |
В5 |
||
|
А1 |
2 |
3 |
4 |
2 |
4 |
140 |
|
А2 |
8 |
4 |
1 |
4 |
1 |
180 |
|
Аз |
9 |
7 |
3 |
7 |
2 |
160 |
|
Потреби |
60 |
70 |
120 |
130 |
100 |
480 |
Розв'язання.
Для кожного рядка і стовпця таблиці умови знайдемо різницю між двома мінімальними тарифами, записаними в даному рядку або стовпці та помістимо їх у відповідному додатковому стовпці або в додатковому рядку таблиці. Так, в рядку А1 мінімальний тариф, який дорівнює 2, зустрічається двічі, тому різниця дорівнює 0.
Аналогічна ситуація і в рядку А2 з мінімальними тарифами 1. В рядку А3 мінімальні тарифи 3 і 2. Отже, різниця між ними: 3-2=1. Аналогічно отримуємо різницю в стовпці В1: 8-2=6. Вирахувавши всі ці різниці, бачимо, що найбільша з них відповідає стовпцю В1. Таким чином, необхідно заповнити клітинку, що знаходиться на перетині рядка А1 та В1, оскільки в ній найменша ціна. Заповнивши її, ми тим самим задовольнили потреби пункту В1, тому виключаємо з розгляду стовпець В1 і будемо вважати запаси пункту А1 рівними 140-60=80 одиниць. Після цього визначимо наступну клітинку для заповнення. Знову знаходимо різниці між мінімальними тарифами в кожному з рядків і стовпців, що залишились. Запишемо їх, відповідно, в другому додатковому стовпці і рядку таблиці 5.14. Як видно з цієї таблиці, найбільша різниця (в даному випадку однакова) буде в стовпцях В3 і В4. Вибираємо той стовпець, в якому найменший мінімальний тариф. Це буде стовпець В3. На перетині цього стовпця з рядком А2 визначаємо клітинку з мінімальним тарифом і заповнюємо її з врахуванням потреб пункту В3 та запасів пункту А2, виключаємо з розгляду стовпець В3 і будемо вважати запаси пункту А2 рівними 180-120=60 одиниць.
Продовжимо ітераційний процес, послідовно заповнюючи клітинки, що знаходяться на перетині рядка А3 і стовпця В5, рядка А1 та стовпця В4, рядка А2 і стовпця В2, рядка А3 і стовпця В4, рядка А3 і стовпця В2.
Таблиця 5.14
|
Пункти відправлення |
Пункти призначення |
Запаси |
Різниці по рядках |
||||||||
|
B1 |
В2 |
В3 |
В4 |
В5 |
|||||||
|
А1 |
2 60 |
3 |
4 |
2 80 |
4 |
140 |
0 |
1 |
1 |
1 |
0 |
|
А2 |
8 |
4 60 |
1 120 |
4 |
1 |
180 |
0 |
0 |
3 |
0 |
0 |
|
А3 |
9 |
7 10 |
3 |
7 50 |
2 100 |
160 |
1 |
1 |
5(3) |
0 |
0 |
|
Потреби |
60 |
70 |
120 |
130 |
100 |
480 |
|||||
|
Різниці по стовпцях |
6(1) |
1 |
2 |
2 |
1 |
||||||
|
- |
1 |
2(2) |
2 |
1 |
|||||||
|
- |
1 |
- |
2 |
1 |
|||||||
|
- |
1 |
- |
2(4) |
- |
|||||||
|
- |
3(5) |
- |
3 |
- |
В результаті отримаємо опорний план:
.
При цьому плані загальна вартість перевезення така:
.
Як правило, використання методу апроксимації Фогеля позволяв отримати опорний план, близький до оптимального, або сам оптимальний план.
5.2.4 Розв’язання транспортної задачі методом потенціалів
Для визначення оптимального плану транспортної задачі розроблено декілька методів. Але найбільш часто використовуються метод потенціалів і метод диференціальних рент.
Загальний принцип визначення оптимального плану транспортної задачі методом потенціалів аналогічний принципу розв’язання задачі лінійного програмування симплексним методом, тобто: спочатку знаходять опорний план транспортної задачі, а потім його послідовно покращують до отримання оптимального плану. Планом транспортної задачі буде матриця
,
в якій спостерігається баланс по рядках і стовпцях. Це відповідає вимогам обмежень задачі. Опорний план транспортної задачі буде містити (т+п-1) базисний компонент. Система обмежень складається з (т+п) обмежень – рівнянь, однак незалежних рівнянь (т+п-1) (вимога робить систему залежною).
Якщо число базисних компонент менше, ніж (т+п-1), то опорний план вироджений, тоді деякі нульові компоненти необхідно вважати базисними. Їх кількість така, що базисних усього буде (т+п-1).
Опорний план може бути визначений різними способами. Наприклад, метод північно-західного кута, метод найменших елементів, метод Фогеля. Алгоритмічно визначити опорний план – це значить знайти місця розташування додатних компонентів і їх величини таким чином, щоб був баланс по рядках і стовпцях (умова плану) і їх кількість була (т+п-1) (умова опорності плану).
Метод потенціалів заснований на ідеї знаходження потенціалів (величин і ) і виборі компоненти, що вводиться в базис, якщо опорний план не оптимальний.
Алгоритм методу потенціалів являє собою покроковий процес, що скла-
дається з перевірки на оптимальність опорного плану, виродженості і переходу на повний опорний план. Процес керується ознакою оптимальності. Опишемо цей процес.
Крок 1. Звести дані в таблицю і визначити опорний план. Якщо план вироджений, то ввести нульові базисні елементи.
Крок 2. Скласти систему рівнянь для і обчислити всі потенціали, відповідно, пунктів відправки і пунктів призначення і . Ці числа знаходять з системи рівнянь
-=, (5.2.6)
де – тарифи, що стоять в заповнених клітинках таблиці умов транспортної задачі.
Оскільки кількість заповнених клітинок п+т-1, то система (5.2.6) з п+т невідомими містить п+т-1 рівнянь. Оскільки кількість невідомих перевищує на одиницю кількість рівнянь, то одне з невідомих можна покласти рівним довільному числу, наприклад, =0, і знайти послідовно з рівнянь (5.2.6.) значення інших невідомих. Після того, як всі потенціали визначені, для кожної вільної клітинки визначають числа
=--.
Якщо серед чисел немає додатних, то знайдений опорний план є оптимальним. В іншому випадку необхідно перейти до нового опорного плану. Для цього розглядають всі вільні клітинки, для яких 0, і серед даних чисел вибирають максимальне. Клітинку, якій це число відповідає, необхідно заповнити.
Заповнюючи вибрану клітинку, необхідно змінити об’єми поставок, записаних в рядку інших зайнятих клітинок.
Процедура переходу від одного опорного плану до іншого виконується таким чином. Нехай вибраний деякий небазисний елемент , який вводиться в базис. Відмітимо його знаком “+”. Далі будується замкнений ланцюжок, вершинами якого будуть базисні елементи опорного плану. Вершини ланцюжка позначаються знаком “+” і “-” по черзі, дотримуючи правило: кожний рядок і стовпець повинні мати стільки елементів, відмічених знаком “-”, скільки “+”. Серед елементів, відмічених знаком “-”, вибирається мінімальний елемент. Позначимо його через . додається до елементів, відмічених знаком “+”, і віднімається від елементів, відмічених знаком “-”.
Крок 3. Після цієї процедури отриманий новий опорний план перевіряємо на оптимальність, тобто знову повторюємо всі дії кроку 2.
Зауважимо, що в процесі розв’язання задачі може бути отриманий вироджений опорний план. Щоб уникнути в цьому випадку зациклювання, необхідно відповідні нульові елементи опорного плану замінити малим додатним числом Е і розв’язувати задачу як невироджену. В оптимальному плані такої задачі необхідно вважати Е рівним нулеві.
Визначимо оптимальний план задачі 1 методом потенціалів.
В першу чергу, використовуючи метод мінімального елемента, знаходимо опорний план задачі. Цей план записаний в таблицю 5.15.
|
Пункти відправлення |
Пункти призначення |
Запас |
|||
|
B1 |
В2 |
В3 |
В4 |
||
|
А1 |
7 |
8 |
1 160 |
2 |
160 |
|
А2 |
4 120 |
5 |
9 |
8 20 |
140 |
|
Аз |
9
|
2 50 |
3 30 |
6 90 |
170 |
|
Потреби |
120 |
50 |
190 |
110 |
470 |
Складемо систему рівнянь:
;
;
;
;
;
.
Нехай Тоді , .
Тепер для кожної вільної клітинки визначимо числа =--:
Поміщаємо знайдені числа в рамки і записуємо їх відповідно в кожну вільну клітинку таблиці 5.16.
|
Пункти відправлення |
Пункти призначення |
Запас |
|||
|
B1 |
В2 |
В3 |
В4 |
||
|
А1 |
7 -5 |
8 -6 |
1 160 |
2 |
160 |
|
А2 |
4 120 |
5 -5 |
9 -10 |
8 20 |
140 |
|
Аз |
9 -7 |
2 50 |
3 30 |
6 90 |
170 |
|
Потреби |
120 |
50 |
190 |
110 |
470 |
План неоптимальний, оскільки 0. Потрібно компонент ввести в базис. Побудуємо ланцюжок.
|
-160 |
+ |
||
|
120 |
20 |
||
|
50 |
+30 |
-90 |
Найменше з чисел в мінусових клітинках дорівнює 90. Отже, клітинка, в якій знаходиться це число, стає вільною в новій таблиці 5.17. Інші числа в табл.5.1 отримуються так: до числа 30, що стоїть в плюсовій клітинці табл.5.16. додамо 90 і віднімемо 90 від числа 160, що знаходиться в мінусовій клітинці табл. 5.16. Клітинка на перетині рядка А3 і стовпця В4 стає вільною.
Після цих перетворень отримуємо новий опорний план (табл. 5.17) .
|
Пункти відправлення |
Пункти призначення |
Запас |
|||
|
B1 |
В2 |
В3 |
В4 |
||
|
А1 |
7 -2 |
8 -19 |
1 70 |
2 |
160 |
|
А2 |
4 120 |
5 -16 |
9 -16 |
8 20 |
140 |
|
Аз |
9 -8 |
2 50 |
3 120 |
6 -8 |
170 |
|
Потреби |
120 |
50 |
190 |
110 |
470 |
Перевіримо цей план на оптимальність. Знову знаходимо потенціали пунктів відправки і призначення. Для цього складемо таку систему рівнянь:
;
;
;
;
;
.
Покладемо Отримаємо: .
Для кожної вільної клітинки обчислимо числа .
В даному випадку всі 0 і опорний план є оптимальний. Сумарні транспортні витрати:
Отже, отримуємо такий оптимальний план:
VI Градієнтний метод розв’язУВАННЯ задач
нелінійного програмування
В загальному випадку задача нелінійного програмування полягає у визначенні максимального (мінімального) значення функції
(6.1)
за умови, що її змінні задовольняють співвідношення:
(6.2)
де f і gi – деякі відомі функції n змінних, а ві – задані числа.
В результаті розв’язання задачі буде визначена точка , координати якої задовольняють (6.2) і така, що для будь-якої іншої точки , що задовольняє умови (6.2), виконується нерівність:
Якщо f і gi – лінійні функції, то задача (6.1), (6.2) є задачею лінійного програмування.
Співвідношення (6.2) утворюють систему обмежень і містять в собі умови невід’ємності змінних, якщо такі умови існують. Умови додатності змінних можуть бути задані окремо.
Для розв’язання вищенаведеної задачі не існує універсальних методів. Але для окремих класів задач, в яких зроблені додаткові обмеження відносно властивостей функцій f і gi, розроблені ефективні методи їх розв’язання. Так, наприклад, за умови, що f – опукла функція, то область допустимих розв’язків хі≥0, що визначається обмеженнями (6.2),– опукла.
Для кращого розуміння вищесказаного наведемо деякі означення. Нехай Е – (лінійний простір) дійсна вісь, або комплексна площина. Відрізком, що з’єднує точки x1, x2 Е, називається сукупність всіх точок x вигляду:
Множина W в лінійному просторі Е називається опуклою, якщо кожний раз з того, що x1, x2 W, випливає: W належить відрізок, який з’єднує x1 і x2.
Введемо тепер поняття опуклої функції (функціонала). Нехай на лінійному просторі Е задана функція, що ставить кожному хЕ у відповідність f(x). В цьому випадку говорять, що на Е заданий функціонал f(x). Якщо всі значення f(x) дійсні, то функціонал f(x) називається дійсним функціоналом. Дійсний функціонал f(x) називається опуклим на множині xЕ, якщо для будь-яких точок x1 і x2 із x і t [0,1]:
. (6.3)
Якщо нерівність (6.3) виконується навпаки, то функціонал f(x) називається вгнутим.
Задача (6.1), (6.2) називається задачею опуклого програмування, якщо функція f(x) є вгнутою, а функції – опуклі.
Використання градієнтних методів можливе для кожної задачі нелінійного програмування. Але в загальному випадку використання цих методів дозволяє знайти точку локального екстремуму. Тому більш доцільно використовувати градієнтні методи для знаходження розв’язку задач опуклого програмування. Процес знаходження розв’язку задачі полягає в тому, що, починаючи з деякої точки , відбувається послідовний перехід до деяких інших точок до тих пір, поки не виявляється прийнятний розв’язок даної задачі.
Найбільш розповсюдженим методом, при якому точки, що досліджуються, не виходять за межі допустимих розв’язків задачі, є метод Франка-Вульфа.
Процес знаходження розв’язку задачі починають з визначення точки, що належить області допустимих розв’язків задачі. Нехай це точка , тоді в цій точці визначають градієнт функції f
і будують лінійну функцію
. (6.4)
Потім знаходять максимальне значення цієї функції при обмеженні (6.2) і Нехай розв’язок задачі визначають точкою z(k). Тоді за новий допустимий розв’язок приймають координати точки :
(6.5)
де – деяке число (0≤ ≤1). Це число беруть довільно або визначають мінімальне значення з умови df/dк=0. Якщо його значення більше одиниці, то покладають =1. Після визначення числа знаходять координати точки , визначають значення цільової функції в ній і виясняють необхідність переходу до нової точки . Якщо така необхідність є, то визначають в точці градієнт цільової функції, переходять до відповідної задачі лінійного програмування і знаходять її розв’язок.
Визначають координати точки і досліджують необхідність проведення подальших розрахунків.
Приклад
Методом Франка-Вульфа знайти розв’язок задачі, що полягає у визначенні максимального значення функції
(6.6)
за умов:
(6.7)
Розв’язання
Знайдемо градієнт функції
, (6.8)
і за початковий допустимий розв’язок задачі виберемо точку x(0)=(0;0), а за отримуваний розв’язок – нерівність
(6.9)
І ітерація
В точці x(0) градієнт . Знаходимо максимум функції за формулою (6.4): F1=2x1+4x2 за умов (6.5) і (6.6):
(6.10)
(6.11)
Запишемо задачу в канонічній формі і визначимо точку z максимуму функції F1:
(6.12)
(6.13)
Нехай x3, x4 – базисні змінні, а х1, х2 – вільні змінні.
Тоді маємо:
Покладемо тоді цільова функція .
Нехай тоді
Збільшуємо x2 з врахуванням границь x3 і x4:
, тоді x3 = 0;
Отже, при x1= 0, x2= 4 маємо:
Значення цільової функції при цьому буде дорівнювати:
Тепер нехай x3, x4 базисні змінні, а х1, х2 – вільні змінні. Тоді
Значення цільової функції не збільшиться, отже, оптимальний варіант буде в точці
Знайдемо нові допустимі значення розв’язку задачі за формулою:
де
Підставимо замість і їх значення, отримаємо:
(6.14)
Визначимо тепер число з врахуванням (6.5) і (6.14):
Знайдемо похідну цієї функції за і прирівняємо її до нуля:
Маємо =0,25. Оскільки знайдене значення <1, то приймаємо його за величину кроку. Таким чином:
(6.15)
ІІ ітерація.
Градієнт цільової функції задачі в точці Знаходимо максимум функції за умов (6.6) і (6.7).
Тобто, знову розв’язуємо задачу лінійного програмування. Оскільки залежить тільки від х1, то визначимо її з врахуванням того, що х3=х4=0. В даному випадку, розв’язуючи систему рівнянь (6.16)
(6.16)
маємо
Оптимальний план
Визначаємо тепер
Останню рівність перепишемо так:
(6.17)
Підставимо тепер в (6.4) значення (6.16), отримаємо:
звідки . Прирівнюємо до нуля. Тоді 2=0,15. Таким чином,
(6.18)
тобто , ,
.
ІІІ ітерація
Градієнт функції f в точці . Знаходимо максимум функції за умов (6.6) і (6.7), тобто, розв’язуючи задачу лінійного програмування, знаходимо оптимальний план .
Визначимо . Маємо:
Розв’язуючи рівняння, =0, .
Отже, x(3)=(0,99528;0,96321), ,
.
Таким чином, x(3)=(0,99528;0,96321) є розв’язком даної задачі згідно з вибраним значенням . Зменшуючи значення , можна показати, що максимальне значення цільової функції буде в точці x*=(1;1).
Vii Література
PAGE 75
Реальна ситуація
Постановка задачі
Модель
Прогноз
Непротиріччя
Перевірка адекватності
G1
σ
σ
η2
А
В
Збурювальні дії середовища ξ
Технологічна система
(матеріал, процес, апарат, цех)
з власними параметрами хs,i
Вихідні параметри
Уj
Вхідні
фактори
Хі
Блок управління, що формує дії Δхі і Δхs,i
Мета управління
Зовнішнє середовище системи
EMBED Excel.Sheet.8
EMBED Excel.Sheet.8
Т
х
Т2
3
Т4
Т1
Т5
Т
х
Т2
Т3
Т4
Т1
2
1
5
4
3
Т5
5
4
3
2
1
13
10
7
4
1
14
15
2
3
11
12
8
9
5
6
XIII
XIV
XV
XVI
IX
X
XI
XII
V
VII
VI
VIII
I
II
III
IV
y
x
EMBED Equation.3