Đề tài 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ụ . 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 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 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 . Hàm truyền trên tương ứng với hệ hai phương trình bậc nhất trong đó các biến số y1(t), y2(t) có nghĩa như sau: y1(t)=y(t)-toạ độ ra; ; 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 tại bước tích phân thứ i được xác định như sau 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); 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 . Hàm truyền trên tương ứng với phương trình vi phân . Để dựng đặc trưng quá độ của khâu tích phân cần bổ xung thêm phương trình (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ó ; ; ; ; Vì vậy (1) 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 (đố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; 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 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 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 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 (2) 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 . Hàm truyền trên tương ứng với phương trình vi phân . Giá trị của hàm y(t) tại bước tích phân thứ i được tính như sau . 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; 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] 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 Với điều kiện ban đầu y1(0)=0; y2(0)=1; y3(0)=1, thời gian tích phân 12 giây. Trước hết, ta tạo file chương trình con function dy = rigid(t,y) dy = zeros(3,1); % a column vector (số 3 chỉ số phương trình) dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2); Giả sử ta sử dụng options=odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]); thuật giải ode45. Cú pháp của thuật giải có dạng sau [T,Y] = ode45('rigid',[0 12],[0 1 1],options); plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.'); Thí dụ 5. Sử dụng thuật giải phương trình vi phân trong Matlab để dựng đặc trưng quá độ của khâu quán tính trong thí dụ 3. File chương trình con có dạng: function dy = rigid1(t,y) T=0.1; K=1; a=1/T; 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]); plot(T,Y(:,1),'-',T,Y(:,2),'-.'); hold on; grid on;