Chương 12: Bây giờ véc-tơ mới thật là véc-tơ

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

Véc-tơ là gì?

Từ “véc-tơ” có thể mang những nghĩa khác nhau đối với từng người. Trong MATLAB, véc-tơ là một ma trận chỉ có một hàng, hoặc một cột. Cho đến giờ, ta đã dùng các véc-tơ của MATLAB để biểu diễn:

dãy:
Dãy là một tập hợp các giá trị được nhận diện bởi các chỉ số nguyên; theo cách làm tự nhiên ta có thể lưu các phần tử của dãy như những phần tử của một véc-tơ trong MATLAB.

véc-tơ trạng thái:
Véc-tơ trạng thái là một tập hợp các giá trị để mô tả trạng thái của một hệ vật lý. Khi gọi ode45, bạn cho nó các điều kiện ban đầu dưới dạng một véc-tơ trạng thái. Sau đó, khi ode45 gọi hàm tốc độ mà bạn lập nên, nó sẽ trả kết quả là một véc-tơ trạng thái khác.

ánh xạ rời rạc:
Nếu có trong tay hai véc-tơ cùng độ dài, bạn có thể hình dung chúng như một phép ánh xạ từ những phần tử của một véc-tơ này sang các phần tử thuộc véc-tơ kia. Chẳng hạn, ở Mục {Chuột}, kết quả thu được từ ode45 là các véc-tơ, TY; chúng biểu diễn một phép ánh xạ từ các giá trị thời gian của T sang các giá trị số lượng chuột có trong Y.

Trong chương này ta sẽ xét đến một công dụng khác của véc-tơ trong MATLAB: để biểu diễn các véc-tơ không gian. Một véc-tơ không gian là một giá trị nhằm biểu diễn một đại lượng vật lý nhiều chiều, như vị trí, vận tốc, gia tốc, hoặc lực1.

Các đại lượng nói trên không thể diễn tả bởi một con số vì chúng có chứa nhiều thành phần. Chẳng hạn, trong không gian tọa độ Đề-các 3 chiều, ta cần phải có 3 con số mới xác định được một vị trí trong không gian, mà thường gọi là các tọa độ x, yz. Một ví dụ khác, trong không gian tọa độ cực 2 chiều, bạn có thể xác định một vận tốc bằng hai con số, một là độ lớn và thứ hai là góc, thường được gọi là rθ.

Biểu diễn các véc-tơ không gian bằng véc-tơ trong MATLAB rất tiện vì MATLAB biết cách thực hiện hầu hết các phép tính với véc-tơ phục vụ cho việc mô phỏng vật lý. Chẳng hạn, giả sử bạn có vận tốc của quả bóng chày theo dạng véc-tơ MATLAB với hai phần tử, vxvy, chính là hai thành phần của vận tốc theo các phương xy.

>> V = [30, 40]        % velocity in m/s

Và giả sử rằng bạn cần phải tính gia tốc tổng hợp của quả bóng dưới tác dụng của lực cản không khí và trọng lực. Theo ký hiệu toán học, lực cản được tính bởi

Fd =  –½ ρ ⋅ v2 ⋅ A ⋅ Cd ⋅ 

trong đó V là véc-tơ không gian biểu diễn vận tốc, v là độ lớn của vận tốc (còn được gọi là “tốc độ”), và là một véc-tơ đơn vị theo hướng của véc-tơ vận tốc. Các đại lượng khác, ρ, ACd đều là số vô hướng.

Độ lớn của véc-tơ thì bằng căn bậc hai của tổng các bình phương của từng phần tử. Bạn có thể tính với hypotenuse ở Mục {Cạnh huyền}, hoặc dùng hàm MATLAB có tên norm (norm, hay chuẩn, là tên gọi khác 2 cho độ lớn của véc-tơ):

>> v = norm(V)

v = 50

là một véc-tơ đơn vị, nghĩa là nó phải có chuẩn bằng 1, và chỉ theo cùng hướng với V. Cách dễ nhất để tính véc-tơ đơn vị này là chia V cho chuẩn của chính nó.

