Числові методи і математичне моделювання

Размещено на http://www.allbest.ru/

1. Побудова математичної моделі та числове її дослідження

1.1 Побудова математичної моделі

Конструктивні параметри: L1=100м, r1=0.05м, k_B=0.005м, d=1м.

Вхідні величини: P1=2000Па, P2=9500Па,T1=293K, T2=350K.

Керуючі величини: l2=0.4, l3=0.7.

Значення збурення: .

Параметри стану системи: h,T.

Згідно з рівнянням збереження маси речовини та ввівши деякі припущення (масообмін на границі розділу фаз рідина-повітря відсутній,=const), запишемо диференційне рівняння, що описує зніму рівня в ємності:

; (1.1)

де - площа ємності ,м;

h- рівень рідини в ємності, м;

Q- об'ємна витрата, м/с.

визначається наступною залежністю:

- для короткого турбулентного трубопроводу з регулюючим органом (P₃=0):

Об'ємна витрата Q1 визначається за формулою:

-для довгого ламінарного трубопроводу:

(1.2)

Об'ємна витрата Q2 визначається наступною залежністю:

- для короткого турбулентного трубопроводу з регулюючим органом :

Для отримання виразу зміни температури рідини в ємності, скористаємося законом збереження тепла. Запишемо:

(1.3)

або

