Chương 9: Hệ các PVT

Trở về Mục lục cuốn sách

Ma trận

Ma trận là dạng hai chiều của một véc-tơ. Cũng như véc-tơ, ma trận gồm có các phần tử được phân biệt bởi chỉ số. Sự khác biệt với véc-tơ là ở chỗ các phần tử ma trận được xếp theo hàng và cột, vì vậy cần có hai chỉ số để xác định được một phần tử.

Một trong những cách tạo ra ma trận là dùng hàm magic, vốn trả lại một ma trận kì ảo với kích cỡ cho trước:

>> M = magic(3)

M =  8     1     6
     3     5     7
     4     9     2

Nếu không biết kích cỡ của ma trận, bạn có thể dùng whos để hiển thị nó:

>> whos
  Name        Size                    Bytes  Class
  M           3x3                        72  double array

Hoặc dùng hàm size, lần này sẽ trả về một véc-tơ:

>> V = size(M)

V = 3     3

Phần tử đầu tiên của véc-tơ là số hàng, và phần tử thứ hai là số cột.

Để đọc một phần tử của ma trận, bạn cần chỉ ra số thứ tự hàng và cột:

>> M(1,2)

ans = 1

>> M(2,1)

ans = 3

Khi bắt đầu làm quen với ma trận, bạn cần nhớ được chỉ số nào đi trước, hàng hay cột. Tôi thì thường nói nhẩm “hàng, cột”, như một câu thần chú. Bạn cũng có thể nhớ kiểu “xuống dưới, sang ngang,” hoặc chữ tắt HC.

Một cách khách để tạo ra ma trận là viết các phần tử giữa cặp ngoặc vuông, cùng dấu chấm phẩy giữa các hàng:

>> D = [1,2,3 ; 4,5,6]

D =  1     2     3
     4     5     6

>> size(D)

ans = 2     3

Véc-tơ hàng và cột

Mặc dù có thể dễ nhớ hơn theo các loại số vô hướng, véc-tơ và ma trận, song theo quan điểm của MATLAB, mọi thứ đều là ma trận. Một số vô hướng chính là một ma trận có một hàng và một cột:

>> x = 5;
>> size(x)

ans = 1     1

Và véc-tơ là ma trận có một hàng:

>> R = 1:5;
>> size(R)

ans = 1     5

Thực ra có hai loại véc-tơ. Loại véc-tơ ta đã gặp cho đến giờ là véc-tơ hàng, vì các phần tử được sắp trên một hàng; còn loại kia là véc-tơ cột, vì các phần tử xếp theo một cột.

Một cách tạo ra véc-tơ cột là lập một ma trận chỉ với một phần tử trên mỗi hàng:

>> C = [1;2;3]

C =

     1
     2
     3

>> size(C)

ans = 3     1

Sự khác biệt giữa véc-tơ hàng và cột là quan trọng trong môn đại số tuyến tính, nhưng với đa số các phép tính với véc-tơ thì điều đó chẳng quan trọng. Khi bạn đánh chỉ số cho véc-tơ, bạn không cần biết đó là loại véc-tơ gì:

>> R(2)

ans = 2

>> C(2)

ans = 2

Toán tử chuyển vị

Toán tử chuyển vị, vốn trông giống như dấu nháy đơn, được dùng để tính chuyển vị của ma trận, vốn là một ma trận mới có tất cả những phần tử của ma trận ban đầu, nhưng với mỗi hàng được dựng thành một cột (hoặc bạn có thể nghĩ ngược lại).

Ở ví dụ sau:

>> D = [1,2,3 ; 4,5,6]

D =  1     2     3
     4     5     6

D có hai hàng, vì vậy ma trận chuyển vị của nó có hai cột:

>> Dt = D'

Dt = 1     4
     2     5
     3     6

Toán tử chuyển vị có tác động gì đến véc-tơ hàng, véc-tơ cột, và số vô hướng?

Lotka-Voltera

Mô hình Lotka-Voltera mô tả sự tương tác giữa hai loài, vật săn và con mồi, trong cùng một hệ sinh thái. Một ví dụ là hai loài thỏ và cáo.

Mô hình được diễn tả bởi hệ phương trình vi phân sau:

Rt = aR – bRF

Ft = ebRF – cF
trong đó

R là số con thỏ,