>> Vhat = V / v

Vhat =  0.6  0.8

Bây giờ ta có thể đảm bảo rằng chuẩn của bằng 1:

>> norm(Vhat)

ans = 1

Để tính Fd ta chỉ cần nhân các đại lượng vô hướng với .

Fd = - 1/2 * C * rho * A * v^2 * Vhat

Tương tự như vậy, ta có thể tính gia tốc bằng cách chia véc-tơ Fd cho số vô hướng m.

Ad = Fd / m

Để biểu diễn gia tốc trọng trường, ta tạo nên một véc-tơ gồm hai phần tử:

Ag = [0; -9.8]

Phần tử x của gia tốc trọng trường bằng 0; phần tử y bằng  –9,8 m/s2.

Cuối cùng ta tính gia tốc tổng cộng bằng cách cộng véc-tơ lại:

A = Ag + Ad;

Một điều hay trong cách tính này là ta không cần phải nghĩ nhều về từng thành phần của véc-tơ. Bằng cách coi véc-tơ không gian như những đại lượng cơ bản, ta có thể diễn tả các phép tính phức tạp một cách ngắn gọn.

Tích vô hướng và tích hữu hướng

Việc nhân một véc-tơ với số vô hướng là quá rõ ràng; cộng hai véc-tơ cũng vậy. Nhưng nhân hai véc-tơ thì tinh vi hơn. Hóa ra rằng có hai phép toán véc-tơ trông tựa như phép nhân, đó là tích vô hướngtích hữu hướng.

Tích vô hướng của hai véc-tơ AB cho kết quả là số vô hướng:

d = acosθ
trong đó a là độ lớn của A, b là độ lớn của B, còn θ là góc giữa hai véc-tơ. Ta đã biết cách tính các độ lớn, và có thể hình dung ra cách tính θ, nhưng không cần thiết phải làm việc đó. MATLAB có một hàm, dot, để tính tích vô hướng.

d = dot(A, B)

dot tính được với số chiều bất kì, miễn là AB có cùng số phần tử.

Nếu một trong hai toán hạng là véc-tơ đơn vị thì bạn có thể dùng tích vô hướng để tính ra thành phần của véc-tơ A theo hướng véc-tơ đơn vị î đó:

s = dot(A, ihat)

Ở ví dụ này, shình chiếu vô hướng của A lên phương î. Còn hình chiếu véc-tơ là một véc-tơ có độ lớn s theo phương î:

V = dot(A, ihat) * ihat

Tích hữu hướng của hai véc-tơ AB là một véc-tơ hướng vuông góc với AB và có độ lớn:

c = asinθ
trong đó (một lần nữa) a là độ lớn của A, b là độ lớn của B, còn θ là góc giữa hai véc-tơ. MATLAB có một hàm, cross, để tính tích hữu hướng.

C = cross(A, B)

cross chỉ tính được với véc-tơ 3 chiều; và kết quả cũng là một véc-tơ 3 chiều khác.

Một cách hay dùng của cross là để tính các mô-men lực. Nếu bạn biểu diễn tay đòn mô-men R và lực F dưới dạng các véc-tơ 3 chiều, thì mô-men lực sẽ được tính một cách đơn giản:

Tau = cross(R, F)

Nếu các thành phần của R được tính bằng mét và các thành phần của F tính bằng newton, thì các thành phần trong mô-men lực Tau được tính bằng newton-mét.

Cơ học thiên thể

Việc mô hình hóa cơ học thiên thể là một dịp tốt để ta thực hiện tính toán với véc-tơ không gian. Hãy tưởng tượng một ngôi sao với khối lượng m1 tại một điểm trong không gian được mô tả bởi véc-tơ P1, và một hành tinh khối lượng m2 tại điểm P2. Độ lớn của lực hấp dẫn3 giữa chúng là

fg = G(m1m2) / r2
trong đó r là khoảng cách giữa chúng, còn G là hằng số hấp dẫn, bằng khoảng 6. 67 × 10 - 11Nm2 / kg2. Nhớ rằng đây là giá trị của G chỉ trong trường hợp ta tính khối lượng theo ki-lô-gam, khoảng cách theo mét, và lực theo newton.

