Chương 10: Các hệ bậc hai

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

Hàm lồng ghép

Trong Mục ví dụ bài toán con vịt, ta đã thấy một ví dụ của tập tin M với hơn một hàm:

function res = duck()
    error = error_func(10)
end

function res = error_func(h)
    rho = 0.3;      % density in g / cm^3
    r = 10;         % radius in cm
    res = ...
end

Vì hàm thứ nhất kết thúc trước hàm thứ hai, chúng có cùng cấp thụt đầu dòng. Những hàm như vậy được gọi là song song, trái với lồng ghép. Một hàm lồng ghép được định nghĩa ở bên trong hàm khác, kiểu như:

function res = duck()
    error = error_func(10)

    function res = error_func(h)
        rho = 0.3;      % density in g / cm^3
        r = 10;         % radius in cm
        res = ...
    end
end

Hàm cấp cao nhất, duck, là hàm ngoàierror_funchàm trong.

Các hàm lồng ghép rất có ích vì các biến của hàm ngoài có thể truy cập đến từ hàm trong. Điều này không thực hiện được với các hàm song song.

Ở ví dụ này, việc dùng hàm lồng ghép khiến cho ta có thể chuyển các tham số rhor ra ngoài error_func.

function res = duck(rho)
    r = 10;
    error = error_func(10)

    function res = error_func(h)
        res = ...
    end
end

Cả rhor đều có thể truy cập đến từ error_func. Bằng cách biến rho trở thành đối số đầu vào, ta đã khiến việc thử nghiệm duck với các giá trị tham số khác nhau được dễ dàng hơn.

Chuyển động Newton

Định luật II Newton về chuyển động thường được viết như sau

F = ma

trong đó F là hợp lực tác dụng lên một vật, m là khối lượng của vật, còn a là gia tốc mà vật thu được. Chỉ trong trường hợp đơn giản khi vật chuyển động theo đường thẳng thì cả Fa đều là số vô hướng, nhưng nhìn chung chúng là véc-tơ.

Thậm chí khái quát hơn nữa, nếu Fa thay đổi theo thời gian thì chúng là một hàm và kết quả của việc lượng giá F(t) là một véc-tơ mô tả hợp lực tại thời điểm t. Như vậy một cách viết tường minh hơn định luật Newton là

∀ t: \vec{F}(t) = \vec{a}(t)

Thứ tự sắp xếp của phương trình này cho thấy nếu ta đã biết ma thì có thể tính được lực. Điều này là đúng, nhưng ở phần lớn các mô phỏng vật lý thì ngược lại. Dựa trên mô hình vật lý, bạn biết Fm, rồi đi tính a.

Vậy nếu đã biết gia tốc, a, dưới dạng hàm theo thời gian, bằng cách nào có thể tìm được vị trí p của vật? À, ta biết rằng gia tốc là đạo hàm bậc hai của vị trí, vậy nên ta có thể viết một phương trình vi phân

ptt = a

trong đó ap là các hàm số theo thời gian và trả lại các véc-tơ, còn ptt là đạo hàm bậc hai theo thời gian của p.

Vì phương trình này bao gồm một đạo hàm bậc hai nên nó được gọi là PVT bậc hai. ode45 không thể giải phương trình dưới dạng này, nhưng bằng cách đưa vào một biến mới, v, đại diện cho vận tốc thì ta có thể viết lại nó như một hệ các PVT bậc nhất.

pt = v
vt = a

Phương trình đầu phát biểu rằng đạo hàm bậc nhất của pv; phương trình thứ hai cho thấy đạo hàm của va.

Hiện tượng rơi tự do

Ta hãy bắt đầu bằng một ví dụ đơn giản, một vật rơi tự do trong chân không (ở đó không có sức cản của không khí). Gần mặt đất, gia tốc trọng trường là g = –9,8 m/s2, trong đó dấu trừ biểu thị rằng trọng lực hướng xuống dưới.

Nếu vật rơi thẳng xuống (cùng theo hướng trọng lực) thì ta có thể mô tả vị trí của nó bằng một giá trị vô hướng là độ cao. Vì vậy, tạm thời bây giờ bài toán là một chiều.

Sau đây là hàm tốc độ mà ta có thể dùng với ode45 để giải bài toán này:

function res = freefall(t, X)
    p = X(1);      % the first element is position
    v = X(2);      % the second element is velocity

    dpdt = v;                          
    dvdt = acceleration(t, p, v);

    res = [dpdt; dvdt];    % pack the results in a column vector
end

function res = acceleration(t, p, v)
    g = -9.8;      % acceleration of gravity in m/s^2
    res = g;
end

Hàm thứ nhất là hàm tốc độ. Nó nhận vào các biến tX, trong đó các phần tử của X được hiểu là vị trí và vận tốc. Giá trị trả về từ freefall là một véc-tơ (cột) chứa những đạo hàm của vị trí và vận tốc, tức lần lượt là vận tốc và gia tốc.

Việc tính pt thật dễ dàng vì ta đã có vận tốc như một phần tử của X. Thứ duy nhất mà ta cần tính là gia tốc, và đó là việc mà hàm thứ hai đảm nhiệm.

acceleration làm nhiệm vụ tính gia tốc như một hàm của thời gian, vị trí và vận tốc. Ở ví dụ này, gia tốc tổng hợp là một hằng số, vì vậy ta không cần phải kể ra tất cả những thông tin trên, nhưng rồi sắp tới ta sẽ phải nhắc đến.

Sau đây là cách chạy ode45 với hàm tốc độ này:

>> ode45(@freefall, [0, 30], [4000, 0])

Như mọi khi, đối số thứ nhất là chuôi của hàm, thứ hai là khoảng thời gian (30 giây) và thứ ba là điều kiện ban đầu: trong trường hợp này, độ cao ban đầu là 4000 mét và vận tốc ban đầu là 0. Vì vậy bạn có thể hình dung “vật” là một người nhảy dù vừa rời khỏi máy bay ở độ cao 12000 feet.

Và kết quả sẽ như sau:

Đồ thị biểu diễn hiện tượng rơi tự do

Dòng cuối cùng cho thấy vận tốc bắt đầu tại 0 và hạ xuống một cách tuyến tính. Dòng đầu tiên cho thấy vị trí bắt đầu tại 4000 m và hạ xuống theo hàm parabol (nhớ rằng đây là hàm parabol theo thời gian chứ không phải đường bay hình parabol).

Lưu ý rằng ode45 không biết vị trí mặt đất ở đâu, vì vậy người nhảy dù vẫn hạ độ cao từ không xuống các giá trị độ cao âm. Ta sẽ giải quyết vấn đề này sau.

Lực cản không khí

Để giúp cho việc mô phỏng được xác thực hơn, ta có thể thêm vào lực cản không khí. Với những vật lớn di chuyển nhanh trong không trung, lực cản không khí tác dụng lên vật sẽ tỷ lệ thuận với v2:

Fdrag = cv2

Trong đó c là hằng số cản, phụ thuộc vào các yếu tố mật độ không khí, diện tích mặt cắt ngang của vật thể và đặc tính bề mặt của vật thể. Với mục đích minh họa của bài tập này, có thể lấy c = 0,2.

Để chuyển đổi từ lực sang gia tốc, ta cần phải biết khối lượng, vì vậy hãy coi rằng người nhảy dù (cùng với thiết bị đi theo người) nặng 75 kg.

Sau đây là một phiên bản của acceleration có tính đến lực cản không khí (bạn không cần phải thay đổi gì trong freefall):

function res = acceleration(t, p, v)
    a_grav = -9.8;              % acceleration of gravity in m/s^2
    c = 0.2;                    % drag constant
    m = 75;                     % mass in kg
    f_drag = c * v^2;           % drag force in N
    a_drag = f_drag / m;        % drag acceleration in m/s^2
    res = a_grav + a_drag;      % total acceleration
end

Dấu của lực cản (và gia tốc) là dương nếu khi vật đang rơi, hướng của lực cản hướng lên trên. Gia tốc tổng hợp là tổng của gia tốc trọng lực và gia tốc cản. Hãy cẩn thận khi bạn làm việc với các lực và gia tốc; hãy chắc rằng bạn chỉ cộng lực với lực và gia tốc với gia tốc. Trong mã lệnh vừa viết, tôi dùng lời chú thích để tự nhắc nhở mình rằng các đơn vị ứng với từng giá trị. Điều đó sẽ giúp tôi tránh khỏi những điều vô nghĩa kiểu như cộng lực và gia tốc với nhau.