F là số con cáo,

a là tốc độ tăng trưởng của thỏ khi không có loài săn đuổi,

c là tốc độ chết của cáo khi không có mồi,

b là tốc độ chết của thỏ ứng với mỗi tương tác với cáo,

e là hiệu suất biến đổi số thỏ bị ăn thịt thành cáo.

Thoạt nhìn, bạn có thể nghĩ rằng có thể giải ngay các phương trình bằng cách gọi ode45 một lần để giải ra R như một hàm theo thời gian và một lần khác để giải F. Vấn đề là mỗi phương trình đều chứa đến cả hai biến, cũng vì lý do này mà ta gọi là hệ phương trình chứ không phải là một loạt các phương trình riêng lẻ. Để giải một hệ, bạn cần phải giải tất cả phương trình cùng lúc.

Thật may là ode45 có thể giải được hệ phương trình. Sự khác biệt ở đây là điều kiện ban đầu có dạng một véc-tơ chứa những giá trị ban đầu là R(0)F(0), và kết quả đầu ra là ma trận có chứa một cột dành cho R và một cột cho F.

Sau đây là các hàm tốc độ được viết với các tham số a = 0. 1, b = 0. 01, c = 0. 1 and e = 0. 2:

function res = lotka(t, V)
    % thao go cac phan tu cua V
    r = V(1);
    f = V(2);

    % dat cac tham so
    a = 0.1;
    b = 0.01;
    c = 0.1;
    e = 0.2;

    % tinh cac dao ham
    drdt = a*r - b*r*f;
    dfdt = e*b*r*f - c*f;

    % xep cac gia tri dao ham vao vec-to
    res = [drdt; dfdt];
end

Như thường lệ, biến đầu vào thứ nhất là thời gian. Biến đầu vào thứ hai là một véc-tơ V gồm hai phần tử, R(t)F(t). Tôi viết nó bằng chữ in để dễ gợi nhớ rằng đó là véc-tơ. Phần thân hàm gồm bốn đoạn, mỗi đoạn đều được chú thích kèm theo.

Đoạn thử nhất nhằm tháo gỡ véc-tơ bằng cách sao chép các phần tử ra các biến vô hướng. Điều này không nhất thiết phải làm, nhưng việc đặt tên cho từng giá trị giúp tôi nhớ được ý nghĩa của chúng. Nó cũng giúp ích ở đoạn thứ ba, khi ta tính đạo hàm, thể hiện ở cách viết giống như biểu thức toán học ban đầu, góp phần tránh được lỗi.

Đoạn thứ hai đặt các tham số miêu tả các tốc độ sinh sản của thỏ và cáo, cùng với các đặc trưng tương tác giữa chúng. Nếu ta đang nghiên cứu một hệ thực sự, thì những giá trị trên sẽ phải được rút ra từ quá trình quan sát đời sống của loài vật; nhưng ở ví dụ này tôi chỉ chọn các giá trị sao cho thu được kết quả thú vị.

Đoạn cuối cùng gói các đạo hàm tính được trở lại vào véc-tơ. Khi ode45 gọi hàm này, nó cung cấp đầu vào là một véc-tơ và trông đợi đầu ra cũng là một véc-tơ.

Bạn đọc tinh mắt có thể nhận ra có điểm khác biệt ở dòng:

    res = [drdt; dfdt];

Dấu chấm phẩy ngăn cách giữa hai phần tử của véc-tơ này không có gì sai. Nó rất cần trong trường hợp này vì ode45 yêu cầu kết quả của hàm này phải là một véc-tơ cột.

Bây giờ ta có thể chạy ode45 như sau:

ode45(@lotka, [0, 365], [100, 10])

Như thường lệ, thông số thứ nhất là một chuôi hàm, thứ hai là khoảng thời gian, và thứ ba là điều kiện ban đầu. Ở đây, điều kiện ban đầu là một véc-tơ: phần tử thứ nhất là số con thỏ lúc t = 0, phần tử thứ hai là số con cáo.

Thứ tự của hai phần tử này (số thỏ và số cáo) là tùy vào bạn, nhưng phải thống nhất. Nghĩa là điều kiện ban đầu được cung cấp khi gọi đến ode45 phải có cùng thứ tự như trong lotka, nơi bạn gỡ véc-tơ đầu vào và gói lại véc-tơ đầu ra. MATLAB không hiểu ý nghĩa các giá trị này; tất cả sự theo dõi đều là trách nhiệm của người lập trình.

