Bai 3 (cont 2)

4. Tính đạo hàm, tích phân số (Numerical Differentiation and Numerical Integration)

• Có thể dùng Matlab để tính đạo hàm số. Ví dụ tính đạo hàm số của hàm sin:

» x=0:0.01:2*pi;

» y=sin(x); % đường đồ thị màu xanh

» dydx=diff(y)./diff(x); % đường đồ thị màu đỏ đạo hàm của hàm y = sin(x)

• Cũng có thể tính toán đạo hàm số trên ma trận như sau:

» mat=[1 3 5;4 8 6];

» dm=diff(mat,1,2) % Tính đạo hàm bậc 1 theo hàng

• Tính gradient 2 chiều của một hàm 2 biến bằng cách nhập giá trị x,y theo dạng ma trận rồi dùng lệnh grafient.

» [dx,dy]=gradient(mat);

Phép tích phân số

• Matlab chứa những hàm tính tích phân số theo các phương pháp phổ biến

  • Phương pháp hình vuông của Simpson (tham số đầu vào là một hàm )

» q=quad('myFun',0,10); % q là tích phân của hàm myFun từ 0 đến 10

» q2=quad(@(x) sin(x)*x,0,pi) % q2 là tích phân của sin(x)*x từ 0 đến pi

• • Phương pháp hình thang (tham số đầu vào là một vector)

» x=0:0.01:pi;

» z=trapz(x,sin(x)); % z là tích phân của sin(x) từ 0 đến pi

» z2=trapz(x,sqrt(exp(x))./x)

5. Phương trình vi phân

• • Giải phương trình vi phân (tích phân phương trình vi phân)

  • Cho một phương trình vi phân, ta có thể dùng phương pháp lặp để tìm nghiệm. Phương pháp này:

    • Tính đạo hàm tại mỗi điểm và xấp xỉ quỹ đạo nghiệm bằng những đoạn thẳng

    • Sẽ xuất hiện sai lệch tích lũy

    • Nếu sử dụng phương pháp tìm bước nhảy linh hoạt (variable step) sẽ giúp giảm số lần lặp cần thiết

• • • Matlab hỗ trợ những công cụ tìm nghiệm phương trình vi phân (các ODE Solvers) theo các thuật toán khác nhau

      • Sử dụng thuật toán giải phương trình vi phân đúng cách có thể giúp bạn tiết kiệm thời gian và cho kết quả chính xác hơn. Các ODE solvers thông dụng trong MATLAB gồm:

» ode23 % Phương pháp lặp bậc thấp, sử dụng hiệu quả khi cần lấy tích phân trên một khoảng nhỏ

% và cần tốc độ nhanh hơn là độ chính xác

» ode45 % Lặp bậc cao (Runge-Kutta). Độ chính xác cao hơn và tốc độ vừa phải. Đây là giải thuật

» ode15s % Thuật toán của Gear, dùng khi có phương trình vi phân có hằng số thời gian thay đổi

• • Cú pháp câu lệnh giải phương trình vi phân

• • • Đầu vào:

      • Tên hàm chứa phương trình vi phân muốn giải (được lưu trong 1 mfile khác). Hàm này lấy đầu vào (t,y) và trả về dy/dt

      • [0,10] : Khoảng muốn tích phân phương trình, bắt đầu từ t = 0 và kết thúc tại t = 10.

      • [1;0] : Vector cột chứa điều kiện đầu của phương trình vi phân (t0 = 1, y(0) = 0)

    • Đầu ra:

      • t

      • y = y(t)

• • ODE function

    • là hàm chứa phương trình vi phân cần giải.

    • trả về giá trị đạo hàm của hàm tại các thời điểm

    • Ví dụ: giả sử có một hệ 2 phương trình vi phân mô tả một phản ứng hóa học A<->B

dA/dt= -10A + 50B

dB/dt= 10A - 50B

• • • Lập ODE function thể hiện hệ phương trình vi phân trên như sau:

• • • Giải phương trình vi phân

» [t,y]=ode45('chem',[0 0.5],[0 1]); % Coi rằng chất B tồn tại lúc đầu

» plot(t,y(:,1),'k','LineWidth',1.5);

» hold on;

» plot(t,y(:,2),'r','LineWidth',1.5);

» legend('A','B');

» xlabel('Time (s)');

» ylabel('Amount of chemical (g)');

» title('Chem reaction');

• • Phương trình vi phân bậc cao hơn

  • Ta phải biến đổi phương trình vi phân bậc cao về hệ phương trình vi phân bậc nhất để có thể sử dụng hàm ode45

    • Chú ý: hàm ode45 có thể giải cả các phương trình vi phân phi tuyến

  • Ví dụ về con lắc

    • Cách đưa bài toán con lắc về dạng ODE function bậc 1 như sau

• • Giải phương trình vi phân và vẽ kết quả thu được

% Tính toán vị trí và vận tốc của con lắc

» [t,x]=ode45('pendulum',[0 10],[0.9*pi 0]);

» plot(t,x(:,1));

» hold on;

» plot(t,x(:,2),'r');

» legend('Position','Velocity');

% Biểu diễn trên mặt phẳng pha

» plot(x(:,1),x(:,2));

» xlabel('Position');

» yLabel('Velocity');% Mặt phẳng pha biểu diễn sự phụ thuộc của biến trạng thái này

% đối với biến trạng thái khác

• • Các lựa chọn với các hàm ode

    • Công cụ ODE của Matlab mặc định dùng bước thời gian có thể thay đổi (variable step), tuy nhiên cũng có thể dùng bước cố định (fixed step)

» [t,y]=ode45('chem',[0:0.001:0.5],[0 1]); % Chỉ rõ bước nhảy bằng cách đưa một vector thời gian

% Giải phương trình vi phân với bước cố định thường chậm hơn so với bước linh hoạt

• • Bạn có thể tùy chỉnh sai số cho phép bằng cách sử dụng odeset . Ví dụ

» options=odeset('RelTol',1e-6,'AbsTol',1e-10);

» [t,y]=ode45('chem',[0 0.5],[0 1],options);

% Câu lệnh trên đảm bảo sai số ở mỗi bước là nhỏ hơn RelTol lần giá trị tại bước nhảy đó

% và đồng thời sai số đó nhỏ hơn AbsTol

• Gõ doc odeset để biết thêm thông tin về các tùy chọn trong giải phương trình vi phân

Bài tập

Dùng lệnh ode45 để tính cho y(t) trong khoảng t = [0 10] với điều kiện khởi tạo y(0) =10 và dy/dt = - t*y/10

Vẽ kết quả thu được.

Hướng dẫn

Cách 1:

Tạo một hàm như sau trong mfile odefun.m:

function dydt=odefun(t,y)

dydt=-t*y/10;

Tính và vẽ kết quả trên Command Window

» [t,y]=ode45(‘odefun’,[0 10],10);

» plot(t,y);xlabel('Time');ylabel('y(t)');

Cách 2: không cần tạo hàm trước, thực hiện như sau:

» [t,y]=ode45(@(t,y) –t*y/10,[0 10],10);

» plot(t,y);xlabel('Time');ylabel('y(t)');

Cả 2 cách đều cho kết quả sau

<< Previous Page Bài 4 >>