(1.3')

Враховуючи рівняння (1.1) і (1.3') зробимо заміну і отримуємо:

(1.4)

Об'єднавши рівняння (1.1), (1.2), (1.4) отримаємо наступну систему звичайних диференціальних рівнянь:

(1.5)

Ця система звичайних диференційних рівнянь разом із початковими умовами є математичною моделлю даного об'єкту.

1.2 Числове дослідження математичної моделі

1. Знайдемо початкові умови.

Знайдемо розв'язок системи нелінійних рівнянь застосовуючи функцію fsolve . Для цього створимо файл даних yos_kr1.m :

%yos_kr1.m

function y=yos_kr1(x);

h=x(1); Q1=x(2); T=x(3);

%-------------------------------

P1=2000;P2=9500; T1=293; T2=350;

L1=100; r1=0.05; kv=0.005;

l2=0.4; l3=0.7; d=1; nu=0.00001;

dz=0.9; ro=1e3; g=9.81;

%-------------------------------

S=pi*d^2/4;

kl=pi*r1^4/8/L1/nu;

T=L1*kl/pi/r1^2;

Q2=kv*l2*sqrt((P2-ro*g*x(1))/ro);

Q3=kv*l3*sqrt(g*x(1));

y=[(x(2)+kv*l2*sqrt((P2-ro*g*x(1))/ro)-kv*l3*sqrt(g*x(1)))/S;

(kl*P1/ro-x(2))/T;

(x(2)*(T1-x(3))+kv*l2*sqrt((P2-ro*g*x(1))/ro)*(T2-x(3)))/S/x(1)];

Потім створюємо файл yos_kr11.m:

x0=[0.1 0.003 320];

y=fsolve('yos_kr1',x0)

Результатами виконання файлу є такі значення параметрів стану:

h0=0.617917411421; Q10=0.004908738521; T0=317.530396154693;

Знайдемо тепер номінальні значення параметрів стану об'єкту числовим методом, використовуючи функції MatLab. Для того, щоб розв'язати систему нелінійних диференціальних рівнянь (1.5) потрібно створити два файли, в одному з яких (файлі-функції yos_kr2.m) будуть записані праві частини системи диференціальних рівнянь, розв'язані відносно перших похідних, а у другому файлі (yos_kr21.m) буде записана функція MatLab ode45, призначена для розв'язування системи диференціальних рівнянь.

Файл-функції yos_kr2.m:

function y=yos_kr2(t,x);

h=x(1); Q1=x(2); T=x(3);

%----------------------------

P1=2400;P2=9500; T1=293; T2=350;

L1=100; r1=0.05; kv=0.005;

l2=0.4; l3=0.7; d=1; nu=0.00001;

dz=0.9; ro=1e3; g=9.81;

%----------------------------

S=pi*d^2/4;

kl=pi*r1^4/8/L1/nu;

T=L1*kl/pi/r1^2;

Q2=kv*l2*sqrt((P2-ro*g*x(1))/ro);

Q3=kv*l3*sqrt(g*x(1));

%--------------------------------------------------------

y=[(x(2)+kv*l2*sqrt((P2-ro*g*x(1))/ro)-kv*l3*sqrt(g*x(1)))/S;

(kl*P1/ro-x(2))/T;

(x(2)*(T1-x(3))+kv*l2*sqrt((P2-ro*g*x(1))/ro)*(T2-x(3)))/S/h];

Файл yos_kr21.m:

x0=[0.617917411421 0.004908738521 317.530396154693];

T=[0 500];

tol=odeset('abstol',[1e-20 1e-20 1e-20]);

[t,y]=ode45('yos_kr2',T,x0,tol) ;

plot(t,y(:,1));grid;

ylabel('h,m');

xlabel('t,sec');pause;

plot(t,y(:,2));grid;

ylabel('Q1,mkub/sec');

xlabel('t,sec');pause;

plot(t,y(:,3));grid;

ylabel('T,K');

xlabel('t,sec');

Результатом роботи цих програм є графіки перехідних процесів параметрів стану досліджуваного об'єкту.

Рисунок 1.1 Перехідний процес h(t) з номінального початкового стану рівноваги до нового стану рівноваги, при нанесенні на об'єкт збурення .



Рисунок 1.2 Перехідний процес Q1(t) з номінального початкового стану рівноваги до нового стану рівноваги, при нанесенні на об'єкт збурення .



Рисунок 1.3 Перехідний процес T(t) з номінального початкового стану рівноваги до нового стану рівноваги, при нанесенні на об'єкт збурення .

2. Дослідження системи шляхом лінеаризації

2.1 Лінеаризації системи відносно стану рівноваги

Запишемо систему (1.5) у вигляді

(2.1)

де

Лінеаризована система диференційних рівнянь матиме вигляд:

(2.2)

або

 (2.3)

вихідні величини : y₁(t)=h (t)

y₂(t)=T(t)

Часткові похідні правих частин системи нелінійних диференційних рівнянь (2.1) беремо по всіх параметрах стану, а також по тих вхідних величинах, відхилення яких від номінального режиму задане в завданні на курсову роботу:

Коефіцієнти матриці стану системи та вектора вхідних величин обчислюються наступним чином:

В матричній формі лінеаризована система диференціальних рівнянь матиме вигляд:

, (2.4)

де - матриця стану системи;

- вектор вхідних величин;

; .

- вектор параметрів стану системи;

- сигнал збурення.

Для обчислення коефіцієнтів власної у матриці системи та вектора вхідних величин складемо програму yos_kr3.m:

h0=0.617917411421; Q10=0.004908738521; T0=317.530396154693;

%--------------------------------------------------------

P1=2000;P2=9500; T1=293; T2=350;

L1=100; r1=0.05; kv=0.005;

l2=0.4; l3=0.7; d=1; nu=0.00001;

dz=0.9; ro=1e3; g=9.81;

S=pi*d^2/4;

kl=pi*r1^4/8/L1/nu;

T=L1*kl/pi/r1^2;

%----------------------------------

a11=-(kv*l2*g/2/sqrt((P2-ro*g*h0)/ro)+kv*l3*g/2/sqrt(g*h0))/S;

a12=1/S;

a21=0;

a22=-1/T;

a31=(Q10*(T1-T0)+kv*l2*(T2-T0)*sqrt((P2-ro*g*h0)/ro))/(-S)/h0^2+kv*l2*(T2--T0)*(-g)/2/sqrt((P2-ro*g*h0)/ro)/S/h0;

a32=(T1-T0)/S/h0;

a33=(-Q10-kv*l2*sqrt((P2-ro*g*h0)/ro))/S/h0;

b2=kl/T/ro;

a=[a11 a12 0 ;

0 a22 0 ;

a31 a32 a33 ]

b=[0 ;

b2;

0]

c=[1 0 0;

0 0 1]

d=[0;0]

Значення коефіцієнтів:

a = -0.0156 1.2732 0

0 -0.0320 0

-0.3540 -50.5457 -0.0178

b = 0

0.7854

0

c =1 0 0

0 0 1

d = 0

0

2.3 Побудова перехідних процесів в лінеаризованій системі та їх порівняння з відповідними перехідними процесами в нелінійній системі

математична модель середовище система

Для побудови перехідних процесів в лінеаризованій системі використаємо функцію MatLab STEP, призначену для знаходження реакції лінійної незбудженої системи на одиничне стрибкоподібне збурення. Результати виконання функції, користуючись властивістю однорідності лінійних систем, необхідно домножити на величину стрибка.

Для накладання графіків перехідних процесів параметрів стану нелінійної та лінеаризованої системи створимо файл yos_kr31.m :

P1=2000;Px=P1*0.2;

h0=0.617917411421 ; Q10=0.004908738521; T0=317.530396151744;

x0=[h0 Q10 T0];

T=[0:8:500];

tol=odeset('Reltol',3e-14);

[t1,y1]=ode45('yos_kr2',T,x0,tol);

yos_kr3;

c=[1 0 0;

0 0 1];

d=[0;0];

t=[0:500];

sys=ss(a,b,c,d);

[y2,t,x2]=step(sys,500);

x2=Px*x2;

plot(t1,y1(:,1),'o',t,x2(:,1)+0.617917411421);grid;

ylabel('h,m');

xlabel('t,sec');pause;

plot(t1,y1(:,2),'o',t,x2(:,2)+0.004908738521);grid;

ylabel('Q1,mkub/sec');

xlabel('t,sec');pause;

plot(t1,y1(:,3),'o',t,x2(:,3)+317.530396154693);grid;

ylabel('T,K');

xlabel('t,sec');

Зауваження: Для порівняння графіків відхилень , та дійсних величин ,, до величин відхилень потрібно додати номінальні значення параметрів h_H=0.617917411421, Q_2H=0.004908738521, T_H=317.530396154693.

Результатом виконання цих програм є наступні графіки перехідних процесів:

Рисунок 2.1 Графік порівняння перехідних процесів з номінального початкового стану рівноваги до нового стану рівноваги по рівню h(t) для нелінійної (-) та лінеаризованої (?) моделей.при нанесенні на об'єкт збурення



Рисунок 2.2 Графік порівняння перехідних процесів з номінального початкового стану рівноваги до нового стану рівноваги по рівню Q1(t) для нелінійної (-) та лінеаризованої (?) моделей.при нанесенні на об'єкт збурення

Рисунок 2.3 Графік порівняння перехідних процесів з номінального початкового стану рівноваги до нового стану рівноваги по рівню T(t) для нелінійної (-) та лінеаризованої (?) моделей.при нанесенні на об'єкт збурення

3. Класичні методи дослідження систем

3.1 Зведення лінеаризованої системи звичайних диференційних рівнянь до одного звичайного диференціального рівняння вищого порядку відносно:

а) рівня h в гідравлічній ємності:

Для зведення лінеаризованої системи до одного рівняння відносно рівня h, запишемо її в операторній формі

(3.1)

Згідно правила Крамера , (3.2)

де ,

, s - оператор диференціювання.

Остаточно після обчислення визначників з врахуванням того, що , , спрощення та групування констант при змінних співвідношеннях (3.2) прийме вигляд :

. (3.3)

Отже, лінійну систему звичайних диференціальних рівнянь запишемо як диференціальне рівняння вищого порядку відносно так:

 (3.4)

де

Маючи залежність параметрів стану від вхідної величини, можна записати функцію передачі системи, яка записується як відношення вхідного оператора до власного оператора системи:

Функцію передачі системи, зведеної відносно рівня в ємності, запишемо у вигляді:

. (3.5)

б) температури в гідравлічній ємності Т:

