В даній роботі розглянено застосування методу Рунге – Кутта та метод Ейлера для дослідження перехідного процесу сельсинної слідкуючою системи. Розв’язок поставленої задачі представлений в середовищі С# і платформі Visual Studio 2010. Графіки уточнень побудовані в середовищі Microsoft Excel.
Зміст
Постановка задачі...................................4
Перетворення рівнянь................................6
Теоретичні відомості................................8
3.1 Метод Рунге – Кутта для розв’язку систем диференціальних рівнянь............................8
3.2 Метод Ейлера...................................12
Лістинг програми....................................14
Результати виконання програм........................22
6.Графіки перехідного процесу.........................29
7.Висновок............................................28
7.Список літератури...................................28
Постановка задачі
СЕЛЬСИННА СЛІДКУЮЧА СИСТЕМА
Схема:
Рівняння ланок :
вимірювальна схема
електронний підсилювач
обмотка збудження ЕМП (електромашинного підсилювача)
двигун
редуктор
короткозамкнута обмотка ЕМП
При початкових параметрах
Параметри
11
Tm. (сек)
0,3
Tk (сек)
0,02
TI (сек)
0,006
C (рад/в.сек)
4
I
350
KI
2
K2
2
KЕП
9
S (в/рад)
60
Звести систему алгебро-диференціальних рівнянь до системи чотирьох диференціальних рівнянь першого порядку, представити її у нормальній формі та розв’язати цю систему вказаними методами. Початкові умови - =1 радіан, решта початкових умов – нульові. Числові значення сталих параметрів, заданих в таблиці, слід зобразити з допомогою одиниць системи СІ.
Побудувати графік зміни величини
Перетворення рівнянь
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. Метод Рунге – Кутта для розв’язку систем
диференціальних рівнянь.
Диференціальне рівняння (ДР), що містить лише одну незалежну змінну і похідні за нею, називають звичайними (ДР). ДР, що містить декілька незалежних змінних і похідні за ними, називають рівняння в частинних похідних.
Порядком ДР називається найвищий порядок похідної (або диференціалу), який входить в рівняння. Звичайне ДР (ЗДР) -го порядку в загальному випадку має незалежну змінну, невідому функцію та її похідні до -го порядку включно:
- незалежна змінна;
- невідома функція (залежна змінна);
- похідні цієї функції.
Диференціальне рівняння -го порядку, розв’язане відносно старшої похідної, може бути записано у вигляді:
Щоб розв’язати ЗДР, необхідно мати значення залежної змінної та (або) її похідних при деяких значення незалежної змінної.
В методі Рунге-Кутта значення функції визначається за формулою
Якщо розкласти функцію в ряд Тейлора і обмежитись членами до включно, то приріст можна записати у вигляді
(1)
Замість того, щоб обчислювати члени ряду за формулою (1) в методі Рунге-Кутта використовують наступні формули.
Похибка на кожному кроці має порядок . Таким чином метод Рунге-Кутта забезпечує високу точність, однак вимагає більшого об’єму обчислень.
Деколи зустрічається інша форма представлення методу Рунге-Кутта 4-го порядку точності.
Методи з автоматичною зміною кроку
Застосовуються в тому випадку, якщо розв’язок потрібно одержати із заданою точністю. При високій точності (похибка ) автоматична зміна кроку забезпечує зменшення загального числа кроків в декілька разів (особливо при розв’язках у вигляді кривих, що сильно відрізняються крутизною).
Метод Рунге-Кутта з автоматичною зміною кроку
Після обчислення з кроком всі обчислення виконуються повторно з кроком . Після цього порівнюються результати, отримані в точці хn+1 з кроком і . Якщо модуль різниці менший , то обчислення продовжуються з кроком , в іншому випадку крок зменшують. Якщо нерівність дуже сильна, то крок збільшують.
Маємо
- значення незалежної змінної в точці
- значення функції в точці
- значення функції в точці , обчислене з кроком
- значення функції в точці , обчислене з кроком
- значення функції , обчислене з кроком
1) Якщо
обчислення повторюються з кроком і т.д., доки не виконається умова .
2) Якщо виконується ця умова, то можливі два варіанти, в залежності від значення K, де K – ознака поділу кроку.
Початкове значенняі залишається таким після першого поділу кроку на два. Надалі, якщо крок ділиться, то K приймає значення одиниці.
а) Якщо , то навіть коли виконалась умова , крок не змінюється, тобто лишається тим самим (обчислення далі проводяться з попереднім кроком).
б) Якщо і виконалась умова , тоді .
В обох випадках а) і б) результат виводиться на друк.
Для розв’язку системи диференціальних рівнянь використовують цей самий метод, за виключенням того, що всі рівняння системи необхідно розв’язувати паралельно.
3.2 Метод Ейлера
для розв’язку систем диференціальних рівнянь
Метод Ейлера є найпростішим методом розв’язування задачі Коші. Він дозволяє інтегрувати ДР першого порядку виду.
(4)
Метод Ейлера базується на розкладі функції в ряд Тейлора в околі точки
(5)
Якщо мале, то, знехтувавши членам розкладу, що містять в собі і т.д. отримаємо
(6)
Похідну знаходимо з рівняння (4), підставивши в нього початкову умову. Таким чином можна знайти наближене значення залежної змінної при малому зміщенні від початкової точки. Цей процес можна продовжувати, використовуючи співвідношення.
,
роблячи як завгодно багато кроків.
Похибка методу має порядок , оскільки відкинуті члени, що містять в другій і вище степенях.
Недолік методу Ейлера - нагромадження похибок, а також збільшення об’ємів обчислень при виборі малого кроку з метою забезпечення заданої точності.
В методі Ейлера на всьому інтервалі тангенс кута нахилу дотичної приймається незмінним і рівним . Очевидно, що це призводить до похибки, оскільки кути нахилу дотичної в точках та різні. Точність методу можна підвищити, якщо покращити апроксимацію похідної.
Метод Ейлера з автоматичною зміною кроку
Після обчислення з кроком всі обчислення виконуються повторно з кроком . Після цього порівнюються результати, отримані в точці хn+1 з кроком і . Якщо модуль різниці менший , то обчислення продовжуються з кроком , в іншому випадку крок зменшують. Якщо нерівність дуже сильна, то крок збільшують.
- значення функції в точці
- значення функції в точці , обчислене з кроком
- значення функції в точці , обчислене з кроком
- значення функції , обчислене з кроком
1) Якщо
То,
обчислення повторюються доки не виконається умова .
2) Якщо виконується ця умова, то можливі два варіанти, в залежності від значення K, де K – ознака поділу кроку.
Початкове значенняі залишається таким після першого поділу кроку на два. Надалі, якщо крок ділиться, то K приймає значення одиниці.
а) Якщо , то навіть коли виконалась умова , крок не змінюється, тобто лишається тим самим (обчислення далі проводяться з попереднім кроком).
б) Якщо і виконалась умова , тоді .
В обох випадках а) і б) результат виводиться на друк.
Для розв’язку системи диференціальних рівнянь використовують цей самий метод, за виключенням того, що всі рівняння системи необхідно розв’язувати паралельно.
4.Лістинг програм
1. Метод Рунге – Кутта з автоматичною зміною кроку
using System;
using System.IO;
namespace RK
{
class Data
{
int Kr = 0;
double[,] Yp = new double[4, 1000000];
double[,] y = new double[4, 1000000];
double[] time = new double[1000000];
double[] t = new double[1000000];
double[] m = new double[4];
double[] n = new double[4];
double[] l = new double[4];
double[] p = new double[4];
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;
public Data()
{
Kp = 0.425;
Tm = 0.3;
Tk = 0.02;
Tl = 0.006;
C = 4;
i = 350;
k1 = 2;
k2 = 2;
S = 60;
Qin = 1;
dt = 0.0001;
e = 0.0001;
b = 5;
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] + 2 * m[1] + 2 * m[2] + m[3]) / 6.0;
y[1, k + 1] = y[1, k] + (n[0] + 2 * n[1] + 2 * n[2] + n[3]) / 6.0;
y[2, k + 1] = y[2, k] + (l[0] + 2 * l[1] + 2 * l[2] + l[3]) / 6.0;
y[3, k + 1] = y[3, k] + (p[0] + 2 * p[1] + 2 * p[2] + p[3]) / 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;
p[0] = f4(t, y1, y2, y3, y4) * h;
m[1] = F1(t + h / 2.0, y1 + m[0] / 2.0, y2 + n[0] / 2.0, y3 + l[0] / 2.0, y4 + p[0] / 2.0) * h;
n[1] = f2(t + h / 2.0, y1 + m[0] / 2.0, y2 + n[0] / 2.0, y3 + l[0] / 2.0, y4 + p[0] / 2.0) * h;
l[1] = f3(t + h / 2.0, y1 + m[0] / 2.0, y2 + n[0] / 2.0, y3 + l[0] / 2.0, y4 + p[0] / 2.0) * h;
p[1] = f4(t + h / 2.0, y1 + m[0] / 2.0, y2 + n[0] / 2.0, y3 + l[0] / 2.0, y4 + p[0] / 2.0) * h;
m[2] = F1(t + h / 2.0, y1 + m[1] / 2.0, y2 + n[1] / 2.0, y3 + l[1] / 2.0, y4 + p[1] / 2.0) * h;
n[2] = f2(t + h / 2.0, y1 + m[1] / 2.0, y2 + n[1] / 2.0, y3 + l[1] / 2.0, y4 + p[1] / 2.0) * h;
l[2] = f3(t + h / 2.0, y1 + m[1] / 2.0, y2 + n[1] / 2.0, y3 + l[1] / 2.0, y4 + p[1] / 2.0) * h;
p[2] = f4(t + h / 2.0, y1 + m[1] / 2.0, y2 + n[1] / 2.0, y3 + l[1] / 2.0, y4 + p[1] / 2.0) * h;
m[3] = F1(t + h, y1 + m[2], y2 + n[2], y3 + l[2], y4 + p[2]) * h;
n[3] = f2(t + h, y1 + m[2], y2 + n[2], y3 + l[2], y4 + p[2]) * h;
l[3] = f3(t + h, y1 + m[2], y2 + n[2], y3 + l[2], y4 + p[2]) * h;
p[3] = f4(t + h, y1 + m[2], y2 + n[2], y3 + l[2], y4 + p[2]) * 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 Provirka(double t, double y1, double y2, double y3, double y4, double h)
{
Zminna(t, y1, y2, y3, y4, dt);
y0 = y[1, k + 1];
Zminna(time[k] + dt, y1, y2, y3, y4, dt);
y2 = y[1, k + 1];
Zminna(time[k] + dt / 2.0, y1, y2, y3, y4, dt / 2.0);
y1 = y[1, k + 1];
Zminna(time[k] + dt, y1, y2, y3, y4, dt / 2.0);
y3 = y[1, k + 1];
Kr = -1;
while (Math.Abs(y3 - y1) > e)
{
dt = dt / 2.0;
Kr++;
y0 = y[1, k + 1];
Zminna(time[k] + dt / 2.0, y1, y2, y3, y4, dt / 2.0);
y2 = y[1, k + 1];
y3 = y0;
y1 = y2;
}
if (Kr != 0)
dt = 2.0 * dt;
y[1, k + 1] = y2;
}
public void Prod()
{
k = 1;
do
{
t[k] = dt;
Pohidni(time[k], y[0, k], y[1, k], y[2, k], y[3, k]);
Zminna(time[k], y[0, k], y[1, k], y[2, k], y[3, k], dt);
time[k + 1] = time[k] + dt;
k++;
} while (time[k] < b);
{
Random rand = new Random();
StreamWriter log_out;
log_out = new StreamWriter("logfile.txt");
Console.SetOut(log_out);
Console.WriteLine("t\ty1\ty1'");
for (int j = 0; j < k; j += rand.Next(1, 7))
Console.WriteLine("{0:0.######}\t{1:0.####}", time[j], y[1, j]);
log_out.Close();
}
public void vuv()
{
Console.WriteLine("t\ty1\ty1'");
for (int j = 0; j < k; j += 4)
{
Console.WriteLine("{0:0.######}\t{1:0.####}", time[j], y[1, j]);
}
}
}
class Program
{
static void Main(string[] args)
{
Data r = new Data();
r.Prod();
r.vuv();
r.vuvid();
Console.ReadLine();
}
}
}
2. Модифікований метод Ейлера
using System;
using System.IO;
class Data
{
long z = 1000000;
double[,] Yp;
double[,] y;
double[] time;
double[] m = new double[4];
double[] n = new double[4];
double[] l = new double[4];
double[] p = new double[4];
double dt;
double k1;
double KU;
double S;
double T1;
double k2;
double TK;
double C;
double TM;
double II;
long k;
int Qin;
double a;
double b;
double y1;
double y2;
double y3;
double y4;
double e;
double d;
public Data()
{
Yp = new double[4, z];
y = new double[4, z];
time = new double[z];
KU = 0.406;
TM = 0.3;
TK = 0.02;
T1 = 0.006;
C = 4;
II = 350;
k1 = 2;
k2 = 2;
S = 60;
Qin = 1;
dt = 0.001;
e = 0.0001;
a = 0;
b = 2.9;
time[0] = a;
}
public double F1(double t, double y1, double y2, double y3, double y4)
{
return (k1 * KU * S * (Qin - y2) - y1) / T1;
}
public double F2(double t, double y1, double y2, double y3, double y4)
{
return y4 / II;
}
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()
{
Yp[0, k] = F1(time[k], y[0, k], y[1, k], y[2, k], y[3, k]);
Yp[1, k] = F2(time[k], y[0, k], y[1, k], y[2, k], y[3, k]);
Yp[2, k] = F3(time[k], y[0, k], y[1, k], y[2, k], y[3, k]);
Yp[3, k] = F4(time[k], y[0, k], y[1, k], y[2, k], y[3, k]);
}
public void Znach()
{
y1 = y[0, k] + dt * F1(time[k], y[0, k], y[1, k], y[2, k], y[3, k]) / 2.0;
y2 = y[1, k] + dt * F2(time[k], y[0, k], y[1, k], y[2, k], y[3, k]) / 2.0;
y3 = y[2, k] + dt * F3(time[k], y[0, k], y[1, k], y[2, k], y[3, k]) / 2.0;
y4 = y[3, k] + dt * F4(time[k], y[0, k], y[1, k], y[2, k], y[3, k]) / 2.0;
y[0, k + 1] = y[0, k] + dt * F1(time[k] + dt / 2.0, y1, y2, y3, y4);
y[1, k + 1] = y[1, k] + dt * F2(time[k] + dt / 2.0, y1, y2, y3, y4);
y[2, k + 1] = y[2, k] + dt * F3(time[k] + dt / 2.0, y1, y2, y3, y4);
y[3, k + 1] = y[3, k] + dt * F4(time[k] + dt / 2.0, y1, y2, y3, y4);
{
Znach();
d = l[0] / 150.0 - l[2] * 3.0 / 100.0 + l[3] * 16.0 / 75.0 + l[4] / 20.0 - l[5] * 6.0 / 25.0;
if (Math.Abs(d) > e)
dt = dt / 2.0;
if (Math.Abs(d) < e / 30.0)
dt = 2.0 * dt;
} while (Math.Abs(d) > e || Math.Abs(d) < e / 30.0);
}
public void Prod()
{
do
{
Pohidni();
Znach();
time[k + 1] = time[k] + dt;
k++;
} while (time[k] < b);
Pohidni();
}
public void vuvid()
{
Random rand = new Random();
StreamWriter log_out;
log_out = new StreamWriter("logfile.txt");
Console.SetOut(log_out);
for (int j = 0; j < k; j += rand.Next(1, 5)) Console.WriteLine("{0:0.####}\t{1:0.###}", time[j], y[1, j]);
}
public void vuv()
{
Console.WriteLine("t\ty1\ty1'");
for (int j = 0; j < k; j += 2)
{
Console.WriteLine("{0:0.####}\t{1:0.###}", time[j], y[1, j]);
}
}
}
static void Main(string[] args)
{
Data d = new Data();
d.Prod();
d.vuv();
d.vuvid();
Console.ReadLine();
}
}
5.Результати виконання програми
1. Виконаної методом Рунге – Кутта з автоматичною зміною кроку t ω(t)
0 0
0 0
0,0001 0
0,0002 0
0,0004 0
0,0007 0
0,001 0
0,0016 0
0,0017 0
0,0019 0
0,0025 0
0,0028 0
0,003 0
0,0034 0
0,004 0
0,0044 0
0,0047 0
0,0049 0
0,005 0
0,0055 0
0,0059 0
0,006 0
0,0061 0
0,0067 0
0,0073 0
0,0076 0
0,0078 0
0,008 0
0,0085 0
0,0091 0
0,0097 0
0,0103 0
0,0106 0
0,0108 0
0,011 0
0,0116 0
0,0119 0
0,0121 0
0,0122 0
0,0128 0
0,013 0
0,0135 0
0,0136 0
0,0142 0
0,0147 0
0,0148 0
0,015 0
0,0154 0
0,0155 0
0,016 0
0,0162 0
0,0165 0,0001
0,0166 0,0001
0,0172 0,0001
0,0178 0,0001
0,0181 0,0001
0,0182 0,0001
0,0184 0,0001
0,0185 0,0001
0,0191 0,0001
0,0192 0,0001
0,0193 0,0001
0,0194 0,0001
0,0196 0,0001
0,0198 0,0001
0,0199 0,0001
0,02 0,0001
0,0206 0,0001
0,0208 0,0001
0,0214 0,0001
0,0215 0,0001
0,0217 0,0001
0,0218 0,0001
0,0219 0,0001
0,0222 0,0001
0,0228 0,0002
0,0234 0,0002
0,0235 0,0002
0,0237 0,0002
0,0243 0,0002
0,0244 0,0002
0,0249 0,0002
0,025 0,0002
0,0251 0,0002
0,0256 0,0002
0,0259 0,0002
0,0264 0,0002
0,0268 0,0003
0,027 0,0003
0,0275 0,0003
0,0281 0,0003
0,0283 0,0003
0,0284 0,0003
0,0286 0