Sau đây là dạng của kết quả khi có lực cản không khí:

Đồ thị biểu thị hiện tượng rơi tự do [2]

Có sự khác biệt lớn. Với lực cản không khí, vận tốc tăng lên đến khi gia tốc cản bằng với g; sau đó, vận tốc là không đổi, được biết đến với tên gọi “vận tốc cuối,” và vị trí sẽ hạ thấp tuyến tính theo thời gian (và sẽ chậm hơn nhiều so với trường hợp rơi trong chân không). Để kiểm tra lại kết quả một cách cẩn thận hơn, ta có thể gán chúng vào các biến

>> [T, M] = ode45(@freefall, [0, 30], [4000, 0]);

rồi đọc các vị trí và vận tốc cuối:

>> M(end,1)

ans = 2.4412e+03          % altitude in meters

>> M(end,2)

ans = -60.6143            % velocity in m/s

Hãy tăng khối lượng của người nhảy dù, và kiểm tra lại để chắc rằng vận tốc cuối tăng lên. Quan hệ này chính là nguồn gốc dẫn đến cảm nhận trực giác rằng những vật nặng thì rơi nhanh hơn; và đúng là trong không khí, chúng rơi nhanh hơn thật!

Nhảy dù!

Ở mục trước, ta đã thấy rằng vận tốc cuối của một người nhảy dù nặng 75 kg là khoảng 60 m/s,vốn là khoảng 130 mph (dặm/giờ). Nếu quả thật phải tiếp đất với vận tốc như vậy, chắc là bạn sẽ chết. Đó là tại sao phải có dù.

Hãy sửa lại acceleration sao cho sau 30 giây tự do, người nhảy mở dù, và (gần như) lập tức tăng hằng số cản lên 2,7.

Vận tốc cuối bây giờ bằng bao nhiêu? Bao lâu (sau khi mở dù), người nhảy sẽ tiếp đất?

Bài toán hai chiều

Đến đây ta đã dùng ode45 giải một hệ phương trình bậc nhất và một phương trình bậc hai. Bước làm hợp lý tiếp theo là một hệ phương trình bậc hai, và ví dụ hợp lý tiếp theo là về một đường bay. “Đường bay” được vạch ra bởi một vật thể bắn đi trong không khí, thường là hướng về nhằm phá hủy mục tiêu.

Nếu đường bay nằm trong một mặt phẳng, ta có thể coi hệ là hai chiều, với x biểu thị khoảng cách theo phương ngang bay được và y biểu thị độ cao. Như vậy bây giờ thay vì một người nhảy dù, bạn có thể hình dung về người trong gánh xiếc bị bắn ra từ nòng pháo.

Theo Wikipedia1, khoảng cách kỉ lục của người bị bắn từ nòng pháo là xa 56,5 mét.

Sau đây là một khung chương trình có dùng ode45 để tính đường bay của “đạn” trong không gian hai chiều:

function res = projectile(t, W)
    P = W(1:2);
    V = W(3:4);

    dPdt = V;                          
    dVdt = acceleration(t, P, V);

    res = [dPdt; dVdt];
end

function res = acceleration(t, P, V)
    g = -9.8;             % acceleration of gravity in m/s^2
    res = [0; g];
end

Đối số thứ hai của hàm tốc độ là một véc-tơ, W, với 4 phần tử. Hai phần tử đầu được gán cho P, để biểu diễn vị trí; hai phần sau gán cho V, để biểu diễn vận tốc. PV là các véc-tơ với các phần tử là những thành phần theo phương xy.

Kết quả thu được từ acceleration cũng là một véc-tơ; (tạm thời) bở qua lực cản không khí, gia tốc theo phương x bằng không, và theo phương y thì bằng g. Ngoài ra thì mã lệnh giống hệt như ta đã thấy ở Mục [Rơi tự do].

Nếu ta phóng người bay từ độ cao ban đầu là 3 mét, với các vận tốc thành phần là 40 m/s và 30 m/s theo các phương xy thì lời gọi hàm ode45 sẽ như sau:

ode45(@projectile, [0,10], [0, 3, 40, 30]);