Nhưng nếu bạn xếp đúng thứ tự thì kết quả sẽ được như sau:

Đồ thị Lotka-Volterra

Trục x là thời gian tính theo ngày; trục y là số cá thể. Đường cong phía trên chỉ số thỏ; đường phía dưới chỉ số cáo. Kết quả này là một trong số những dạng mẫu mà hệ này có thể rơi vào, tùy theo các điều kiện ban đầu và các tham số. Bạn hãy tự thử nghiệm với các giá trị khác nhau, như một bài tập cho riêng mình.

Điều gì có thể gây trục trặc?

Véc-tơ đầu ra của hàm tốc độ phải là một véc-tơ cột; nếu không bạn sẽ nhận được

??? Error using ==> funfun/private/odearguments
LOTKA must return a column vector.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs,
odeFcn, ...

vốn là một lời thông báo lỗi khá rõ ràng. Chỉ có điều là tại sao phải cần một véc-tơ cột, song đó không phải là điều mà ta quan tâm.

Một lỗi khác có thể xảy ra khác là việc đảo lộn thứ tự của các phần tử trong số các điều kiện đầu, hay là các véc-tơ bên trong lotka. Một lần nữa, MATLAB không biết rằng các phần tử có ý nghĩa gì, cho nên nó không thể phát hiện ra lỗi kiểu này, mà chỉ cho kết quả sai.

Ma trận kết quả

Như ta đã thấy ở trên, nếu bạn gọi ode45 mà không gán kết quả vào biến thì nó sẽ vẽ đồ thị kết quả đó. Nếu gán kết quả vào biến thì việc làm này sẽ ngăn cản việc hiển thị hình vẽ. Câu lệnh cụ thể như sau:

>> [T, M] = ode45(@lotka, [0, 365], [100, 10]);

Bạn có thể hình dung vế trái của lệnh gán này như một véc-tơ chứa các biến.

Cũng như trong các ví dụ trước, T là véc-tơ chứa các giá trị thời điểm mà tại đó ode45 thực hiện việc ước tính. Nhưng khác với các ví dụ trước, biến đầu ra thứ hai là một ma trận chứa một cột cho mỗi biến (trong trường hợp này là RF) và một hàng cho mỗi giá trị thời gian.

>> size(M)

ans = 185     2

Cấu trúc này—một cột cho mỗi biến—là cách thông thường để dùng ma trận. plot hiểu được cấu trúc này, nên khi bạn viết:

>> plot(T, M)

MATLAB sẽ hiểu rằng nó cần vẽ đồ thị cho mỗi cột của M theo T.

Bạn có thể sao chép các cột của M ra các biến khác như sau:

>> R = M(:, 1);
>> F = M(:, 2);

Ở đây, dấu hai chấm biểu thị cho khoảng từ 1 đến end, vì vậy M(:, 1) nghĩa là “tất cả các hàng thuộc cột 1” và M(:, 2) nghĩa là “tất cả các hàng thuộc cột 2.”

>> size(R)

ans = 185     1

>> size(F)

ans = 185     1

Như vậy RF là các véc-tơ cột.

Nếu vẽ đồ thị các véc-tơ này theo nhau, kiểu như:

>> plot(R, F)

thì bạn sẽ thu được một đồ thị pha giống hình dưới đây:

đồ thị pha

Mỗi điểm trên đồ thị này biểu diễn một số thỏ nhất định (theo trục x) và một số cáo nhất định (theo trục y).

Vì chỉ có hai biến trong hệ này, mỗi điểm trên mặt phẳng biểu thị cho trạng thái đầy đủ của hệ.

Theo thời gian, trạng thái sẽ di chuyển quanh mặt phẳng: hình vẽ này cho thấy đường đi vạch nên bởi trạng thái trong khoảng thời gian tính toán. Đường đi này được gọi là quỹ đạo.

Vì biểu hiện của hệ thống này có tính chu kì nên quỹ đạo thu được là một đường cong khép kín.

