МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ
НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ
„ЛЬВІВСЬКА ПОЛІТЕХНІКА”
ІКТА
Кафедра ЗІ
Курсова робота
з курсу:
" Комп'ютерні методи дослідження інформаційних процесів і систем "
на тему:
" СЕЛЬСИННА СЛІДКУЮЧА СИСТЕМА "
Тема 5, варіант 10
Зміст
1. Завдання 3
2. Перетворення рівнянь 5
3. Теоретичні відомості 7
3.1 Модифікований метод Ейлера для розв’язку систем диференціальних рівнянь 7
3.2 Метод Рунге-Кутто-Мерсона 10
4. Лістинг програм 11
4.1 Модифікований метод Ейлера 11
4.2 Метод Рунге-Кутто-Мерсона з автоматичною зміною кроку 14
5. Результати виконання програми 17
5.1 Модифікований метод Ейлера 17
5.2 Метод Рунге-Кутто-Мерсона з автоматичною зміною кроку 18
6. Графік перехідного процесу 20
7. Список літератури 21
1. Завдання
Схема:
Рівняння ланок :
вимірювальна схема
електронний підсилювач
обмотка збудження ЕМП (електромашинного підсилювача)
двигун
редуктор
короткозамкнута обмотка ЕМП
При початкових параметрах
Параметри
Tm. (сек)
0.2
Tk (сек)
0.02
TI (сек)
0.005
C (рад/в.сек)
20
i
400
KI
2
K2
3
KЕП
2
S (в/рад)
50
Звести систему алгебро-диференціальних рівнянь до системи чотирьох диференціальних рівнянь першого порядку, представити її у нормальній формі та розв’язати цю систему вказаними методами. Початкові умови - =1 радіан, решта початкових умов – нульові. Числові значення сталих параметрів, заданих в таблиці, слід зобразити з допомогою одиниць системи СІ.
Побудувати графік зміни величини
2. Перетворення рівнянь
U = S ( - )
+ = ( - )
+ = C =
= i
+ = =
=
=
+ = C =
=
= ; = ; ; = ;
= K1; = KU; S = S; = T1;
= K2; = TK; C = C; = TM;
= II;
(0) = 0; (0) = 0;(0) = 0;(0) = 0 ; = 1
= F[1]; - Y[1];
= F[2]; - Y[2];
= F[3]; - Y[3];
= F[4]; - Y[4];
3. Теоретичні відомості
3.1 Модифікований метод Ейлера для розв’язку систем диференціальних рівнянь
Диференціальне рівняння (ДР), що містить лише одну незалежну змінну і похідні за нею, називають звичайними (ДР). ДР, що містить декілька незалежних змінних і похідні за ними, називають рівняння в частинних похідних.
Порядком ДР називається найвищий порядок похідної (або диференціалу), який входить в рівняння. Звичайне ДР (ЗДР) -го порядку в загальному випадку має незалежну змінну, невідому функцію та її похідні до -го порядку включно:
- незалежна змінна;
- невідома функція (залежна змінна);
- похідні цієї функції.
Диференціальне рівняння -го порядку, розв’язане відносно старшої похідної, може бути записано у вигляді:
Щоб розв’язати ЗДР, необхідно мати значення залежної змінної та (або) її похідних при деяких значення незалежної змінної.
Метод Ейлера є найпростішим методом розв’язування задачі Коші. Він дозволяє інтегрувати ДР першого порядку виду.
(4)
Метод Ейлера базується на розкладі функції в ряд Тейлора в околі точки
(5)
Якщо мале, то, знехтувавши членам розкладу, що містять в собі і т.д. отримаємо
(6)
Похідну знаходимо з рівняння (4), підставивши в нього початкову умову. Таким чином можна знайти наближене значення залежної змінної при малому зміщенні від початкової точки. Цей процес можна продовжувати, використовуючи співвідношення.
,
роблячи як завгодно багато кроків.
Похибка методу має порядок , оскільки відкинуті члени, що містять в другій і вище степенях.
Недолік методу Ейлера - нагромадження похибок, а також збільшення об’ємів обчислень при виборі малого кроку з метою забезпечення заданої точності.
В методі Ейлера на всьому інтервалі тангенс кута нахилу дотичної приймається незмінним і рівним . Очевидно, що це призводить до похибки, оскільки кути нахилу дотичної в точках та різні. Точність методу можна підвищити, якщо покращити апроксимацію похідної.
Модифікований метод Ейлера з перерахунком
Обчислення за методом Ейлера з перерахунком робляться в два етапи.
Прогноз:
.
Корекція:
.
Модифікований метод Ейлера з перерахунком має другий порядок точності, проте для його реалізації необхідно двічі обчислювати праву частину функції. Зауважимо, що метод Ейлера з перерахунком являє собою різновид методів Рунге-Кутта (предиктор-коректор).
Метод Ейлера з автоматичною зміною кроку
Після обчислення з кроком всі обчислення виконуються повторно з кроком . Після цього порівнюються результати, отримані в точці хn+1 з кроком і . Якщо модуль різниці менший , то обчислення продовжуються з кроком , в іншому випадку крок зменшують. Якщо нерівність дуже сильна, то крок збільшують.
- значення функції в точці
- значення функції в точці , обчислене з кроком
- значення функції в точці , обчислене з кроком
- значення функції , обчислене з кроком
1) Якщо
То,
обчислення повторюються доки не виконається умова .
2) Якщо виконується ця умова, то можливі два варіанти, в залежності від значення K, де K – ознака поділу кроку.
Початкове значенняі залишається таким після першого поділу кроку на два. Надалі, якщо крок ділиться, то K приймає значення одиниці.
а) Якщо , то навіть коли виконалась умова , крок не змінюється, тобто лишається тим самим (обчислення далі проводяться з попереднім кроком).
б) Якщо і виконалась умова , тоді .
В обох випадках а) і б) результат виводиться на друк.
Для розв’язку системи диференціальних рівнянь використовують цей самий метод, за виключенням того, що всі рівняння системи необхідно розв’язувати паралельно.
3.2 Метод Рунге-Кутто-Мерсона з автоматичною зміною кроку
Метод дозволяє оцінити похибку на кожному кроці інтегрування. При цьому не потрібно зберігати в пам’яті обчислення значень функцій на кроці і для оцінки похибки.
Алгоритм методу
1. Задаємо число рівнянь , похибку , початковий крок інтегрування , початкові умови .
2. За допомогою п’яти циклів з керуючою змінною обчислюємо коефіцієнти
3. Знаходимо значення
та похибку
4. Перевіряємо виконання умов
Можливі випадки:
а) Якщо перша умова не виконується, тобто , то ділимо крок на 2 та повторюємо обчислення з п.2, встановивши початкові значення .
б) Якщо виконується перша та друга умови, значення та виводяться на друк.
Якщо друга умова не виконується, крок збільшується вдвічі і тоді обчислення знову повторюється з п.2.
Треба відмітити, що похибка на кожному кроці методу Рунге-Кутта-Мерсона оцінюється приблизно. При розв’язуванні нелінійних ДР істинна похибка може відрізнятися в декілька разів від заданої .
, де .
- крок поділити на 2 і повернутися на початок.
для всіх рівнянь: виводимо на друк , а крок збільшуємо удвічі.
4. Лістинг програм
4.1 Модифікований метод Ейлера
using System;
using System.IO;
namespace RK
{
class Data
{
int Kr = 0;
double[,] Y = new double[4, 10000000];
double[,] y = new double[4, 10000000];
double[] t = new double[10000000];
double[] m = new double[6];
double[] n = new double[6];
double[] l = new double[6];
double[] p = new double[6];
double dt;
double k1;
double Kp;
double S;
double Tl;
double k2;
double Tk;
double C;
double Tm;
double i;
long k;
int Qin;
double y0;
double e;
double b;
double d;
double z1;
double z2;
double z3;
double z4;
double z5;
double z6;
double z7;
double z8;
public Data()
{
Kp = 2;
Tm = 0.2;
Tk = 0.02;
Tl = 0.005;
C = 15;
i = 300;
k1 = 2;
k2 = 3;
S = 50;
Qin = 1;
dt = 0.001;
e = 0.0001;
b = 2;
t[0] = 0.0;
}
public double F1(double t, double y1, double y2, double y3, double y4)
{
return (k1 * Kp * S * (Qin - y2) - y1) / Tl;
}
public double F2(double t, double y1, double y2, double y3, double y4)
{
return y4 / i;
}
public double F3(double t, double y1, double y2, double y3, double y4)
{
return (k2 * y1 - y3) / Tk;
}
public double F4(double t, double y1, double y2, double y3, double y4)
{
return (C * y3 - y4) / Tm;
}
public void Pohidni()
{
Y[0, k] = F1(t[k], y[0, k], y[1, k], y[2, k], y[3, k]);
Y[1, k] = F2(t[k], y[0, k], y[1, k], y[2, k], y[3, k]);
Y[2, k] = F3(t[k], y[0, k], y[1, k], y[2, k], y[3, k]);
Y[3, k] = F4(t[k], y[0, k], y[1, k], y[2, k], y[3, k]);
}
public void Znach(double x, double y1, double y2, double y3, double y4)
{
y[0, k + 1] = y[0, k] + dt * F1(x + dt / 2.0, y1, y2, y3, y4);
y[1, k + 1] = y[1, k] + dt * F2(x + dt / 2.0, y1, y2, y3, y4);
y[2, k + 1] = y[2, k] + dt * F3(x + dt / 2.0, y1, y2, y3, y4);
y[3, k + 1] = y[3, k] + dt * F4(x + dt / 2.0, y1, y2, y3, y4);
}
public void Mod(double x, double y1, double y2, double y3, double y4)
{
Kr = 0;
Znach(x, y1, y2, y3, y4);
z1 = y[0, k + 1];
z2 = y[1, k + 1];
z3 = y[2, k + 1];
z4 = y[3, k + 1];
z5 = y[0, k] + dt * (F1(x, y1, y2, y3, y4) + F1(x + dt, z1, z2, z3, z4)) / 2.0;
z6 = y[1, k] + dt * (F2(x, y1, y2, y3, y4) + F2(x + dt, z1, z2, z3, z4)) / 2.0;
z7 = y[2, k] + dt * (F3(x, y1, y2, y3, y4) + F3(x + dt, z1, z2, z3, z4)) / 2.0;
z8 = y[3, k] + dt * (F4(x, y1, y2, y3, y4) + F4(x + dt, z1, z2, z3, z4)) / 2.0;
while (Math.Abs(z1 - z5) > e)
{
z1 = z5;
z2 = z6;
z3 = z7;
z4 = z8;
z5 = y[0, k] + dt * (F1(x, y1, y2, y3, y4) + F1(x + dt, z1, z2, z3, z4)) / 2.0;
z6 = y[1, k] + dt * (F2(x, y1, y2, y3, y4) + F2(x + dt, z1, z2, z3, z4)) / 2.0;
z7 = y[2, k] + dt * (F3(x, y1, y2, y3, y4) + F3(x + dt, z1, z2, z3, z4)) / 2.0;
z8 = y[3, k] + dt * (F4(x, y1, y2, y3, y4) + F4(x + dt, z1, z2, z3, z4)) / 2.0;
if (Kr % 4 == 0)
dt = dt / 2.0;
Kr++;
}
y[0, k + 1] = z5;
y[1, k + 1] = z6;
y[2, k + 1] = z7;
y[3, k + 1] = z8;
}
public void Prod()
{
do
{
Mod(t[k], y[0, k], y[1, k], y[2, k], y[3, k]);
t[k + 1] = t[k] + dt;
k++;
} while (t[k] < b);
Pohidni(t[k], y[0, k], y[1, k], y[2, k], y[3, k]);
}
public void vuvid()
{
int j;
Random rand = new Random();
StreamWriter log_out;
log_out = new StreamWriter("logfile.txt");
Console.SetOut(log_out);
Console.WriteLine("t\ty1'");
for (j = 0; j < k; j += rand.Next(1, 15))
{
Console.WriteLine("{0:0.####}\t\t{1:0.####}", t[j], y[2, j]);
}
log_out.Close();
}
public void vuv()
{
Console.WriteLine("t\ty1\ty1'");
for (int j = 0; j < k; j += 20)
{
Console.WriteLine("{0:0.####}\t\t\t{1:0.####}", t[j], y[2, j]);
}
}
}
class Program
{
static void Main(string[] args)
{
Data d = new Data();
d.Prod();
d.vuv();
d.vuvid();
Console.ReadLine();
}
}
}
4.2 Метод Рунге-Кутта-Мерсона
using System;
using System.IO;
namespace RK
{
class Data
{
double[,] Yp = new double[4, 10000000];
double[,] y = new double[4, 10000000];
double[] time = new double[10000000];
double[] t = new double[10000000];
double[] m = new double[5];
double[] n = new double[5];
double[] l = new double[5];
double[] s = new double[5];
double dt;
double k1;
double Kp;
double S;
double Tl;
double k2;
double Tk;
double C;
double Tm;
double i;
long k;
int Qin;
double e;
double b;
double z1;
double z2;
double z0;
public Data()
{
Kp = 2;
Tm = 0.2;
Tk = 0.02;
Tl = 0.005;
C = 20;
i = 400;
k1 = 2;
k2 = 3;
S = 50;
Qin = 1;
dt = 0.001;
e = 0.0001;
b = 2;
time[0] = 0.0;
}
public double F1(double t, double y1, double y2, double y3, double y4)
{
return (k1 * Kp * S * (Qin - y2) - y1) / Tl;
}
public double F2(double t, double y1, double y2, double y3, double y4)
{
return y4 / i;
}
public double F3(double t, double y1, double y2, double y3, double y4)
{
return (k2 * y1 - y3) / Tk;
}
public double F4(double t, double y1, double y2, double y3, double y4)
{
return (C * y3 - y4) / Tm;
}
public void Pohidni(double t, double y1, double y2, double y3, double y4)
{
Yp[0, k] = F1(t, y1, y2, y3, y4);
Yp[1, k] = F2(t, y1, y2, y3, y4);
Yp[2, k] = F3(t, y1, y2, y3, y4);
Yp[3, k] = F4(t, y1, y2, y3, y4);
}
public void Zmin()
{
y[0, k + 1] = y[0, k] + (m[0] + 4.0 * m[3] + m[4]) / 6.0;
y[1, k + 1] = y[1, k] + (n[0] + 4.0 * n[3] + n[4]) / 6.0;
y[2, k + 1] = y[2, k] + (l[0] + 4.0 * l[3] + l[4]) / 6.0;
y[3, k + 1] = y[3, k] + (s[0] + 4.0 * s[3] + s[4]) / 6.0;
}
public void Kof(double t, double y1, double y2, double y3, double y4, double h)
{
m[0] = F1(t, y1, y2, y3, y4) * h;
n[0] = F2(t, y1, y2, y3, y4) * h;
l[0] = F3(t, y1, y2, y3, y4) * h;
s[0] = F4(t, y1, y2, y3, y4) * h;
m[1] = F1(t + h / 3.0, y1 + m[0] / 3.0, y2 + n[0] / 3.0, y3 + l[0] / 3.0, y4 + s[0] / 3.0) * h;
n[1] = F2(t + h / 3.0, y1 + m[0] / 3.0, y2 + n[0] / 3.0, y3 + l[0] / 3.0, y4 + s[0] / 3.0) * h;
l[1] = F3(t + h / 3.0, y1 + m[0] / 3.0, y2 + n[0] / 3.0, y3 + l[0] / 3.0, y4 + s[0] / 3.0) * h;
s[1] = F4(t + h / 3.0, y1 + m[0] / 3.0, y2 + n[0] / 3.0, y3 + l[0] / 3.0, y4 + s[0] / 3.0) * h;
m[2] = F1(t + h / 3.0, y1 + m[1] / 6.0 + m[0] / 6.0, y2 + n[1] / 6.0 + n[0] / 6.0, y3 + l[1] / 6.0 + l[0] / 6.0, y4 + s[1] / 6.0 + s[0] / 6.0) * h;
n[2] = F2(t + h / 3.0, y1 + m[1] / 6.0 + m[0] / 6.0, y2 + n[1] / 6.0 + n[0] / 6.0, y3 + l[1] / 6.0 + l[0] / 6.0, y4 + s[1] / 6.0 + s[0] / 6.0) * h;
l[2] = F3(t + h / 3.0, y1 + m[1] / 6.0 + m[0] / 6.0, y2 + n[1] / 6.0 + n[0] / 6.0, y3 + l[1] / 6.0 + l[0] / 6.0, y4 + s[1] / 6.0 + s[0] / 6.0) * h;
s[2] = F4(t + h / 3.0, y1 + m[1] / 6.0 + m[0] / 6.0, y2 + n[1] / 6.0 + n[0] / 6.0, y3 + l[1] / 6.0 + l[0] / 6.0, y4 + s[1] / 6.0 + s[0] / 6.0) * h;
m[3] = F1(t + h / 2.0, y1 + 3.0 * m[2] / 8.0 + m[0] / 8.0, y2 + 3.0 * n[2] / 8.0 + n[0] / 8.0, y3 + 3.0 * l[2] / 8.0 + l[0] / 8.0, y4 + 3.0 * s[2] / 8.0 + s[0] / 8.0) * h;
n[3] = F2(t + h / 2.0, y1 + 3.0 * m[2] / 8.0 + m[0] / 8.0, y2 + 3.0 * n[2] / 8.0 + n[0] / 8.0, y3 + 3.0 * l[2] / 8.0 + l[0] / 8.0, y4 + 3.0 * s[2] / 8.0 + s[0] / 8.0) * h;
l[3] = F3(t + h / 2.0, y1 + 3.0 * m[2] / 8.0 + m[0] / 8.0, y2 + 3.0 * n[2] / 8.0 + n[0] / 8.0, y3 + 3.0 * l[2] / 8.0 + l[0] / 8.0, y4 + 3.0 * s[2] / 8.0 + s[0] / 8.0) * h;
s[3] = F4(t + h / 2.0, y1 + 3.0 * m[2] / 8.0 + m[0] / 8.0, y2 + 3.0 * n[2] / 8.0 + n[0] / 8.0, y3 + 3.0 * l[2] / 8.0 + l[0] / 8.0, y4 + 3.0 * s[2] / 8.0 + s[0] / 8.0) * h;
m[4] = F1(t + h, y1 + m[0] / 2.0 - 3.0 * m[2] / 2.0 + 2.0 * m[3], y2 + n[0] / 2.0 - 3.0 * n[2] / 2.0 + 2.0 * n[3], y3 + l[0] / 2.0 - 3.0 * l[2] / 2.0 + 2.0 * l[3], y4 + s[0] / 2.0 - 3.0 * s[2] / 2.0 + 2.0 * s[3]) * h;
n[4] = F2(t + h, y1 + m[0] / 2.0 - 3.0 * m[2] / 2.0 + 2.0 * m[3], y2 + n[0] / 2.0 - 3.0 * n[2] / 2.0 + 2.0 * n[3], y3 + l[0] / 2.0 - 3.0 * l[2] / 2.0 + 2.0 * l[3], y4 + s[0] / 2.0 - 3.0 * s[2] / 2.0 + 2.0 * s[3]) * h;
l[4] = F3(t + h, y1 + m[0] / 2.0 - 3.0 * m[2] / 2.0 + 2.0 * m[3], y2 + n[0] / 2.0 - 3.0 * n[2] / 2.0 + 2.0 * n[3], y3 + l[0] / 2.0 - 3.0 * l[2] / 2.0 + 2.0 * l[3], y4 + s[0] / 2.0 - 3.0 * s[2] / 2.0 + 2.0 * s[3]) * h;
s[4] = F4(t + h, y1 + m[0] / 2.0 - 3.0 * m[2] / 2.0 + 2.0 * m[3], y2 + n[0] / 2.0 - 3.0 * n[2] / 2.0 + 2.0 * n[3], y3 + l[0] / 2.0 - 3.0 * l[2] / 2.0 + 2.0 * l[3], y4 + s[0] / 2.0 - 3.0 * s[2] / 2.0 + 2.0 * s[3]) * h;
}
public void Zminna(double t, double y1, double y2, double y3, double y4, double h)
{
Kof(t, y1, y2, y3, y4, h);
Zmin();
}
public void Perevirka(double t, double y1, double y2, double y3, double y4, double h)
{
A1: Zminna(t, y1, y2, y3, y4, dt);
z0 = y[2, k + 1];
z1 = y[2, k] + (l[0] - 3.0 * l[2] + 4.0 * l[3]) / 2.0;
z2 = Math.Abs(z0 - z1) / 2.0;
if (z1 > e)
{
dt = dt / 2.0;
goto A1;
}
if (z1 < e / 30.0)
{
dt = dt * 2.0;
goto A1;
}
y[2, k + 1] = z0;
}
public void Prod()
{
k = 1;
do
{
Pohidni(t[k], y[0, k], y[1, k], y[2,