Physical Process Modeling BG

 

Числени методи за решаване на системи диференциални уравнения

 

 

 

Темата разглежда методите за решаване на диференциални уравнения (ДУ) от първи ред с начални условия:

 

(1)

(2)

 

Методиките за числено решаване на ДУ се прилагат за системи ДУ от първи ред, което налага уравнения от по – висок ред да се сведат до системи уравнения от първи ред. Например ДУ от втори ред:

 

 

се представя като:

 

 

 

където z е нова функция на x, дефинирана от второто уравнение. Това е система от две уравнения за функциите y(x) и z(x), като нейното решение дава едновременно у и производната и.

Механизма на численото решаване на уравнения (1) и (2) се свежда до следното: уравнение (1) дефинира непълно крива в равнината Оxy и дава стойността на производната във всяка точка от кривата като функция на x и y. Уравнение (2) уточнява конкретна крива от фамилията криви удовлетворяващи (1), минаваща през дадена точка. Решаването на задача (1) – (2) представлява функция на х. За да бъдат намерени нейните стойности трябва да се зададат конкретни стойности в израза за x и да се пресметнат стойностите на у.

В най – общия смисъл численото решаване протича по следният ред:

-                        ДУ задава наклона на кривата в произволна точка като функция на x и у. От условието е известна само една точка, т.е. точката с координати (x0, y0)

-                        От тази точка се пресмята наклона на кривата х = х0 и у = у0.

-                        Определя се стъпка на предвижване (или такава е зададена предварително) h, като координатите се предвижват спрямо нея по допирателната.

-                        В новата точка се определя абсциса х1 = х0 + х и ордината y = y1

-                        Определят се следващите стъпки с приближение h.

Описаният метод е известен още като метода на Ойлер. Характерно за него е известна неустойчивост на решението изразяваща се във възможността за получаване на резултат с недопустима грешка. Този проблем може да се избегне при отчитане на кривината на точното решение. Методите за това се делят на едностъпкови и многостъпкови

Едностъпковите методи се основават на използване на информацията за кривата в една точка и не извършват итерации. Към този метод се причислява методът на Тейлър, който не е удобен за употреба. С цел практическа реализация на числените методи за решаване на ДУ по-подробно е разгледан методът на Рунге – Кута. Тази група методи са директни – не използват итерации, поради което изглеждат икономични, но в действителност изискват многократно пресмятане на функцията. Основен недостатък е трудното осъществяване на оценка на допуснатата грешка.

При Многостъпковите методи намирането на следващата точка става с по-малко пресмятания, но изискват итерации за постигане на по-голяма точност. В голямата си част тези методи се наричат предикторно – коректорни методи. Характерни за многостъпковите методи са някои трудности при организирането на итерациите, но основно предимство са получените оценки за допуснатата грешка заедно с резултатите.

 

Метод на Тейлър

 

Методът на Тейлър теоретически дава решението на произволно ДУ, но представлява малък интерес за изчислителната практика. Неговата важност се състои в това, че той е база за сравняване и оценяване на различните практически удобни методи.

Удобно е излагането на метода да започне с допускане, че вече е получено приближение у(х) в интервала х0 ≤ х ≥ хm Формулата на Тейлър за функцията y(x) в точката х = хm :

 

(3)

 

 

 

Където ym(j) e j–тата производна на у(х) в точката х = хm a ξ e в интервала

 

 

 

При заместване х в (3) с хm + h се получава приближена стойност на у(х) в точката хm+1 = xm + h :

 

 

 

 

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

 

 

Диференцирайки относно х:

 

 

 

 

Където частните производни по х и у са означени с f със съответен долен индекс:

 

 

 

Като пример ще бъде разгледан случай при който j = 2. Получава се:

 

 

 

 

 

 

 

При пренебрегване на последния член и пресмятане на уm+1:

 

 

 

 

Грешка от прекъсване:

 

 

 

Ако третата производна не се мени много се записва:

 

 

Където К е константа.

Както беше отбелязано по-горе методът е разгледан с допускане на някакво произволно приближение. Стойността на h може да се намери по следният начин: полага се m = 0 с което се пресмята у1; от това се определя приближената стойност на решението в точката х = х0 + h; полага се m = 1 и с получената стойност у1 и х1 = х0 + h се пресмята y2; по същият начин се определят y3, y4, ... , ym, ym+1 ... Грешката от прекъсване се натрупва на всяка стъпка.

От практическа гледна точка методът на Тейлор е свързан с някои трудности при намирането на fx и fy. По - малка грешка се получава при пресмятане на уm’’’ която е равна на:

 

 

 

Методи на Рунге - Кута

 

Методите на Рунге – Кута притежават следните отличителни свойства:

1.                                                        Те са едностъпкови: за определяне на уm+1 e необходима информация само за предишната точка (хm, ym).

2.                                                        Съгласуват се с тейлоровото развитие до члена с hp, където степента p е различна за различните методи и се нарича ред на метода.

3.                                                        Изискват само пресмятане на функцията f(x,y), но не и на нейни производни.

 

Методът на Ойлер се използва рядко, но разглеждането му е отправна точка, за по-нататъшното разглеждане на методите от този клас. Графичното представяне е показано на фигура 1.

 