Аналогічно, як в п.3.1. приведемо систему 2.3 до вигляду для застосування правила Крамера, яке в даному випадку буде мати вигляд:

, (3.6)

де ,

Остаточно після обчислення визначників з врахуванням того, що , , спрощення та групування констант при змінних співвідношеннях (3.6) прийме вигляд:

(3.7)

Отже, лінійну систему звичайних диференціальних рівнянь запишемо як диференціальне рівняння вищого порядку відносно так:

(3.8)

де

Маючи залежність параметрів стану від вхідної величини, можна записати функцію передачі системи, яка записується як відношення вхідного оператора до власного оператора системи:

Функцію передачі системи, зведеної відносно рівня в ємності, запишемо у вигляді:

. (3.9)

3.2 Одержання аналітичних виразів перехідних та імпульсних перехідних функцій систем, отриманих в п.3.1.

а). Знайдемо імпульсну перехідну функцію системи (3.4):

 (3.10)

Початкові умови , знайдемо із формул :

Для знаходження розв`язку рівняння (3.10), знайдемо корені його характеристичного рівняння. Для цього створимо скрипт-файл (yos_kr4.m):

yos_kr3;

A0=a11*a22;B0=a12*b2;

A1=-a11-a22;A2=1;

h0=0; h0p=B0;

lb=roots([A2 A1 A0])

Результати виконання програми:

lb =

-0.03200000000000

-0.01561419577680

Оскільки корені характеристичного рівняння системи (3.10) дійсні та різні, то розв'язок має вигляд: , (3.11)

де - константи, які знайдемо з початкових умов.

(3.12)

Для знаходження констант використаємо функцію LINSOLVE з пакету MATLAB. Створимо файл yos_kr41.m:

A=[1 1;

lb(1) lb(2)];

H=[h0;h0p];

C=linsolve(A,H);

C=double(C)

Одержимо константи:

C =1.0e-005 *

-0.61028435734888

0.61028435734888

Для знаходження імпульсної перехідної функції запишемо рівняння:

.Підставимо h(t):

. (3.13)

б) Знайдемо імпульсну перехідну функцію системи (3.8):

(3.14)

Початкові умови , змінені за рахунок дії одиничного імпульсного вхідного сигналу , знайдемо застосовуючи відомі залежності:



Для знаходження розв`язку рівняння (3.14), знайдемо корені його характеристичного рівняння. Для цього створимо скрипт-файл (yos_kr5.m):