Hướng của lực lên ngôi sao tại P1 là chỉ về phía P2. Ta có thể tính hướng tương đối giữa chúng bằng cách trừ véc-tơ; nếu ta tính R = P2 - P1, thì hướng của R bây giờ sẽ chỉ từ P1 đến P2.

Khoảng cách giữa hành tinh và ngôi sao là chiều dài của R:

r = norm(R)

Hướng của lực tác dụng lên ngôi sao là :

rhat = R / r

Hãy viết một dãy lệnh MATLAB để tính F12 là véc-tơ biểu diễn lực mà hành tinh tác dụng lên ngôi sao, và F21, là lực do ngôi sao tác dụng lên hành tinh.

Hãy gói các câu lệnh sau vào một hàm có tên gravity_force_func nhận vào các biến P1, m1, P2, và m2 rồi trả lại F12.

Hãy viết một chương trình mô phỏng quỹ đạo của Sao Mộc quanh Mặt Trời. Khối lượng của Mặt Trời vào khoảng 2,0 × 1030 kg. Bạn có thể lấy số liệu khối lượng của Sao Mộc cũng như khoảng cách đến Mặt Trời và vận tốc của Sao Mộc trên quỹ đạo từ trang http://en.wikipedia.org/wiki/Jupiter. Chạy chương trình để đảm bảo rằng Sao Mộc quay trọn 1 vòng quỹ đạo quanh Mặt Trời hết khoảng 4332 ngày.

Tạo hình chuyển động

Hình động là một công cụ hữu ích để kiểm tra kết quả tính toán từ một mô hình vật lý. Nếu có gì trục trặc xảy ra, mọi sự sẽ rõ ràng trên hình động. Có hai cách làm hình động trong MATLAB. Một cách là dùng getframe để chụp một loạt các ảnh và movie để phát lại chúng. Một cách làm khác, không chính thức, là vẽ một loạt các hình đồ thị. Sau đây là một ví dụ tôi viết cho Bài tập {Sao Mộc}:

function animate_func(T,M)
    % animate the positions of the planets, assuming that the
    % columns of M are x1, y1, x2, y2.
    X1 = M(:,1);
    Y1 = M(:,2);
    X2 = M(:,3);
    Y2 = M(:,4);

    minmax = [min([X1;X2]), max([X1;X2]), min([Y1;Y2]), max([Y1;Y2])];

    for i=1:length(T)
        clf;
        axis(minmax);
        hold on;
        draw_func(X1(i), Y1(i), X2(i), Y2(i));
        drawnow;
    end
end

Các biến đầu vào là kết quả đầu ra từ ode45, một véc-tơ T và một ma trận M. Các cột của M chứa các vị trí và vận tốc của Mặt Trời và Sao Mộc, vì vậy X1Y1 nhận các tọa độ của Mặt Trời; còn X2Y2 nhận các tọa độ của Sao Mộc.

minmax là một véc-tơ gồm 4 phần tử được dùng bên trong vòng lặp để thiết lập các trục trên biểu đồ. Thông tin này cần thiết bởi vì nếu không, MATLAB sẽ tự co dãn hình mỗi lần qua một lượt lặp, và các trục sẽ luôn thay đổi, dẫn đến khó theo dõi hình động.

Qua mỗi lượt lặp, animate_func dùng clf để xóa hình và axis để đặt lại các trục. hold on giúp ta có thể vẽ được nhiều đường đồ thị lên cùng một hệ trúc (nếu không MATLAB sẽ tự xóa hình mỗi lần bạn gọi câu lệnh plot mới).

Cũng qua mỗi lượt lặp, ta phải gọi drawnow để MATLAB thực sự hiển thị từng hình vẽ. Còn nếu không, nó sẽ đợi đến tận khi bạn vẽ xong hết các hình rồi mới cập nhật quá trình hiển thị.

draw_func là hàm số đảm nhiệm việc vẽ đồ thị:

function draw_func(x1, y1, x2, y2)
    plot(x1, y1, 'r.', 'MarkerSize', 50);
    plot(x2, y2, 'b.', 'MarkerSize', 20);
end

Các biến đầu vào là vị trí của Mặt Trời và Sao Mộc. draw_func dùng plot để vẽ Mặt Trời như một dấu đỏ lớn và Sao Mộc như một dấu xanh lam nhỏ hơn.

Để chắc rằng bạn hiểu được cách hoạt động của animate_func, hãy thử chú thích để loại ra một số dòng lệnh xem điều gì xảy ra.

Một hạn chế của kiểu hình động này là tốc độ của hình động phụ thuộc vào máy của bạn thực hiện việc vẽ hình lúc nào. Vì kết quả từ ode45 thường không xuất ra đều theo thời gian nên hình động có thể sẽ chậm lại khi ode45 nhận bước thời gian ngắn và nhanh lên khi bước thời gian dài hơn.

Có hai cách sửa vấn đề này:

  1. Khi gọi ode45 bạn có thể chỉ định một véc-tơ chứa các thời điểm cần phải tính kết quả. Sau đây là một ví dụ:
    end_time = 1000;
    step = end_time/200;
    [T, M] = ode45(@rate_func, [0:step:end_time], W);
    

    Đối số thứ hai là một véc-tơ khoảng số chạy từ 0 đến 1000 với bước chạy quy định bởi step. Vì step bằng end_time/200, nên sẽ có khoảng 200 hàng trong TM (chính xác là 201).

    Tùy chọn này sẽ không ảnh hưởng đến độ chính xác của kết quả; ode45 vẫn dùng biến thời gian để thực hiện ước tính, nhưng rồi nó sẽ nội suy các giá trị trước khi trả lại kết quả.

  2. Bạn có thể dùng pause để chạy hình động theo đúng thời gian thật. Sau khi vẽ từng hình và gọi nó drawnow, bạn có thể tính thời gian từ đó đến hình kế tiếp và dùng pause để đợi:
    dt = T(i+1) - T(i);
    pause(dt);
    

    Một hạn chế của phương pháp này là nó bỏ qua thời gian cần để vẽ hình, vì vậy nó có xu hướng chạy chậm, đặc biệt nếu hình vẽ phức tạp hoặc bước thời gian quá ngắn.

Hãy dùng animate_funcdraw_func để hiển thị kết quả mô phỏng Sao Mộc. Sửa các hàm này để sao cho thời gian một ngày mô phỏng tương đương với 0,001 giây đồng hồ—mỗi vòng quay sẽ mất khoảng 4,3 giây.

Bảo toàn năng lượng

Một cách làm tiện dụng để kiểm tra độ chính xác của một chương trình giải PVT là xem nó có bảo toàn được năng lượng của hệ hay không. Với chuyển động của hành tinh, hóa ra là ode45 không thể bảo toàn được.

Động năng của một vật chuyển động bằng mv2/2; động năng của hệ Mặt Trời thì bằng tổng các động năng của từng hành tinh và của Mặt Trời. Thế năng của Mặt Trời với khối lượng m1 và một hành tinh có khối lượng m2 cách nhau một khoảng r thì bằng:

U = – G(m1m2) / r

Hãy viết một hàm có tên energy_func để nhận vào các kết quả TM từ chương trình mô phỏng Sao Mộc, rồi tính tổng năng lượng (động năng và thế năng) của hệ với mỗi vị trí và vận tốc ước tính được. Vẽ đồ thị của kết quả này như một hàm theo thời gian và đảm bảo rằng nó giảm dần trong quá trình mô phỏng. Hàm cần viết phải tính được sự thay đổi tương đối về năng lượng, độ chênh lệch giữa năng lượng ban đầu và lúc kết thúc, và biểu thị độ chênh này theo phần trăm so với năng lượng ban đầu.

Bạn có thể giảm tốc độ thất thoát năng lượng bằng cách giảm tùy chọn dung sai của ode45 qua việc dùng odeset (xem Mục {Sự kiện}):