Nếu có 3 biến trong hệ, ta sẽ cần 3 chiều không gian để biểu thị trạng thái của hệ, vì vậy quỹ đạo là một đường cong 3 chiều. Bạn có thể dùng plot3 để vạch nên đường quỹ đạo trong không gian 3 chiều, nhưng khi có từ 4 biến trở lên thì bạn chẳng còn cách biểu thị trực quan nào nữa.

Thuật ngữ

véc-tơ hàng:
Trong MATLAB, một ma trận chỉ có một hàng.

véc-tơ cột:
Một ma trận chỉ có một cột.

chuyển vị:
Phép tính để chuyển các hàng của ma trận thành cột (hoặc ngược lại).

hệ phương trình:
Tập hợp các phương trình chứa một tập hợp các biến sao cho các phương trình liên hệ qua lại lẫn nhau.

đoạn:
Một loạt các lệnh liên tiếp hình thành nên một phần của hàm, thường được kèm theo bởi một lời giải thích rõ ràng.

tháo gỡ:
Sao chép các phần tử trong một véc-tơ vào một tập hợp các biến.

đóng gói:
Sao chép các giá trị từ một tập hợp các biến vào một véc-tơ.

trạng thái:
Nếu một hệ có thể được miêu tả bởi một tập hợp các biến thì giá trị của các biến đó được gọi là trạng thái của hệ.

đồ thị pha:
Đồ thị cho thấy trạng thái của hệ như một điểm trên không gian các trạng thái hiện có.

quỹ đạo:
Đường đi của đồ thị pha cho thấy trạng thái của một hệ thay đổi như thế nào theo thời gian.

Bài tập

Dựa vào các ví dụ ta đã thấy, bạn có thể nghĩ rằng tất cả PVT mô tả số cá thể đều có dạng hàm theo thời gian, nhưng thực ra không phải vậy.

Theo Wikipedia1, “điểm hấp dẫn Lorenz, được Edward Lorenz đề xuất vào năm 1963, là một hệ động lực tất định ba chiều phi tuyến được suy từ các phương trình giản hóa đặc trưng cho các cuộn đối lưu từ phương trình động lực khí quyển. Với một tập hợp nhất định các tham số, hệ thống sẽ có biểu hiện hỗn loạn và cho thấy cái mà ngày nay ta gọi là điểm hấp dẫn dị thường…”

Hệ được mô tả bởi hệ phương trình vi phân sau:
xt = σ(y – x)

yt = x(r – z) – y

zt = xy – bz
Các giá trị thông dụng của tham số gồm σ = 10, b = 8/3r = 28.

Hãy dùng ode45 để ước tính một nghiệm của hệ phương trình này.

  1. Bước đầu tiên là viết một hàm có tên là lorenz nhận tV làm các biến đầu vào, trong đó các thành phần của V được hiểu là các giá trị hiện thời của x, yz. Nó cần phải tính được các đạo hàm tương ứng và trả chúng về dưới dạng một véc-tơ cột.
  2. Bước tiếp theo là thử nghiệm hàm của bạn bằng cách gọi nó từ dòng lệnh với các giá trị như t = 0, x = 1, y = 2z = 3. Một khi hàm chạy được, bạn cần biến nó thành một hàm lặng trước khi gọi ode45.
  3. Giả định rằng Bước 2 đã xong, bạn có thể dùng ode45 để ước tính nghiệm cho khoảng thời gian từ t0 = 0 đến te = 30 với điều kiện đầu x = 1, y = 2z = 3.
  4. Hãy dùng plot3 để vẽ đồ thị quỹ đạo theo x, yz.
Advertisements

4 phản hồi

Filed under Mô hình hóa, Sách

4 responses to “Chương 9: Hệ các PVT

  1. Pingback: Chương 1. Các biến và giá trị | Blog của Chiến

  2. Pingback: Chương 11: Tối ưu hóa và nội suy | Blog của Chiến

  3. Pingback: Chương 6: Tìm nghiệm | Blog của Chiến

  4. Pingback: Mô hình hóa hiện tượng vật lý bằng MATLAB | Blog của Chiến

Trả lời

Mời bạn điền thông tin vào ô dưới đây hoặc kích vào một biểu tượng để đăng nhập:

WordPress.com Logo

Bạn đang bình luận bằng tài khoản WordPress.com Đăng xuất / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Đăng xuất / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Đăng xuất / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Đăng xuất / Thay đổi )

Connecting to %s