kr3;

A3=1; A2=-a11-a22-a33;

A1=a11*a33+a22*a33+a11*a22;

A0=-a11*a22*a33;B2=0;

B1=a32*b2; B0=a12*a31*b2-a32*a11*b2;

Т0=0;Т0p=B1;Т0p2=B0-A2*B1;

lb=roots([A3 A2 A1 A0])

Результати виконаня:

lb = -0.03200000000000

-0.01775609460462

-0.01561419577680

Розв'язок системи має вигляд:

, (3.15)

де ,C₃- константи, які знайдемо з початкових умов.

 (3.16)

Для знаходження констант створимо файл yos_kr51.m:

A=[1 1 1;

lb(1) lb(2) lb(3);

lb(1)^2 lb(2)^2 lb(3)^2];

H=[Т0;Т0p;Т0p2];

C=linsolve(A,H);

C=double(C);

Результати:

C = 0.00012704856739

0.00088148763791

-0.00100853620530

Для знаходження імпульсної перехідної функції запишемо рівняння:

.Підставимо h(t):

(3.17)

3.3 Порівняння графіків перехідних функцій, отриманих за аналітичним розв'язком та числовим методом.

а) по рівню h

Створимо файл yos_kr6.m:

t=[0:3:500];

yos_kr3;

yos_kr4;yos_kr41;