фигура 1

Геометрично представяне на метода на Ойлер

 

Известна е точка с координати (хm, ym) лежаща на търсената крива. През нея е прекарана крива с наклон

 

 

Може да се положи ym+1 равно на ординатата на пресечната точка на L1 и правата x = xm+1 = x+h. Уравнението на правата L1 е:

 

 

и понеже

 

 

то

 

 

Грешката при х=хm+1 e означена на фигурата с e.

Последната формула задава метода на Ойлер. Характерни за него са голямата грешка от прекъсване и неустойчивост в някой случаи т.е. малка грешка от закръгляне нараства с увеличаването на х.

Модифицираният метод на Ойлер се основава на намирането на средно аритметичното на наклоните на допирателните в точките (хm, ym) и (xm+h, ym+hym). Графично методът е представен на фигура 2.

 

фигура. 2

Геометрично представяне на модифицирания метод на Ойлер

 

Разглежданият метод може да се опише със следният алгоритъм:

1.                                                                                   Намира се тока с координати (хm, ym) и (xm+h, ym+hym), лежаща на права L1.

2.                                                                                   Пресмята се наклонът на кривата през същата точка посредством намиране на стойността на f.

3.                                                                                   От извършените изчисления се намира правата L2.

4.                                                                                   Средно аритметично на двата наклона се определя пунктираната права L

5.                                                                                   Пресечната точка на тази права с вертикалната права x=xm+1 = xm+h e търсената точка (хm+1, ym+1)

 

Наклонът на правата L e

 

 

 

където:

 

 

Уравнението за L се задава като:

 

 

 

 

Последните уравнения изразяват модифицирания метод на Ойлер.

 

При обобщаването на изложените методи може да се изведе методът на Рунге – Кута за решаване на уравнения от трети и четвърти ред. Класическият метод се задава с равенствата:

 

 

 

където:

 

 

 

 

 

 

 

 

 

Грешката от прекъсване е:

 

 

 

където К е константа.

 

Методи на Адамс – Башфорд и Адамс – Маултън

 

Двата вида методи са базирани върху формулата на Нютон за интерполиране назад. Методите от първи, втори и трети порядък изглеждат така:

 

Адамс – Башфорд

Адамс – Маултън

 

Методите са от явен тип и за да стартират имат нужда от акумулирано количество решения за минали моменти от време.

 

Методи на Гриар

 

 

Методите на Рунге – Кута са едностъпкови и самостартиращи се, което ги прави удобни като предиктор или за натрупване на комплект начални решения преди да се приложи някой многостъпков метод за интегриране. При тях трудно може да се предскаже големината на стъпката, поради което промяната порядъка на интегриране в хода на решаване на диференциалните уравнения се прави трудно. При решаване на системи ДУ тези методи увеличават изключително много количеството изчислителна работа. При методите на Адамс – Бошфорд, Адамс – Маултън и Гриар сравнително лесно се определя големината на стъпката, а с това има възможност за лесна промяна на порядъка на интегриране в хода на изчислителният процес. При решаване на системи ДУ тези методи са удобни за работа, защото не увеличават несъразмерно количеството изчислителна работа.

 

В настоящата разработка “математическо моделиране на електротехнологични процеси”, системата диференциални уравнения се представя във нормализирана форма:

 

 

където [x(t)] е квадратна матрица с неизвестни функции, а [f[x(t)],t] е матрица – стълб съдържаща нелинейните функции на неизвестните величини. В програма Matlab са реализирани функции за числено решаване на системи диференциални уравнения във файлове ode15s.m, ode45.m, ode23.m, ode23s.m, ode113.m. Всички получени резултати предложени при разглеждането на различните модели и регулиране са получени чрез използването им.

Освен чрез посочените функции, след съставянето на системите диференциални уравнения поставените задачите на проекта се реализират и чрез функциите за намиране на собствени вектори на матрици, матрични експоненти, метод на променливите на състоянието.

Собствените числа и собствените вектори са в основата на редица числени алгоритми за изчисляване на матрични експоненти определяне на системите ДУ, определяне на решенията на системите ДУ във временната или честотната област и др. Това определя задачата като една от най-важните в изчислителната математика. Собствените числа lI и собствените вектори [xi] на квадратната матрица [A] се определят от следното уравнение:

 

 

В програма Matlab основните функции са:

eig (A) – връща като резултат вектор със собствените числа на квадратна матрица [A]

[V,D] = eig(A) – произвежда диагонална матрица [D], съдържаща собствените числа на матрицата [A] и матрицата [V], чиито колони съдържат компонентите на собствените вектори на матрицата [A]. Матричната експонента:

 

 

представлява специален тип матрица, която съдържа на мястото на всеки свой елемент група от експоненциални функции. Функцията за изчисляване на матрични експоненти е expm(A), която използва апроксимацията на Паде.

Съставените уравнения се решават по формулата на Коши:

 

 

[x(t)] – променливи на състоянието;

[f(t)] – входни въздействия;

[y(t)] – включва всички величини които не са променливи на състоянието

[A], [B], [C], [D] – матрици, съдържащи параметрите на елементите на веригата

 

 

 

Physical Process Modeling BG