Nội dung text: 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ử
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ụ y f ( y , x , t ) . 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 y f ( y , x , t ) tại bước tích phân thứ i được xác định như sau
y(i) = y(i-1) + f[y(i-1), x(i-1)]Tk,
trong đó Tk-bước tích phân.
Để tính hàm giá trị của y tại bước tích phân đầu tiên ta sử dụng điều kiện ban
đầu y(t0)=y0. Phương pháp này có độ chính xác lớn nhất khi được sử dụng để để dựng
đặc tính quá độ của khâu tích phân. Phương pháp này có độ chính xác thấp hơn so với
các phương pháp khác khi được sử dụng để mô phỏng các khâu động học khác có
quán tính.
2. PHƯƠNG PHÁP TASTIN HAY PHƯƠNG PHÁP “HÌNH THANG”
Đây là phương pháp có độ chính xác rất cao. Theo phương pháp Tastin, lời giải
của phương trình y f ( y , x , t ) tại bước tích phân thứ i được xác định như sau
y(i) = y(i-1) + (Tk/2){3f[y(i-1), x(i-1)]- f[y(i-2), x(i-2)]}.
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
y (p ) K
K (p ) .
x ( p ) p (T p 1 )
Hàm truyền trên tương ứng với hệ hai phương trình bậc nhất
y 1 (t ) y 2 (t ) ;
1 K
2
y (t ) y (t ) x (t ) ,
T 2 T
trong đó các biến số y1(t), y2(t) có nghĩa như sau: y1(t)=y(t)-toạ độ ra; y 2 (t ) y 1 (t ) ;
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
y1(i) = y1(i-1) + Tky2(i-1);
y2(i) = (1-Tk/T) y2(i-1) + (Tk/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
y1(i) = y1(i-1) + (Tk/2)[3y2(i-1) - y2(i-2)];
y2(i) = y2(i-1) + (Tk/2T)[y2(i-2) + 3Kx(i-1) - 3y2(i-1) - Kx(i-2)].
3. PHƯƠNG PHÁP RUNGE - KUTTA
Đây là phương pháp có độ chính xác rất cao. Theo đó lời giải của hệ phương
trình
y 1 f 1 ( y1 , y 2 , t ) ;
y 2 f 2 ( y1 , y 2 , t ) ,
tại bước tích phân thứ i được xác định như sau
y1 (i ) y1 (i 1 ) (k11 2k 21 2k 31 k 41) / 6 ;
y 2 (i ) y 2 (i 1 ) (k12 2k 22 2k 32 k 42) / 6 ,
trong đó
k11 = Tkf1(y1(i-1), y2(i-1));
k12 = Tkf2(y1(i-1), y2(i-1));
k21 = Tkf1(y1(i-1) + k11/2, y2(i-1) + k12/2);
k22 = Tkf2(y1(i-1) + k11/2, y2(i-1) + k12/2);
k31 = Tkf1(y1(i-1) + k21/2, y2(i-1) + k22/2);
k32 = Tkf2(y1(i-1) + k21/2, y2(i-1) + k22/2);
2
k41 = Tkf1(y1(i-1) + k31, y2(i-1) + k32);
k42 = Tkf1(y1(i-1) + k31, y2(i-1) + k32).
Thí dụ 2: Dựng đặc trưng quá độ của khâu tích phân
y (p ) 1
K (p ) ( ) .
x p p
Hàm truyền trên tương ứng với phương trình vi phân
y (t ) x (t ) .
Để dựng đặc trưng quá độ của khâu tích phân cần bổ xung thêm phương trình
x (t ) 0 (vì x(t) = 1).
Khi sử dụng phương pháp Runge - kutta để tích phân phương trình trên ta có
k11 x (i 1 ) T k ;
2
k 21 x (i 1 )(T k T k / 2 ) ;
2 3
k 31 x (i 1 )(T k T k / 2 T k / 4 ) ;
2 3 4
k 41 x (i 1 )(T k T k T k / 2 T k / 4 ) ;
Vì vậy
x (i ) 1 ;
2 3 4 (1)
y (i ) y (i 1 ) x (i 1 )(T k T k / 2 T k / 6 T k / 24 ) .
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
2 3 4
(T k T k / 2 T k / 6 T k / 24 ) (đố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);
t(n)=n1*Tk;
n=n+1;
end
plot(t(1:n1),x(1:n1),'k'); hold on;
plot(t(1:n1),y(1:n1),'k'); grid on;
Đồ thị đặc trưng quá độ của khâu tích phân trên được thể hiện trên H.pl-1.
Khi áp dụng phương pháp Tastin thì chương trình dựng đặc trưng quá độ của
khâu tích phân trên trong môi trường Matlab có dạng
clear all;
clc;
Tk=0.01;
Tk1=Tk/2;
x(1)=1;
x(2)=1;
y(1)=0;
y(2)=0;
t(1)=0;
t(2)=0;
n=3;
while n<=5003
n1=n-1;
n2=n-2;
x(n)=1;
y(n)=y(n1)+Tk1*(3*x(n1)-x(n2));
t(n)=n2*Tk;
n=n+1;
end
4
plot(t(2:n2),x(2:n2),'k'); hold on;
plot(t(1:n1),y(1:n1),'k--'); grid on;
Đồ thị đặc trưng quá độ của khâu tích phân trên được thể hiện trên H.pl-1.
Phương pháp Runge - kutta
y(t)
Phương pháp Tastin
t, s
Hình PL-1. Đặc trưng quá độ của khâu tích phân
Trong tài liệu [2] các hệ số k41 và k42 có dạng
k41 = Tkf1(y1(i-1) + k31/2, y2(i-1) + k32/2);
k42 = Tkf1(y1(i-1) + k31/2, y2(i-1) + k32/2),
theo nghiên cứu của tác giả, nó không làm thay đổi đáng kể kết quả tích phân. Nếu ta
sử dụng giá trị k41 và k42 này thì kết quả tích phân bước thứ i của thí dụ 2 có dạng
x (i ) 1 ;
2 3 4 (2)
y (i ) y (i 1 ) x (i 1 )(T k 5 T k /12 T k / 8 T k / 48 ) .
5
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
y (p ) K
K (p ) .
x ( p ) Tp 1
Hàm truyền trên tương ứng với phương trình vi phân
1 K
y (t ) y (t ) x (t ) a y (t ) b x (t ) .
T T
Giá trị của hàm y(t) tại bước tích phân thứ i được tính như sau
y (i ) y (i 1 ) [ a y (i 1 ) b x (i 1 )](T k a T 2k / 2 a2 T 3k / 6 a3 T 4k / 24 ) .
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
n1=n-1;
x(n)=1;
y(n)=y(n1)+(-a*y(n1)+b*x(n1))*(a1+a2+a3+a4);
t(n)=n1*Tk;
n=n+1;
6
end
plot(t(1:n1),x(1:n1),'b'); hold on;
plot(t(1:n1),y(1:n1),'r'); grid on;
Khi sử dụng phương pháp Runge-Kutta để giải các phương trình vi phân, thời
gian tính toán hết ít hơn so với phương pháp Tastin (với cùng bước tích phân). Tuy
nhiên khi thời gian tích phân càng dài thì độ chính xác của phương pháp Tastin càng
cao hơn (bạn đọc có thể dễ dàng kiểm tra, thí dụ H.PL-1). Tuy nhiên, đây là hai
phương pháp tích phân có độ chính xác rất cao. Trong thí dụ 2, tại thời điểm 45 s sai
số của phương pháp Runge-Kutta vào khoảng 0,45 %, còn với phương pháp Tastin sai
số đạt tới ≈ 0 %.
4. MỘT SỐ PHƯƠNG PHÁP GIẢI HỆ PHƯƠNG TRÌNH VI PHÂN CÓ SẴN
TRONG PHẦN MỀM MATLAB.
Trong phần mềm Matlab có sẵn rất nhiều thuật giải hệ phương trình vi phân, thí
dụ ode45; ode23; ode113; ode15s; ode23s; ode23t; ode23tb. Đây là các thuật giải
phương trình vi phân có độ chính xác rất cao và dễ sử dụng. Khi sử dụng chúng, ta
không cần đưa ra bước tích phân, mà thuật giải tự chọn bước tích phân phù hợp để
đảm bảo độ chính xác. Đối với xác phương trình vi phân tuyến tính, thuật giải tự chọn
bước tích phân thay đổi. Cú pháp của các thuật giải như sau:
[T,Y] = solver(odefun,tspan,y0);
hoặc
[T,Y] = solver(odefun,tspan,y0,options);
hoặc
[T,Y] = solver(odefun,tspan,y0,options,p1,p2...);
hoặc
[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