h=C(1)*exp(lb(1)*t)+C(2)*exp(lb(2)*t);

sys=ss(a,b,c,d);

[y,tt,x]=impulse(sys,500),

plot(t,h,'.k',tt,x(:,1));grid;

title('impulsna perexidna f');

ylabel('h-h0,m');

xlabel('time,sec');

pause;

r=C(1)/lb(1)*(exp(lb(1)*t)-1)+C(2)/lb(2)*(exp(lb(2)*t)-1);

sys=ss(a,b,c,d);

[y,tt,x]=step(sys,500),

plot(t,r,'.k',tt,x(:,1));

grid;

title('perehidna f');

ylabel('h-h0,m');

xlabel('time,sec');

В результаті отримуємо:

Рисунок 3.1. Графіки порівняння імпульсної перехідної функції,одержаної:

«-»- за аналітичною залежністю; «?» - із застосуванням функції IMPULSE.



Рисунок 3.2. Графіки порівняння перехідної функції,одержаної: «-»- за аналітичною залежністю; «?» - із застосуванням функції STEP

б)по температурі Т:

Створимо файл yos_kr61.m:

t=[0:8:500];

yos_kr3;

yos_kr5;yos_kr51;

h=C(1)*exp(lb(1)*t)+C(2)*exp(lb(2)*t)+C(3)*exp(lb(3)*t);

sys=ss(a,b,c,d);

[y,tt,x]=impulse(sys,500),

plot(t,h,tt,x(:,3),'.k');grid;

title('impulsna perexidna f');

ylabel('T-T0,K');

xlabel('time,sec');pause;

r=C(1)/lb(1)*(exp(lb(1)*t)-1)+C(2)/lb(2)*(exp(lb(2)*t)-1)+C(3)/lb(3)*(exp(lb(3)*t)-1);

sys=ss(a,b,c,d);

[y,tt,x]=step(sys,500),

plot(t,r,tt,x(:,3),'.k');grid;

title('perehidna f');

ylabel('T-T0,K');

xlabel('time,sec');

В результаті отримуємо:

Рисунок.3.3. Графіки порівняння імпульсної перехідної функції,одержаної:

«-»- за аналітичною залежністю; «?» - із застосуванням функції IMPULSE.



Рисунок 3.4 Графіки порівняння перехідної функції,одержаної: «-»- за аналітичною залежністю; «?» - із застосуванням функції STEP.

3.4 Аналітичні вирази для визначення реакції системи на вхідний синусоїдальний сигнал () за допомогою інтеграла згортки

Запишемо інтеграл згортки:

а) по рівню h.

Для порівняння графіків за допомогою аналітичних виразів створимо файл yos_kr7.m:

yos_kr3;

c=[1 0 0]; d=[0];

sys = ss(a, b, c, d);

w=0.2;

t=[0:1:500];

U=sin(w*t);

[y,t] = lsim(sys, U,t);

yos_kr4;yos_kr41;

s=C(1)*(w*exp(lb(1)*t)-(w*cos(w*t)+lb(1)*sin(w*t)))/(w^2+lb(1)^2)+C(2)*(w*exp(lb(2)*t)-(w*cos(w*t)+lb(2)*sin(w*t)))/(w^2+lb(2)^2);

figure(1)

plot(t,s,t,y,'.k');grid;

ylabel('h,m');

xlabel('time,sec');

yos_kr3;

w=0.2;

t=[0:1:500];

yos_kr5;yos_kr51;