options = odeset('RelTol', 1e-5);
[T, M] = ode45(@rate_func, [0:step:end_time], W, options);

Tên của tùy chọn cần chỉnh là RelTol, viết tắt của “relative tolerance” (dung sai tương đối). Giá trị mặc định là 1e-3 hay 0,001. Các giá trị nhỏ hơn sẽ khiến cho ode45 thêm “khắt khe”, vì vậy cần hàm sẽ cần tính nhiều hơn để giảm nhỏ sai số.

Hãy chạy ode45 với một loạt các giá trị của RelTol và đảm bảo rằng khi dung sai càng nhỏ thì tốc độ thất thoát năng lượng cũng chậm lại.

Hãy chạy chương trình mô phỏng bạn viết với một trong các hàm giải ODE khác có trong MATLAB và xem chúng có bảo toàn năng lượng không.

Mô hình dùng để làm gì?

Trong mục {Mô hình hóa}, tôi đã định nghĩa một “mô hình” như một sự mô tả được giản hóa về một hệ vật lý, và đề cập rằng một mô hình tốt rất thích hợp cho việc phân tích và mô phỏng, và khiến cho các ước tính đủ độ tin cậy để dùng cho mục đích định trước.

Từ đó, ta đã thấy một số ví dụ; bây giờ ta có thể nói thêm những công dụng của mô hình. Các mục đích mà mô hình cần đạt đến thường có xu hướng rơi vào ba dạng sau.

dự đoán:
Một số mô hình có thể dự đoán về các hệ vật lý. Lấy một ví dụ đơn giản, mô hình con vịt ở Bài tập {Vịt} dự đoán về độ cao mà con vịt sẽ nổi. Ở phía đầu kia xét về độ phức tạp, mô hình khí tượng toàn cầu dự đoán tình hình thời tiết hàng mười năm hoặc trăm năm trong tương lai.

thiết kế:
Mô hình rất có ích trong thiết kế kĩ thuật, đặc biệt là để kiểm tra tính khả thi của một phương án thiết kế, và để tối ưu hóa. Chẳng hạn, trong Bài tập {Golf} bạn được yêu cầu phải thiết kế cú đánh golf với sự kết hợp hoàn hảo giữa góc vụt, vận tốc và độ xoáy.

giải thích:
Mô hình có thể giải thích các câu hỏi khoa học. Chẳng hạn, mô hình Lotka-Volterra trong Mục {Lotka-Volterra} đề xuất một lời giải thích cho cơ chế biến động của hệ quần thể động vật dưới hình thức các mối tương tác giữa loài ăn thịt và vật mồi.

Các ví dụ ở cuối chương này sẽ gồm mỗi mô hình theo từng loại kể trên.

Thuật ngữ

véc-tơ không gian:
Một giá trị dùng để biểu thị một đại lượng vật lý gồm nhiều chiều; chẳng hạn vị trí, vận tốc, gia tốc, và lực.

chuẩn:
Độ lớn của một véc-tơ. Đôi khi còn được gọi là “độ dài,” nhưng không được nhầm với số phần tử thuộc một véc-tơ trong MATLAB.

véc-tơ đơn vị:
Véc-tơ có chuẩn bằng 1, dùng để biểu thị hướng.

tích vô hướng:
Số vô hướng, là kết quả phép nhân hai véc-tơ, nó tỉ lệ với độ lớn của từng véc-tơ và cô-sin của góc giữa chúng.

tích hữu hướng:
Véc-tơ là kết quả phép nhân hai véc-tơ khác, với chuẩn tỉ lệ với chuẩn của từng véc-tơ và với sin của góc giữa chúng, và có hướng vuông góc với cả hai véc-tơ.

hình chiếu:
Thành phần của một véc-tơ xét theo hướng của véc-tơ kia (có thể được dùng với nghĩa “hình chiếu vô hướng” hay “hình chiếu véc-tơ”).

Bài tập

Nếu bạn đặt hai bát nước giống hệt nhau vào ngăn đá tủ lạnh, nước ở trong một bát thì ở nhiệt độ bình thường còn ở bát kia thì đang sôi, liệu bát nước nào sẽ đông đá trước?