Và kết quả sẽ như sau:

Đường bay của vật thể

Bạn có thể sẽ phải nghĩ một lát để hình dung ra đường nào ứng với đại lượng gì. Dường như là thời gian bay vào khoảng 6 giây.

Hãy lọc ra các giá trị thành phần theo các phương xy của vị trí, vạch ra đường bay của người, rồi ước tính khoảng cách bay được.

Hãy thêm sức cản không khí vào bài toán mô phỏng này. Ở trường hợp người nhảy dù, ta giả sử rằng hằng số cản bằng 0,2; nhưng nó được dựa trên giả thiết là người nhảy rơi có mặt tiếp xúc phẳng. Nhưng người bị bắn pháo, bay hướng đầu về phía trước, có lẽ sẽ ứng với hằng số cản chừng 0,1. Vận tốc ban đầu cần để đạt được khoảng cách bay kỷ lục 65,6 sẽ phải bằng bao nhiêu? Gợi ý: góc bắn tối ưu bằng bao nhiêu?

Điều gì trục trặc có thể xảy ra?

Điều gì trục trặc có thể xảy ra? À, chẳng hạn như vertcat. Để giải thích từ này nghĩa là gì, tôi sẽ bắt đầu bằng phép ghép nối, phép toán để ghép hai ma trận thành một ma trận lớn hơn. “Vertical catenation” (ghép theo phương thẳng đứng) sẽ ghép hai ma trận bằng cách đặt chồng lên nhau (cột nối cột); “horizontal catenation” (ghép theo phương ngang) đặt hai ma trận cạnh nhau (hàng nối hàng).

Sau đây là một ví dụ ghép ngang với các véc-tơ hàng:

>> x = 1:3

x = 1     2     3

>> y = 4:5

y = 4     5

>> z = [x, y]

z = 1     2     3     4     5

Bên trong cặp ngoặc vuông, dấu phẩy làm nhiệm vụ ghép ngang. Còn toán tử ghép dọc là dấu chấm phẩy. Sau đây là một ví dụ với các ma trận:

>> X = zeros(2,3)

X =  0     0     0
     0     0     0

>> Y = ones(2,3)

Y =  1     1     1
     1     1     1

>> Z = [X; Y]

Z =  0     0     0
     0     0     0
     1     1     1
     1     1     1

Các toán tử này chỉ có tác dụng khi các ma trận có cùng kích thước ở chiều được ghép với nhau. Nếu không, bạn sẽ nhận được:

>> a = 1:3

a = 1     2     3

>> b = a'

b =  1
     2
     3

>> c = [a, b]
??? Error using ==> horzcat
All matrices on a row in the bracketed expression must have the 
 same number of rows.

>> c = [a; b]
??? Error using ==> vertcat
All rows in the bracketed expression must have the same 
number of columns.

Ở ví dụ này, a là véc-tơ hàng còn b là véc-tơ cột, vì vậy chúng chẳng thể ghép được theo bất cứ chiều nào.

Bằng việc đọc lời thông báo lỗi, bạn có thể đã đoán ra được rằng horzcat là hàm thực hiện ghép nối theo phương ngang, và tương tự với vertcat cùng việc ghép nối theo phương dọc.

Các phép toán này có liên quan đến projectile vì dòng cuối cùng, vốn để gói dPdtdVdt vào biến đầu ra:

function res = projectile(t, W)
    P = W(1:2);
    V = W(3:4);

    dPdt = V;                          
    dVdt = acceleration(t, P, V);

    res = [dPdt; dVdt];
end

Miễn là cả dPdtdVdt đều là véc-tơ cột thì dấu chấm phẩy sẽ thực hiện ghép nối theo chiều dọc, va kết quả là véc-tơ cột có 4 phần tử. Nhưng nếu một trong hai véc-tơ đó là véc-tơ hàng thì sẽ trục trặc.

ode45 trông đợi kết quả của projectile như một véc-tơ cột, vì vậy nếu bạn dùng ode45, có thể sẽ tốt hơn nếu để mọi đại lượng dưới dạng véc-tơ cột.

Nhìn chung, nếu bạn gặp trục trặc với horzcatvertcat, hãy dùng size để hiển thị kích cỡ của các toán hạng, và đảm bảo chắc chắn véc-tơ của bạn được sắp xếp hàng-cột cho đúng.