s=C(1)*(w*exp(lb(1)*t)-(w*cos(w*t)+lb(1)*sin(w*t)))/(w^2+lb(1)^2)+C(2)*(w*exp(lb(2)*t)-(w*cos(w*t)+lb(2)*sin(w*t)))/(w^2+lb(2)^2)+C(3)*(w*exp(lb(3)*t)-(w*cos(w*t)+lb(3)*sin(w*t)))/(w^2+lb(3)^2);

c=[0 0 1]; d=[0];

sys = ss(a, b, c, d);

y = lsim(sys, sin(w*t),t);

figure(2)

plot(t,s,t,y,'.k');grid;

ylabel('T,K');

xlabel('time,sec');

Отримуємо графіки:

Рисунок 3.5. Графіки порівняння перехідних процесів по рівню h при дії вхідного сигналу u=sin?t, одержані: «-»- за аналітичною залежністю; «?» - із застосуванням функції LSIM.



Рисунок 3.6.Графіки порівняння перехідних процесів по температурі Т при дії вхідного сигналу u=sin?t, одержані: «-»- за аналітичною залежністю; «?» - із застосуванням функції LSIM.

4. Частотні методи аналізу системи

Амплітудочастотна характеристика (А(w)) - це залежність відношення амплітуди вихідного періодичного сигналу до вхідної амплітуди цього сигналу як функції частоти в усталеному режимі роботи.

Фазочастотна характеристика ((w)) - це залежність зміни фази вихідного сигналу по відношенню до фази вхідного сигналу від частоти в усталеному режимі роботи.

а) по рівню h

В усталеному режимі коливань складова прямує до нуля. Отже, при ,

Звідси визначаємо АЧХ і ФЧХ:

Запишимо рівняння (3.5) замінивши ,отримаємо:

(4.1)

З (4.1) виділимо дійсну та уявну частини:

Для порівняння графіків АЧФ і ФЧХ по рівню h створимо файл yos_kr8.m:

yos_kr3;

w=[0:0.01:0.5];

yos_kr4;yos_kr41;

p=j*w;

Hh=B0./(A2.*p.^2+A1.*p+A0);

Aw=sqrt(real(Hh).^2+imag(Hh).^2);

[Avuh,fivuh]=bode([B0],[A2 A1 A0],w);

plot(w,Aw,w,Avuh,'.k');grid;

ylabel('Аvuh/Аvh');

xlabel('w,rad/s');

pause;

Fiw=atan(imag(Hh)./real(Hh));

Fiw=Fiw*180/pi;

plot(w,Fiw,w,fivuh,'.k');grid;

ylabel('fivuh/fivh');

xlabel('w,rad/s');

Одержуємо такі графіки:

Рисунок 4.3 Порівняння графіків амплітудо-частотної характеристики системи : «-» - за аналітичною залежністю; «?» - із застосуванням функції MatLab BODE

Рисунок 4.4 Порівняння графіків фазо-частотної характеристики системи: «-» - за аналітичною залежністю; «?» - із застосуванням функції MatLab BODE.

б) по температурі Т:

В усталеному режимі коливань складова прямує до нуля. Отже, при ,

Звідси визначаємо АЧХ і ФЧХ:

Запишимо рівняння (3.5) замінивши ,отримаємо:

(4.2)

З (4.2) виділимо дійсну та уявну частини:

Для порівняння графіків АЧФ і ФЧХ по рівню T створимо файл yos_kr81.m:

yos_kr3;

w=[0:0.007:0.5];

yos_kr5;yos_kr51;

p=j*w;

Hh=B1.*p+B0./(A3.*p.^3+A2.*p.^2+A1.*p+A0);

Aw=sqrt(real(Hh).^2+imag(Hh).^2);

[Avuh,fivuh]=bode([B1 B0],[A3 A2 A1 A0],w);

plot(w,Aw,w,Avuh,'.k');grid;

ylabel('Аvuh/Аvh');

xlabel('w,rad/s');

pause;

Fiw=atan(imag(Hh)./real(Hh));

Fiw=Fiw*180/pi;