Gợi ý: bạn có thể sẽ cần tìm hiểu thêm về hiệu ứng Mpemba.

Bạn được yêu cầu phải thiết kế một dường dốc để trượt ván; khác với những dốc trượt thông thường, cái này có thể xoay (kiểu bập bênh) quanh một chốt cố định. Người trượt ván tiến đến đường dốc, trước đó còn đi trên đất bằng, rồi lập tức leo dốc; họ không được phép đặt chân xuống. Nếu họ trượt đủ nhanh thì đường dốc sẽ xoay và họ sẽ chuyển tư thế thành trượt xuống một cách nhẹ nhàng. Sẽ có ban giám khảo chấm điểm kĩ thuật và nghệ thuật cho động tác trượt.

Việc của bạn là thiết kế một dốc trượt để cho phép người trượt hoàn thành bài biểu diễn này, và tạo ra một mô hình vật lý của hệ, một chương trình mô phỏng để tính ra chuyển động của người trượt trên dốc, và hình động biểu diễn kết quả tính được.

Một hệ hai sao bao gồm hai ngôi sao quay xung quanh lẫn nhau và đôi khi có cả những hành tinh quay quanh từng ngôi sao 4. Trong hệ hai sao, một số quỹ đạo là “bền vững” theo nghĩa một hành tinh có thể bay trong quỹ đạo đó mà không bị đâm vào một trong hai sao, và không bị bay tuốt vào không trung.

Mô phỏng là công cụ có ích phục vụ cho việc nghiên cứu bản chất của các quỹ đạo này, như đề cập đến trong Holman, M.J. and P.A. Wiegert, 1999, “Long-Term Stability of Planets in Binary Systems,” Astronomical Journal 117, có thể tải về từ http://citeseer.ist.psu.edu/358720.html.

Hãy đọc bài báo này rồi sửa lại chương trình mô phỏng hành tinh đã viết để lặp lại hoặc mở rộng kết quả.


  1. Xem http://en.wikipedia.org/wiki/Vector_(spatial).
  2. Độ lớn đôi khi cũng được gọi là độ dài (“length”) nhưng tôi tránh dùng từ này vì sẽ lẫn với hàm length để trả lại số phần tử có trong một véc-tơ của MATLAB.
  3. See http://en.wikipedia.org/wiki/Gravity
  4. Xem http://en.wikipedia.org/wiki/Binary_star.

8 phản hồi

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

8 responses to “Chương 12: Bây giờ véc-tơ mới thật là véc-tơ

  1. CHÚ Ý: Đường link đến bài báo trong Bài tập thứ hai (Hệ hai sao) là:
    http://arxiv.org/pdf/astro-ph/9809315

  2. ===== HẾT =====
    Chương 12 này đã kết thúc cuốn sách: “Mô hình hóa hiện tượng vật lý bằng MATLAB” (Physical Modeling with MATLAB) của tác giả Allen B. Downey. Tôi rất vui khi dịch cuốn sách này; hi vọng nó sẽ giúp bạn ít nhiều, không chỉ về lập trình mà còn về cách hiểu được cách thiết lập một mô hình toán.
    –NQC

  3. Bạn đọc có thể tải về một bản PDF của cuốn sách, với nội dung gần giống hệt các chương đã được post lên đây, chỉ khác ở chỗ dùng ngôn ngữ Octave (có cú pháp giống MATLAB, nhưng phần mềm là tự do / miễn phí). Hãy bấm vào Menu “MATLAB” rồi xem “Ebook mới …”

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

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

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

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

  8. Huệ

    cám ơn thầy đã dịch quyển sách này, nó rất dễ hiểu và hữu ích cho một du học sinh vừa bắt đầu học matlab như em! Nếu có thể mong thầy sẽ post nhiều bài về mô hình hóa vật lý trong Matlab hơn nữa ạ.
    Trân trọng cám ơn thầy và chúc thầy nhiều sức khỏe

Gửi phản hồ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 Log Out / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Log Out / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Log Out / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Log Out / Thay đổi )

Connecting to %s