Thuật ngữ

các hàm song song:
Những hàm được định nghĩa cạnh nhau, và một hàm kết thúc trước khi hàm khác bắt đầu.

hàm lồng ghép:
Hàm được định nghĩa bên trong hàm khác.

hàm ngoài:
Hàm có chứa một định nghĩa hàm khác.

hàm trong:
Hàm được định nghĩa bên trong một lời định nghĩa hàm khác. Hàm trong có thể truy cập các biến của hàm ngoài.

ghép nối:
Phép toán nối tiếp hai ma trận để tạo nên một ma trận mới.

Bài tập

Đường bay của quả bóng chày bị chi phối bởi ba lực: trọng lực, lực cản không khí, và lực Magnus do chuyển động quay. Nếu ta bỏ qua gió và lự cản Magnus thì đường bay của quả bóng chày luôn nằm trong một mặt phẳng, và ta có thể mô phỏng nó như vật được bắn đi theo không gian hai chiều.

Một mô hình đơn giản cho lực cản không khí tác dụng lên quả bóng là: A simple model of the drag of a baseball is:
Fd = –½ ρv2ACd\hat{V}
trong đó Fd là một véc-tơ biểu diễn lực tác dụng lên quả bóng chày gây bởi lực cản, Cd là hệ số cả (0,3 là một giá trị phù hợp), ρ là mật độ không khí (1,3 kg/m3 ở mực nước biển), A là diện tích mặt cắt ngang của quả bóng (0,0042 m2), v là độ lớn của véc-tơ vận tốc, còn \hat{V} là một véc-tơ đơn vị theo hướng của véc-tơ vận tốc. Khối lượng của quả bóng là 0,145 kg.

Bạn có thể đọc thêm về lực cản ở trang http://en.wikipedia.org/wiki/Drag_(physics).

  • Hãy viết một hàm nhận vào các biến gồm vận tốc ban đầu và góc ném của quả bóng chày, dùng ode45 để tính đường bay, và trả lại tầm xa (khoảng cách phương ngang bay được) là các biến đầu ra.
  • Hãy viết một hàm nhận vào biến vận tốc ban đầu của quả bóng, tính góc ném sao cho khiến cho tầm xa là lớn nhất, và trả lại góc ném tối ưu nêu trên cùng với tầm xa như những biến đầu ra. Góc ném tối ưu thay đổi theo vận tốc ban đầu như thế nào?
  • Khi đội Red Sox đoạt giải World Series năm 2007, họ gặp đội Colorado Rockies tại sân nhà ở Denver, Colorado. Hãy ước tính mật độ không khí ở Mile High City. Nó có ảnh hưởng gì đến sức cản không khí? Hãy dự đoán về ảnh hưởng của yếu tố này [lực cản] đến góc ném tối ưu, và dùng mô phỏng của bạn để kiểm tra lại dự đoán đó.
  • Sân Green Monster ở Fenway Park ở độ cao khoảng 12 m và chiều dài khoảng 97 m. Vận tốc tối thiểu của quả bóng ngay sau khi vụt phải bằng bao nhiêu để bay hết chiều dài sân (giả thiết rằng bóng bay theo góc tối ưu)? Bạn có nghĩ rằng một người có thể ném bóng hết chiều dài sân được không?
  • Lực cản của không khí lên quả bóng chày thực tế còn phức tạp hơn rất nhiều so với tính toán bằng mô hình đơn giản nêu trên. Đặc biệt, hệ số cản còn thay đổi theo vận tốc. Bạn có thể đọc thêm từ cuốn The Physics of Baseball2; và cũng tìm thông tin trên web. Sau đó, hãy thiết lập một mô hình thực tế hơn về lực cản và sửa đổi chương trình tính theo lực cản mới này. Bằng mô hình mới thì tầm xa tính được sẽ thay đổi bao nhiêu? Góc ném tối ưu sẽ thay đổi bao nhiêu?

  1. http://en.wikipedia.org/wiki/Human_cannonball
  2. Robert K. Adair, Harper Paperbacks, 3rd Edition, 2002.
Advertisements

4 phản hồi

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

4 responses to “Chương 10: Các hệ bậc hai

  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