plot(w,Fiw,w,fivuh,'.k');grid;ylabel('fivuh/fivh');xlabel('w,rad/s');

Одержуємо такі графіки:

Рисунок. 4.5 Порівняння графіків амплітудо-частотної характеристики системи: «-» - за аналітичною залежністю; «?» - із застосуванням функції MatLab BODE



Рисунок 4.6 Порівняння графіків фазо-частотної характеристики системи :

«-» - за аналітичною залежністю; «?» - із застосуванням функції MatLab BODE.

5. Дослідження моделі в середовищі SimuLink

Дослідимо отриману математичну модель в середовищі SimuLink. Побудуємо, використовуючи блоки SimuLink, структурну модель двоємнісного обєкту.

Для запуску моделі в SimuLink попередньо потрібно ініціалізувати початкові значення параметрів стану та констант з допомогою програми yos_kr9.m

P1=2000;P2=9500;

T1=293; T2=350;

L1=100; r1=0.05;

kv=0.005;

l2=0.4; l3=0.7;

d=1; nu=0.00001;

dz=0.9; ro=1e3; g=9.81;

S=pi*d^2/4;

kl=pi*r1^4/8/L1/nu;

A=L1*kl/pi/r1^2;

h0=0.617917411421; Q10=0.004908738521; T0=317.530396154693;

Рисунок 5.1 Модель системи у вікні SIMULINK

В даній моделі встановимо для блоків наступні параметри:

MATLAB Fcn: kv*l2*sqrt((P2-ro*g*u)/ro)-kv*l3*sqrt(g*u);

MATLAB Fcn1: kl*u/ro;

MATLAB Fcn2: (u(1)*(T1-u(3))+kv*l2*sqrt((P2-ro*g*u(2))/ro)*(T2-u(3)))/u(2);

Gain: 1/S;

Gain1: 1/A;

Gain1: 1/S;

Встановимо початкові умови для інтеграторів:

Integrator: h0;

Integrator1: Q10;

Integrator2: T0;

Параметри блоку STEP:

Step time: 0;

Initial value: 2000;

Final value: 2400;

Параметри симуляції:

Вкладка Solver:

Start time: 0;

Stop time: 500;

Save to workspace:

Time: tout;

States: xout;

Порівняння одержаних перехідних процесів в середовищах Matlab та Simulink.

Для накладання графіків перехідних процесів, отриманих в середовищі SimuLink і за допомогою розв'язування системи рівнянь функцією ode45 потрібно запустити на виконання файл yos_kr10.m :

x0=[0.617917411421 0.004908738521 317.530396154693];

T=[0 500];

tol=odeset('RelTol',3e-14);

[t,y]=ode45('yos_kr2',T,x0,tol);

plot(t,y(:,1),tout,xout(:,1),'.k');grid;

ylabel('h,m');

xlabel('t,sec');pause;

plot(t,y(:,2),tout,xout(:,2),'.k');grid;

ylabel('Q1,mkub/sec');

xlabel('t,sec');pause;

plot(t,y(:,3),tout,xout(:,3),'.k');grid;

ylabel('T,K');

xlabel('t,sec');

Результатом виконання програми є наступні графіки:

Рисунок 5.1 Графіки порівняння перехідних процесів для h, отримані:

«-» - шляхом розв'язку системи диференціальних рівнянь за допомогою функції ode45; «?» - шляхом побудови схеми моделі у вікні SIMULINK.



Рисунок 5.2 Графіки порівняння перехідних процесів для Q1, отримані: «-» - шляхом розв'язку системи диференціальних рівнянь за допомогою функції ode45; «?» - шляхом побудови схеми моделі у вікні SIMULINK

Рисунок 5.3 Графіки порівняння перехідних процесів для T, отримані: «-» - шляхом розв'язку системи диференціальних рівнянь за допомогою функції ode45; «?» - шляхом побудови схеми моделі у вікні SIMULINK.

Размещено на Allbest.ru