Phụ lục
MỘT SỐ PHƯƠNG PHÁP GIẢI CÁC
PHƯƠNG TRÌNH VI PHÂN TRÊN MÁY TÍNH ĐIỆN TỬ
Ngày nay có rất nhiều phần mềm ứng dụng cho phép nghiên cứu thuận lợi các
hệ thống động học, thí dụ Control system toolbox, Simulink. Tuy nhiên, việc áp dụng
chúng trong một số trường hợp cũng gặp khó khăn. Control system toolbox chỉ được
áp dụng cho các hệ tuyến tính. Trong Simulink, nhiều khâu phi tuyến trong thực tế
chưa có sẵn trong thư viện. Khi sử dụng Simulink để mô phỏng các hệ thống tự dẫn
KCB, ta không biết chính xác thời gian tự dẫn, vì vậy, không thể đặt thời gian mô
phỏng chính xác. Ta sẽ gặp nhiều khó khăn khi tối ưu tham số của hệ thống điều khiển
tự động trong Simulink. Vì vậy, trong nhiều trường hợp, việc chúng ta tự giải hệ
phương trình vi phân mô tả động học của hệ thống là thuận lợi hơn cả.
Để giải các phương trình vi phân trên máy tính điện tử, ta cần viết chúng dưới
dạng hệ các phương trình vi phân bậc nhất (dạng Kôsi), thí dụ
)(
,,
tx
yfy
=

. Sau
đó ta có thể sử dụng một trong các phương pháp sau để giải chúng trên máy tính điện
tử.
1. PHƯƠNG PHÁP ƠLE
Đây là một trong các phương pháp đơn giản nhất, theo đó lời giải của phương
trình
)(
,,
tx
yfy
=

Thí dụ 1: giả sử hàm truyền của phần tuyến tính của hệ truyền động có dạng
1
)()(
)(
)(
1
+
==
p
T
p
K
p
x
py
p
K
.
Hàm truyền trên tương ứng với hệ hai phương trình bậc nhất







+−=
=
,
1

t
y

=
;
x(t)-toạ độ vào.
Khi sử dụng phương pháp Ơle, ta nhận được lời giải của hệ hai phương trình
trên ở bước tích phân thứ i như sau
y
1
(i) = y
1
(i-1) + T
k
y
2
(i-1);
y
2
(i) = (1-T
k
/T) y
2
(i-1) + (T
k
/T)Kx(i-1).
Khi sử dụng phương pháp Tastin, ta nhận được lời giải của hệ hai phương trình
trên ở bước tích phân thứ i như sau
y
1

,,,
;,,
)(
)(
2122
2111
t
yyfy
t
yyfy


tại bước tích phân thứ i được xác định như sau





++++−=
++++−=
,
622
1
;
622
1
/)()()(
/)()()(
42322212
22

(i-1), y
2
(i-1));
k
21
= T
k
f
1
(y
1
(i-1) + k
11
/2, y
2
(i-1) + k
12
/2);
k
22
= T
k
f
2
(y
1
(i-1) + k
11
/2, y
2

(i-1) + k
22
/2);
2
k
41
= T
k
f
1
(y
1
(i-1) + k
31
, y
2
(i-1) + k
32
);
k
42
= T
k
f
1
(y
1
(i-1) + k
31
, y

ixk
k
)(
1
11
−=
;
)/)((
21
2
21
TT
ixk
kk
+−=
;
)//)((
421
32
31
TTT
ixk
kkk
++−=
;
)//)((
421
432
41
TTTT

Thông thường, khi tích phân trên máy tính điện tử ta thường sử dụng vòng lặp
while-end; hoặc if-end. Khi đó, bằng việc sơ bộ tính biểu thức
)///(
24
6
2
432
TTTT
kkkk
+++
(đối với thí dụ trên, còn với các phương trình khác ta sẽ có
biểu thức dưới dạng khác) bên ngoài vòng lặp while-end; hoặc if-end ta sẽ tiết kiệm
được thời gian tính toán so với hai phương pháp Ơle và Tastin. Chương trình tích phân
thí dụ trên được viết trong môi trường Matlab với x(t)=1 có dạng
clear all;
clc;
Tk=0.01;
a1=Tk;
a2=Tk^2/2;
a3=Tk^3/6;
a4=Tk^4/24;
3
x(1)=1;
y(1)=0;
t(1)=0;
n=2;
while n<=5002
n1=n-1;
x(n)=1;
y(n)=y(n1)+x(n1)*(a1+a2+a3+a4);

Đồ thị đặc trưng quá độ của khâu tích phân trên được thể hiện trên H.pl-1.
Trong tài liệu [2] các hệ số k
41
và k
42
có dạng
k
41
= T
k
f
1
(y
1
(i-1) + k
31
/2, y
2
(i-1) + k
32
/2);
k
42
= T
k
f
1
(y
1
(i-1) + k

ixi
y
i
y
ix
kkkk
(2)
5
Phương pháp Tastin
Hình PL-1. Đặc trưng quá độ của khâu tích phân
t, s
y(t)
Phương pháp Runge - kutta
Biểu thức của y(i) trong (1) và (2) chỉ khác nhau không đáng kể các hệ số của
các số hạng bậc cao của bước tích phân. Với bước tích phân nhỏ (nhỏ hơn rất nhiều
so với 1) thì kết quả của hai cách tính không khác nhau đáng kể.
Thí dụ 3: Dựng đặc trưng quá độ của khâu khuyếch đại quán tính
1
)(
)(
)(
+
==
Tp
K
p
x
py
p
K

a
T
a
T
ixbi
y
ai
y
i
y
kkkk
−+−−+−−+−=
.
Chương trình được viết trong môi trường Matlab có dạng như sau
clear all;
clc;
T=0.1;
K=1;
a=1/T;
b=K/T;
Tk=0.01;
a1=Tk;
a2=-a*Tk^2/2;
a3=a^2*Tk^3/6;
a4=-a^3*Tk^4/24;
x(1)=0;
y(1)=0;
t(1)=0;
n=2;
while n<=1000

[T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)
hoặc
sol = solver(odefun,[t0 tf],y0 ),
trong đó, solver là một trong các thuật giải ode45, ode23, ode113, ode15s,
ode23s, ode23t, or ode23tb.
[T, Y] là các ma trận cột thời gian và các hàm cần tìm;
odefun là tên file chương trình con mô tả các phương trình vi phân; thí dụ, tên
file chương trình con là “rigid” sẽ được viết như sau 'rigid'
tspan là véc tơ chỉ khoảng thời gian tích phân; thí dụ, khoảng thời gian tích
phân từ 0 tới 10 s được thể hiện như sau [0 10]
7
y0 là véc tơ điều kiện ban đầu của các hàm cần tìm;
options cho phép thay đổi các tham số tích phân.
Thí dụ 4. Giải hệ phương trình vi phân







−=
−=
=
,51,0
;
;
213
312
321

b=K/T;
dy = zeros(2,1); % a column vector
dy(1)=-a*y(1)+b*y(2);
dy(2) = 0;
Cú pháp của thuật giải có dạng
[T,Y] = ode45('rigid1',[0 10],[0 1]);
8
plot(T,Y(:,1),'-',T,Y(:,2),' '); hold on; grid on;
9