МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ
НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ
„ЛЬВІВСЬКА ПОЛІТЕХНІКА”
ІКТА
Кафедра БІТ
Курсова робота
з курсу:
" Комп'ютерні методи дослідження інформаційних процесів і систем "
на тему:
" АВТОМАТИЧНИЙ ПОТЕНЦІОМЕТР З МАГНІТНИМ ПІДСИЛЮВАЧЕМ "
Тема 2, варіант 11
Львів- 2013
В даній роботі розглянено застосування методу Рунге-Кутта та Рунге-Кутта-Мерсона для дослідження перехідного процесу автоматичного потенціометра з диференціюючим контуром. Розв’язок поставленої задачі представлений в середовищі С# і платформі Visual Studio 2010. Графіки уточнень побудовані в середовищі С#.
Зміст
1. Постановка задачі……………………………………………………………….. 4
2. Перетворення рівнянь………………………………………………………….. 6
3. Теоретичні відомості……………………………………………………………. 8
3.1 Метод Рунге-Кутта для розв’язку систем диференціальних рівнянь….. 8
3.2 Метод Рунге-Кутта-Мерсона для розв’язку диференціальних рівнянь…………………………………………………………………………… 11
4. Лістинг програм………………………………………………………………… 12
4.1 Метод Рунге-Кутта з автоматичною зміною кроку………………………. 12
4.2 Метод Рунге-Кутта-Мерсона……….…………………………………………. 17
5. Результати виконання програми……………………………………………... 22
5.1 Метод Рунге-Кутта з автоматичною зміною кроку………………………… 22
5.2 Метод Рунге-Кутта-Мерсона……………...…………………………………… 26
6. Графік перехідного процесу…………………………………………………… 30
7. Список літератури………………………………………………………………. 31
1. Постановка задачі
Схема:
Рівняння ланок :
вимірювальна схема
електронний підсилювач
магнітний підсилювач
двигун
редуктор
При початкових параметрах
Параметри
11
(m (рад)
5
Un (мв)
50
Cu (г.см.в)
9
C( (г.см.сек/рад)
3
IД (г.см.сек2)
0,03
І (г.см.сек2)
7
КМ
10
Т (сек)
0,05
і
30
Звести систему алгебро-диференціальних рівнянь до системи трьох диференціальних рівнянь першого порядку, представити її у нормальній формі та розв’язати цю систему вказаними методами. Початкові умови - =1 радіан, решта початкових умов – нульові. Числові значення сталих параметрів, заданих в таблиці, слід зобразити з допомогою одиниць системи СІ.
Побудувати графік зміни величини
2. Перетворення рівнянь
МП Т + e2 = км e1
ЕП U = кп ·e2
2х фазний двигун змінного струму
I = CuU — C ω
Редуктор
Вимірювальна схема
е1 = (вх — вих)
Необхідно звести ці рівняння до системи ЗДР I-го порядку
T + e2 = км(вх — вих)
Розв ‘ язуємо відносно
= (вх — вих) —
Це перше рівняння системи
У рівняння (3) підставляємо (2) та (4), при цьому з рівняння (4) знаходимо
ω = = i =
I = CuU — C ω ;
I = Cu· кпe2 — C ω· ω
I=Iд+ Iн/ i2 = (1— ) —
=
=
Таблиця ідентифікаторів
I — I1
Ѳ — QM
Un — UN
T — T
Cu — CU KU – вибір експериментів
C ω — CW
Id — ID
In — IN
E2 —Y(1)
Ѳвих — Y(2)
ω — Y(3)
— F(1)
— F(2)
— F(3) II=ID+IN/(I1*I1)
Систему звичайних диф. рівнянь запишемо у вигляді підпрограми
F(1)=KM*UN/(T*QM)*(1-Y(2))-Y(1)/T
F(2)=Y(3)/I1
F(3)= (CU*KU*Y(1)-CW*Y(3))/
3. Теоретичні відомості
3.1 Метод Рунге – Кутта для розв’язку систем
диференціальних рівнянь.
Диференціальне рівняння (ДР), що містить лише одну незалежну змінну і похідні за нею, називають звичайними (ДР). ДР, що містить декілька незалежних змінних і похідні за ними, називають рівняння в частинних похідних.
Порядком ДР називається найвищий порядок похідної (або диференціалу), який входить в рівняння. Звичайне ДР (ЗДР) -го порядку в загальному випадку має незалежну змінну, невідому функцію та її похідні до -го порядку включно:
- незалежна змінна;
- невідома функція (залежна змінна);
- похідні цієї функції.
Диференціальне рівняння -го порядку, розв’язане відносно старшої похідної, може бути записано у вигляді:
Щоб розв’язати ЗДР, необхідно мати значення залежної змінної та (або) її похідних при деяких значення незалежної змінної.
В методі Рунге-Кутта значення функції визначається за формулою
Якщо розкласти функцію в ряд Тейлора і обмежитись членами до включно, то приріст можна записати у вигляді
(1)
Замість того, щоб обчислювати члени ряду за формулою (1) в методі Рунге-Кутта використовують наступні формули.
Похибка на кожному кроці має порядок . Таким чином метод Рунге-Кутта забезпечує високу точність, однак вимагає більшого об’єму обчислень.
Деколи зустрічається інша форма представлення методу Рунге-Кутта 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;
class Data
{
double i;
double Qm;
double Cu;
double Cw;
double Id;
double In;
double T;
double I;
double Kp;
double Qin;
double a;
double b;
double dt;
double Un;
double Km;
int k;
int Kr;
double[,] y = new double[3, 1000000];
double[,] Y = new double[3, 1000000];
double[] t = new double[1000000];
double[] p = new double[1000000];
double[] m = new double[5];
double[] l = new double[5];
double[] n = new double[5];
double z0;
double z1;
double z2;
double z3;
double e;
public Data()
{
dt = 0.0001;
e = 0.00001;
a = 0;
b = 1.5;
t[0] = a;
Qm = 5.0;
Un = 0.05;
Cu = 9.0 * Math.Pow(10.0, -5.0);
Cw = 3.0 * Math.Pow(10.0, -5.0);
Id = 0.03 * Math.Pow(10.0, -5.0);
In = 7.0 * Math.Pow(10.0, -5.0);
Km = 10.0;
T = 0.05;
i = 30.0;
Kp = 450.0;
Qin = 1.0;
I = Id + In / Math.Pow(i, 2);
}
public double F2(double t, double y1, double y2, double y3)
{
return Km * Un * (Qin - y1) / (Qm * T) - y2 / T;
}
public double F1(double t, double y1, double y2, double y3)
{
return y3 / i;
}
public double F3(double t, double y1, double y2, double y3)
{
return Cu * Kp * y2 / I - Cw * y3 / I;
}
public void Pohidni(double t, double y1, double y2, double y3)
{
Y[0, k] = F1(t, y1, y2, y3);
Y[1, k] = F2(t, y1, y2, y3);
Y[2, k] = F3(t, y1, y2, y3);
}
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;
}
public void Kof(double t, double y1, double y2, double y3, double h)
{
m[0] = F1(t, y1, y2, y3) * h;
n[0] = F2(t, y1, y2, y3) * h;
l[0] = F3(t, y1, y2, y3) * h;
m[1] = F1(t + h / 2.0, y1 + m[0] / 2.0, y2 + n[0] / 2.0, y3 + l[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) * h;
l[1] = F3(t + h / 2.0, y1 + m[0] / 2.0, y2 + n[0] / 2.0, y3 + l[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) * h;
n[2] = F2(t + h / 2.0, y1 + m[1] / 2.0, y2 + n[1] / 2.0, y3 + l[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) * h;
m[3] = F1(t + h, y1 + m[2], y2 + n[2], y3 + l[2]) * h;
n[3] = F2(t + h, y1 + m[2], y2 + n[2], y3 + l[2]) * h;
l[3] = F3(t + h, y1 + m[2], y2 + n[2], y3 + l[2]) * h;
}
public void Zminna(double t, double y1, double y2, double y3, double h)
{
Kof(t, y1, y2, y3, h);
Zmin();
}
public void Provirka(double t, double y1, double y2, double y3, double h)
{
Zminna(t, y1, y2, y3, dt);
z0 = y[0, k + 1];
Zminna(t + dt, y1, y2, y3, dt);
z1 = y[0, k + 1];
Zminna(t + dt / 2.0, y1, y2, y3, dt / 2.0);
z2 = y[0, k + 1];
Zminna(t + dt, y1, y2, y3, dt / 2.0);
z3 = y[0, k + 1];
Kr = -1;
A1: if ((Math.Abs(z3 - z1) > e))
{
z3 = z0;
z1 = z2;
dt = dt / 2.0;
Zminna(t, y1, y2, y3, dt);
z0 = y[0, k + 1];
Zminna(t + dt / 2.0, y1, y2, y3, dt / 2.0);
z2 = y[0, k + 1];
Kr++;
goto A1;
}
Zminna(t, y1, y2, y3, dt);
if (Kr == 0 || Kr == -1)
dt = 2.0 * dt;
}
public void Prod()
{
k = 1;
do
{
p[k] = dt;
Pohidni(t[k], y[0, k], y[1, k], y[2, k]);
Zminna(t[k], y[0, k], y[1, k], y[2, k], dt);
t[k + 1] = t[k] + dt;
k++;
} while (t[k] < b);
Pohidni(t[k], y[0, k], y[1, k], y[2, 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(10, 15))
{
Console.WriteLine("{0:0.######}\t\t{1:0.####}", t[j], y[0, 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[0, 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;
class Data
{
double i;
double Qm;
double Cu;
double Cw;
double Id;
double In;
double T;
double I;
double Kp;
double Qin;
double a;
double b;
double dt;
double Un;
double Km;
int k;
double[,] y = new double[3, 1000000];
double[,] Y = new double[3, 1000000];
double[] t = new double[1000000];
double[] p = new double[1000000];
double[] m = new double[5];
double[] l = new double[5];
double[] n = new double[5];
double z0;
double z1;
double z2;
double e;
public Data()
{
dt = 0.0001;
e = 0.00001;
a = 0;
b = 2;
t[0] = a;
Qm = 5.0;
Un = 0.05;
Cu = 9.0 * Math.Pow(10.0, -5.0);
Cw = 3.0 * Math.Pow(10.0, -5.0);
Id = 0.03 * Math.Pow(10.0, -5.0);
In = 7.0 * Math.Pow(10.0, -5.0);
Km = 10.0;
T = 0.05;
i = 30.0;
Kp = 450.0;
Qin = 1.0;
I = Id + In / Math.Pow(i, 2);
}
public double F2(double t, double y1, double y2, double y3)
{
return Km * Un * (Qin - y1) / (Qm * T) - y2 / T;
}
public double F1(double t, double y1, double y2, double y3)
{
return y3 / i;
}
public double F3(double t, double y1, double y2, double y3)
{
return Cu * Kp * y2 / I - Cw * y3 / I;
}
public void Pohidni(double t, double y1, double y2, double y3)
{
Y[0, k] = F1(t, y1, y2, y3);
Y[1, k] = F2(t, y1, y2, y3);
Y[2, k] = F3(t, y1, y2, y3);
}
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;
}
public void Kof(double t, double y1, double y2, double y3, double h)
{
m[0] = F1(t, y1, y2, y3) * h;
n[0] = F2(t, y1, y2, y3) * h;
l[0] = F3(t, y1, y2, y3) * h;
m[1] = F1(t + h / 3.0, y1 + m[0] / 3.0, y2 + n[0] / 3.0, y3 + l[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) * h;
l[1] = F3(t + h / 3.0, y1 + m[0] / 3.0, y2 + n[0] / 3.0, y3 + l[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) * 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) * 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) * 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) * 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) * 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) * 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]) * 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]) * 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]) * h;
}
public void Zminna(double t, double y1, double y2, double y3, double h)
{
Kof(t, y1, y2, y3, h);
Zmin();
}
public void Perevirka(double t, double y1, double y2, double y3, double h)
{
A1: Zminna(t, y1, y2, y3, 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, k]);
Zminna(t[k], y[0, k], y[1, k], y[2, k], dt);
t[k + 1] = t[k] + dt;
k++;
} while (t[k] < b);
Pohidni(t[k], y[0, k], y[1, k], y[2, 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, 13))
{
Console.WriteLine("{0:0.######}\t\t{1:0.####}", t[j], y[0, j]);
}
log_out.Close();
}
public void vuv()
{
Console.WriteLine("t\ty1\ty1'");
for (int j = 0; j < k; j += 30)
{
Console.WriteLine("{0:0.######}\t\t\t{1:0.#####}", t[j], y[0, j]);
}
}
}
class Program
{
static void Main(string[] args)
{
Data d = new Data();
d.Prod();
d.vuv();
d.vuvid();
Console.ReadLine();