MATLAB đã sử dụng thuật toán nào de giải phương trình vi phân bằng phương pháp hình thức

Home - Video - Giải phương trình vi phân trong MATLAB

Prev Article Next Article

source

Xem ngay video Giải phương trình vi phân trong MATLAB

“Giải phương trình vi phân trong MATLAB “, được lấy từ nguồn: https://www.youtube.com/watch?v=YuLwnon_1jU

Tags của Giải phương trình vi phân trong MATLAB: #Giải #phương #trình #phân #trong #MATLAB

Bài viết Giải phương trình vi phân trong MATLAB có nội dung như sau:

Từ khóa của Giải phương trình vi phân trong MATLAB: vi phân

Thông tin khác của Giải phương trình vi phân trong MATLAB:
Video này hiện tại có lượt view, ngày tạo video là 2021-11-08 02:34:58 , bạn muốn tải video này có thể truy cập đường link sau: https://www.youtubepp.com/watch?v=YuLwnon_1jU , thẻ tag: #Giải #phương #trình #phân #trong #MATLAB

Cảm ơn bạn đã xem video: Giải phương trình vi phân trong MATLAB.

Prev Article Next Article

Home - Video - hướng dẫn giải phương trình vi phân bằng matlab – matlab tutorial 06

Prev Article Next Article

source

Xem ngay video hướng dẫn giải phương trình vi phân bằng matlab – matlab tutorial 06

“hướng dẫn giải phương trình vi phân bằng matlab – matlab tutorial 06 “, được lấy từ nguồn: https://www.youtube.com/watch?v=z3ahCL3jjdc

Tags của hướng dẫn giải phương trình vi phân bằng matlab – matlab tutorial 06: #hướng #dẫn #giải #phương #trình #phân #bằng #matlab #matlab #tutorial

Bài viết hướng dẫn giải phương trình vi phân bằng matlab – matlab tutorial 06 có nội dung như sau:

Từ khóa của hướng dẫn giải phương trình vi phân bằng matlab – matlab tutorial 06: vi phân

Thông tin khác của hướng dẫn giải phương trình vi phân bằng matlab – matlab tutorial 06:
Video này hiện tại có lượt view, ngày tạo video là 2022-01-15 23:02:59 , bạn muốn tải video này có thể truy cập đường link sau: https://www.youtubepp.com/watch?v=z3ahCL3jjdc , thẻ tag: #hướng #dẫn #giải #phương #trình #phân #bằng #matlab #matlab #tutorial

Cảm ơn bạn đã xem video: hướng dẫn giải phương trình vi phân bằng matlab – matlab tutorial 06.

Prev Article Next Article

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

Phương trình vi phân

Phương trình vi phân là phương trình mô tả các đạo hàm của một hàm số chưa biết. “Giải phương trình vi phân” nghĩa là tìm một hàm số có các đạo hàm thỏa mãn phương trình đã cho.

Chẳng hạn, khi vi phuẩn sống trong môi trường đặc biệt thuận lợi thì tốc độ sinh trưởng tại bất kì thời điểm nào cũng tỉ lệ thuận với số vi khuẩn lúc đó. Có lẽ điều ta quan tâm là số vi khuẩn được biểu diễn dưới dạng hàm theo thời gian. Ta hãy định nghĩa f là hàm chiếu từ thời gian, t, đến số vi khuẩn, y. Dù không biết y bằng bao nhiêu, nhưng ta vẫn có thể viết một phương trình để mô tả nó:

df / dt = af

trong đó a là hằng số đặc trưng cho mức độ vi khuẩn tăng nhanh bao nhiêu.

Lưu ý rằng cả hai vế của phương trình đều là hàm số. Khi ta nói hai hàm số bằng nhau nghĩa là các giá trị của chúng luôn luôn bằng nhau. Nói cách khác:

∀ t: (df / dt)(t) = af(t)

Đây là một phương trình vi phân thường (PVT) vì tất cả các đạo hàm đều được lấy theo cùng một biến. Nếu như phương trình liên hệ các đạo hàm theo nhiều biến khác nhau (đạo hàm từng phần), thì ta sẽ có phương trình đạo hàm riêng.

Phương trình này là bậc nhất vì nó chỉ có chứa những đạo hàm cấp một. Nếu có mặt đạo hàm bậc hai, phương trình sẽ là bậc hai, và cứ như vậy.

Phương trình này là tuyến tính vì mỗi số hạng chứa t, f hoặc df / dt đều chỉ với bậc lũy thừa bằng một; nếu bất kì số hạng nào có chứa tích các đại lượng trên hoặc các lũy thừa của t, f và df / dt thì phương trình sẽ là phi tuyến.

Các PVT bậc nhất, tuyến tính đều có thể giải được theo cách giải tích; nghĩa là ta có thể biểu diễn nghiệm dưới dạng một hàm số của t. Riêng PVT này có vô số nghiệm, nhưng tất cả nghiệm đó đều có chung dạng sau:

f(t) = beat

Với giá trị b bất kì, hàm số này đều thỏa mãn PVT. Nếu bạn không tin điều này, hãy tự lấy đạo hàm và kiểm tra.

Nếu ta biết số vi khuẩn tại một thời điểm nhất định thì ta có thể dùng thông tin nói trên để xác định trong số vô hạn các nghiệm, đâu là nghiệm duy nhất thỏa mãn điều kiện trên để mô tả diễn biến của số vi khuẩn thực tế theo thời gian.

Chẳng hạn, nếu đã biết rằng f(0) = 5 tỷ tế bào thì ta có thể viết:

f(0) = 5 = bea × 0

và giải tìm ra b bằng 5. Từ đó ta tìm được hàm mong muốn:

f(t) = 5eat

Thông tin thêm nói trên, vốn quyết định giá trị của b, được gọi là điều kiện ban đầu (cho dù không phải lúc nào nó cũng được chỉ định tại t = 0).

Không may là phần lớn các hệ vật lý thú vị đều được mô tả bởi các phương trình vi phân phi tuyến, mà đa số đều không thể giải được theo cách giải tích. Một cách khác là giải bằng phương pháp số.

Phương pháp Euler

Phương pháp số đơn giản nhất để giải PVT là phương pháp Euler. Sau đây là một bài kiểm tra nhỏ để xem bạn có thông minh được như Euler không. Giả sử như ở thời điểm t bạn đo được số vi khuẩn, y, và tốc độ tăng, r. Theo bạn thì số vi khuẩn sẽ là bao nhiêu sau một thời gian Δ t trôi qua?

Nếu bạn nói rằng y + rΔ t thì xin chúc mừng! Bạn vừa phát minh ra phương pháp Euler (nhưng bạn vẫn chưa giỏi bằng Euler).

Cách ước tính này dựa trên giả sử rằng r là hằng số, nhưng nhình chung thì không phải vậy, do đó ta chỉ trông đợi là ước lượng sẽ tốt nếu r thay đổi từ từ và Δ t nhỏ.

Nhưng hãy (tạm) giả sử rằng PVT đang xét có thể được viết sao cho

(df / dt)(t) = g(t, y)

trong đó g là một hàm nào đó chiếu từ (t, y) đến r; nghĩa là, cho trước thời điểm và số vi khuẩn, hàm này sẽ tính tốc độ thay đổi. Sau đó ta có thể tiến từ một thời điểm đến thời điểm tiếp theo bằng các phương trình sau đây:

Tn+1 = Tn + Δt

Fn+1 = Fn + g(t,y) Δt

Ở đây {Ti} là một dãy các thời điểm mà tại đó ta cần tính các giá trị của f, còn {Fi} là một dãy các giá trị ước tính được. Với mỗi chỉ số i, Fi là một ước lượng của f(Ti). Khoảng Δ t được gọi là bước thời gian.

Giả sử rằng ta khởi đầu tại t = 0 và có một điều kiện ban đầu f(0) = y0 (trong đó y0 là một giá trị cụ thể đã biết); đặt T1 = 0 và F1 = y0, sau đó dùng các PT {euler1} và {euler2} để tính các giá trị của Ti và Fi đến khi Ti nhận giá trị của t mà ta quan tâm.

Nếu tốc độ tăng thay đổi không quá nhanh và bước thời gian không quá dài thì phương pháp Euler cũng đủ chính xác với nhiều mục đích tính toán. Một cách kiểm tra là chạy nó lần đầu với bước thời gian Δ t và lần sau với bước thời gian Δt / 2. Nếu các kết quả như nhau, thì có lẽ kết quả là đúng; nếu không, hãy thử rút ngắn bước thời gian một lần nữa.

Phương pháp Euler là bậc nhất, nghĩa là mỗi lần bạn chia đôi bước tính, bạn trông đợi sai số của giá trị ước tính sẽ giảm bớt một nửa. Với một phương pháp bậc hai, bạn trông đợi sai số giảm đi 4 lần; với phương pháp bậc ba giảm 8 lần, v.v. Cái giá phải trả cho phương pháp bậc cao là chúng phải lượng giá g nhiều lần hơn trong cùng một bước thời gian.

Lưu ý thêm về cách viết

Có rất nhiều kí hiệu toán học trong chương này, vì vậy tôi sẽ tạm dừng ở đây để ôn lại những gì chúng ta đã học được đến giờ. Dưới đây là các biến số, ý nghĩa, và kiểu của chúng:

t	thời gian	biến vô hướng
Δ t	bước thời gian	hằng vô hướng
y	số cá thể	biến vô hướng
r	tốc độ thay đổi	biến vô hướng
f	Hàm chưa biết được chỉ định	hàm t → y
	dưới hình thức ẩn trong PVT.
df / dt	Đạo hàm bậc nhất theo thời gian của f	hàm t → r
g	“Hàm tốc độ,” được suy ra từ
	PVT, để tính tốc độ thay đổi với	hàm t, y → r
	giá trị bất kì của t, y.
T	một loạt các thời điểm, t, tại đó	dãy
	ta ước tính f(t)
F	một loạt các giá trị ước tính của f(t)	dãy

Như vậy f là một hàm để tính số cá thể (vi trùng), có dạng một hàm theo thời gian, f(t) là giá trị của hàm được tính ở một thời điểm nhất định, và nếu ta gán f(t) cho một biến, ta thường gọi biến đó là y.

Tương tự, g là một “hàm tốc độ” để tính tốc độ thay đổi dưới hình thức một hàm phụ thuộc vào thời gian và số cá thể. Nếu ta gán g(t, y) cho một biến, ta thường gọi nó là r.

df / dt là đạo hàm bậc nhất của f, nó chiếu từ t đến một giá trị tốc độ. Nếu ta gán df / dt(t) cho một biến, ta gọi nó là r.

Thường thì df / dt rất dễ bị lẫn với g, nhưng lưu ý rằng cả về kiểu chúng đã khác nhau rồi. g có tính khái quát hơn: nó có thể tính được tốc độ thay đổi của bất kì số cá thể (giả định) nào tại thời điểm bất kì; df / dt thì cụ thể hơn: đó là tốc độ thay đổi thực sự tại thời điểm t, khi biết số cá thể là f(t).

ode45

Một hạn chế của phương pháp Euler là bước thời gian không đổi qua các lượt lặp khác nhau. Nhưng có những phần của lời giải lại khó ước tính hơn phần khác; nếu bước thời gian đủ ngắn để giải được phần khó thì cũng bước thời gian này sẽ khiến việc tính toán trở nên quá nhiều đối với những phần dễ. Cách làm lý tưởng nhất là điều chỉnh bước thời gian trong khi giải. Các phương pháp như vậy được gọi là thích ứng, và một trong những phương pháp thích ứng tốt nhất là phương pháp Dormand-Prince dùng cặp công thức Runge-Kutta. Bạn chưa cần biết chi tiết về công thức này, vì những người phát triển MATLAB đã đưa nó vào một hàm có tên là ode45. Trong đó chữ ode có nghĩa là “ordinary differential equation [solver];” ([chương trình giải] PVT) còn số 45 là để chỉ rằng đã dùng cách kết hợp các công thức bậc 4 và bậc 5.

Để dùng ode45, bạn phải viết một hàm MATLAB để ước tính g như một hàm của t và y.

Sau đây là một ví dụ minh họa. Giả sử rằng tốc độ tăng trưởng của chuột phụ thuộc vào số con chuột hiện thời và mức độ sẵn có của thức ăn, vốn thay đổi theo thời gian trong năm. Phương trình cơ bản sẽ có dạng

(df / dt)(t) = af(t)[1 + sin(ωt)]
trong đó t là thời gian tính theo ngày và f(t) là số con chuột tại thời điểm t, còn a và ω là các tham số. Tham số là một giá trị đặc trưng cho một khía cạnh vật lý của hệ được mô phỏng. Chẳng hạn, ở Bài tập {duck} ta đã dùng các tham số rho và r để định lượng khối lượng riêng và bán kính của một con vịt. Các tham số thường không đổi, nhưng cũng có thể thay đổi theo thời gian trong một số mô hình.

Ở ví dụ này, a đặc trưng cho tốc độ sinh sản, còn ω là tần số của một hàm tuần hoàn để mô tả ảnh hưởng của mức độ cung cấp thức ăn đến sự sinh sản.

Phương trình này đặc trưng cho quan hệ giữa một hàm với đạo hàm của nó. Để ước tính được các giá trị của f bằng cách số trị, ta phải chuyển nó về một hàm tốc độ.

Bước đầu tiên là giới thiệu một biến, y, là một tên gọi khác của f(t)

(df / dt)(t) = ay[1 + sin(ωt)]

Phương trình này có nghĩa là nếu cho trước các giá trị t và y, ta có thể tính df / dt(t), vốn là tốc độ thay đổi của f. Bước tiếp theo là biểu diễn phép tính đó đưới dạng một hàm gọi là g:

g(t, y) = ay[1 + sin(ωt)]

Cách viết này rất có ích vì ta có thể dùng hàm với phương pháp Euler hoặc ode45 để ước tính các giá trị của f. Tất cả những việc ta cần làm chỉ là viết một hàm MATLAB để lượng giá g. Sau đây là mã lệnh ứng với các giá trị a = 0. 01 và ω = 2π / 365 (một chu kì mỗi năm):

function res = rats(t, y) a = 0.01; omega = 2 * pi / 365; res = a * y * (1 + sin(omega * t)); end

Bạn có thể kiểm tra hàm này từ Command Window bằng cách gọi nó với các giá trị khác nhau của t và y; kết quả sẽ là tốc độ thay đổi (đơn vị là số chuột mỗi ngày):

>> r = rats(0, 2) r = 0.0200

Như vậy nếu có 2 con chuột vào ngày 1/1, ta sẽ trông đợi chúng đẻ thêm với một tốc độ sao cho cứ 2 con chuột được sinh ra trong vòng 100 ngày. Nhưng đến tháng 4, tốc độ này đã tăng gấp đôi:

>> r = rats(120, 2) r = 0.0376

Vì tốc độ tăng thì luôn thay đổi, nên không dễ dự đoán số con chuột trong tương lai; song đó chính là việc mà ode45 đảm nhiệm. Sau đây là cách dùng nó:

>> ode45(@rats, [0, 365], 2)

Đối số đầu tiên là chuôi của hàm để tính g. Đối số thứ hai là khoảng thời gian mà ta quan tâm, ở đây là 1 năm. Đối số thứ ba là số con chuột ban đầu, f(0) = 2.

Khi bạn gọi ode45 mà không gán kết quả cho một biến, MATLAB sẽ biểu thị kết quả dưới dạng hình vẽ:

Trục x chỉ thời gian từ 0 đến 365 ngày; trục y chỉ số con chuột, bắt đầu là 2 rồi tăng dần đến gần 80. Tốc độ tăng khá chậm vào mùa đông và mùa hè, và nhanh hơn vào mùa xuân và mùa thu, đồng thời cũng tăng nhanh hơn khi số chuột nhiều lên.

Nhiều biến đầu ra

ode45 chỉ là một trong số nhiều hàm MATLAB trả lại kết quả dưới dạng nhiều biến. Cú pháp của lời gọi nó và lưu giữ kết quả là

>> [T, Y] = ode45(@rats, [0, 365], 2);

Giá trị trả lại thứ nhất được gán cho T; giá trị thứ hai được gán cho Y. Mỗi phần tử của T là một thời điểm t, mà tại đó ode45 ước tính số chuột; mỗi phần tử của Y là một giá trị ước đoán của f(t).

Nếu bạn gán các giá trị đầu ra cho biến, ode45 sẽ không vẽ hình; phạn phải tự làm điều đó:

>> plot(T, Y, 'bo-')

Nếu vẽ đồ thị các phần tử của T, bạn sẽ thấy rằng khoảng cách giữa các điểm rất không bằng nhau. Ở đầu khoảng thời gian các điểm sát nhau và về sau này thì tách rời ra.

Để biết được số chuột lúc cuối năm, bạn có thể hiển thị phần tử cuối của mỗi véc-tơ:

>> [T(end), Y(end)] ans = 365.0000 76.9530

end là một từ đặc biệt trong MATLAB; khi xuất hiện ở vị trí một chỉ số, nó có nghĩa là “chỉ số của phần tử cuối cùng.” Bạn có thể dùng nó trong một biểu thức, vì vậy Y(end-1) là phần tử áp chót của Y.

Số con chuột lúc cuối năm sẽ thay đổi bao nhiêu nếu bạn tăng gấp đôi số chuột ban đầu? Nó sẽ thay đổi bao nhiêu nếu bạn tính đến thời gian sau hai năm? Và nếu bạn tăng a gấp đôi?

Giải tích hay số trị?

Khi bạn giải PVT theo cách giải tích, kết quả là một hàm, f, cho phép ta tính số cá thể, f(t), với bất kì giá trị nào của t. Khi bạn giải PVT theo cách số trị, bạn nhận được hai véc-tơ. Bạn có thể hình dung hai véc-tơ này như xấp xỉ của hàm f có tính liên tục: “rời rạc” là bởi vì nó chỉ xác định với những giá trị cụ thể của t, và “xấp xỉ” là vì mỗi giá trị Fi chỉ là một ước đoán cho giá trị thật f(t).

Trên đây là những hạn chế của nghiệm số trị. Còn ưu điểm chính là bạn có thể tính nghiệm số trị cho cả những PVT mà không có nghiệm giải tích, vốn chiếm phần lớn các PVT phi tuyến.

Nếu bạn còn tò mò về cơ chế hoạt động của ode45, bạn có thể sửa rats để hiển thị các điểm, (t, y), tại đó ode45 được dùng để tính g. Sau đây là một phiên bản đơn giản:

function res = rats(t, y) plot(t, y, 'bo') a = 0.01; omega = 2 * pi / 365; res = a * y * (1 + sin(omega * t)); end

Mỗi lần rats được gọi, nó chấm một điểm số liệu; để thấy được tất cả các điểm, bạn phải dùng hold on.

>> clf; hold on >> [T, Y] = ode45(@rats, [0, 10], 2);

Hình vẽ này cho thấy một phần của kết quả đầu ra, được phóng to trong khoảng những ngày từ 100 đến 170:

Các vòng tròn cho thấy những điểm tại đó ode45 gọi đến rats. Các đường thẳng đi qua vòng tròn cho thấy độ dốc (tốc độ thay đổi) tính được ở mỗi điểm. Hình chữ nhật cho thấy các vị trí được ước tính (Ti, Fi). Lưu ý rằng ode45 thường lượng giá g vài lần trong mỗi lần ước tính. Bằng cách này, không những giá trị ước tính được cải thiện mà còn phát hiện được những chỗ mà sai số đang tăng, tại đó cần rút ngắn bước thời gian (hoặc ngược lại).

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

Đừng quên dấu @ ở chuôi của hàm. Nếu quên không viết nó, MATLAB sẽ coi đối số thứ nhất như một lời gọi hàm, và sẽ gọi rats mà không cấp cho đối số nào.

>> ode45(rats, [0,365], 2) ??? Input argument "y" is undefined. Error in ==> rats at 4 res = a * y * (1 + sin(omega * t));

Một lần nữa, thông báo lỗi có vẻ rắc rối, bởi vì có vẻ như vấn đề nảy sinh ở rats. Tôi đã bảo bạn cẩn thận rồi nhé!

Mặt khác, cũng cần nhớ rằng hàm bạn viết ra sẽ được gọi bởi ode45, nghĩa là nó phải có dấu của hàm mà ode45 trông đợi: cụ thể là nhận vào hai biến, t và y, theo đúng thứ tự đó, và trả lại một biến đầu ra, r.

Nếu bạn làm việc với một hàm tốc độ kiểu như:

g(t, y) = ay

Thì bạn có thể đã muốn viết như sau:

function res = rate_func(y) % SAI a = 0.1 res = a * y end

Nhưng cách làm này sai. Tại sao? Vì khi ode45 gọi đến rate_func, nó cấp cho hai đối số. Nếu bạn chỉ lấy có một biến đầu vào, bạn sẽ nhận được lỗi. Vì vậy bạn phải viết một hàm nhận cả t làm biến đầu vào, ngay cả khi bạn không dùng đến nó.

function res = rate_func(t, y) % DU'NG a = 0.1 res = a * y end

Một lỗi thường gặp khác là việc viết một hàm mà không gán kết quả cho một biến đầu ra. Nếu bạn viết lệnh kiểu như:

function res = rats(t, y) a = 0.01; omega = 2 * pi / 365; r = a * y * (1 + sin(omega * t)) % SAI end

rồi gọi nó từ ode45, bạn sẽ nhận được

>> ode45(@rats, [0,365], 2) ??? Output argument "res" (and maybe others) not assigned during call to "/home/downey/rats.m (rats)". Error in ==> rats at 2 a = 0.01; Error in ==> funfun/private/odearguments at 110 f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0. Error in ==> ode45 at 173 [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Thông báo này có vẻ đáng sợ, nhưng nếu bạn đọc dòng thứ nhất và quên đi các dòng sau, bạn sẽ nắm được vấn đề. (Đối số đầu ra “res” (và có thể các đối số khác) không được gán trong quá trình gọi [tới tập tin và hàm].)

Một lỗi khác mà có thể mắc phải khi dùng ode45 là quên mất cặp ngoặc vuông ở đối số thứ hai. Trong trường hợp đó, MATLAB sẽ nghĩ rằng có bốn đối số, và bạn nhận được

>> ode45(@rats, 0, 365, 2) ??? Error using ==> funfun/private/odearguments When the first argument to ode45 is a function handle, the tspan argument must have at least two elements. Error in ==> ode45 at 173 [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Một lần nữa, nếu đọc dòng thứ nhất, bạn có thể hình dung ra được vấn đề (tspan là viết tắt của “time span”; nó tương ứng với khoảng thời gian mà ta chỉ định khi gọi hàm).

Đặc tính cứng

Còn một vấn đề khác mà bạn có thể đối mặt, nhưng nếu nó làm bạn cảm thấy khỏe hơn, thì có lẽ không phải lỗi của bạn: bài toán bạn đang giải có thể thuộc loại cứng1.

Tôi sẽ không trình bày về khía cạnh kĩ thuật liên quan đến độ cứng ở đây, ngoại trừ phải nói rằng với một số bài toán (ở những khoảng thời gian nào đó và điều kiện đầu nào đó) thì bước thời gian cần thiết để kiểm soát sai số là rất ngắn; điều đó nghĩa là thời gian tính toán sẽ rất lâu. Sau đây là một ví dụ:

df / dt = f2 - f3

Nếu bạn giải PVT với điều kiện ban đầu f(0) = δ trên khoảng từ 0 đến 2 / δ, với δ = 0,01, bạn sẽ thấy kết quả sau:

Sau bước chuyển dần từ 0 tới 1, bước thời gian trở nên rất ngắn và việc tính toán rất chậm. Với những giá trị nhỏ hơn của δ, tình hình còn tệ hơn.

Trong trường hợp này, vấn đề khá dễ khắc phục: thay vì ode45 bạn có thể dùng ode23s, một hàm giải PVT với xu hướng hoạt động tốt với bài toán có tính cứng (stiff) (cũng vì vậy mà có chữ “s” trong tên gọi của hàm).

Nói chung, nếu thấy hàm ode45 tính rất lâu thì bạn có thể thử chuyển sang dùng một trong số các hàm giải cho hệ cứng. Tôi không thường xuyên gặp vấn đề như vậy, nhưng nếu bài toán có tính cứng thì sự cải thiện về thời gian giải cỏ vẻ là đáng kể.

Hãy viết một hàm tốc độ cho PVT này rồi dùng ode45 để giải với điều kiện ban đầu và khoảng tính cho trước. Bắt đầu với δ = 0. 1 và giảm dần nó theo 10 lần. Nếu quá mệt mỏi chờ đợi phép tính thực hiện chưa xong, bạn có thể ấn nút Stop ở cửa sổ Figure hoặc ấn Control-C từ trong Command Window.

Bây giờ hãy thay thế ode45 bằng ode23s và thử lại!

Thuật ngữ

phương trình vi phân: Phương trình liên hệ các đạo hàm của một hàm số chưa biết. phương trình vi phân thường: Phương trình vi phân thường mà tất cả đạo hàm được lấy theo cùng một biến. phương trình vi phân riêng: Phương trình vi phân có bao gồm các đạo hàm lấy theo nhiều biến. bậc nhất (PVT): Phương trình vi phân chỉ chứa các đạo hàm bậc nhất. tuyến tính: Phương trình vi phân không có chứa tích hoặc lũy thừa của hàm cùng các đạo hàm của nó. bước thời gian: Khoảng thời gian giữa hai lần ước tính kế tiếp trong lời giải số trị của một phương trình vi phân. bậc nhất (phương pháp số): Phương pháp theo đó sai số xấp xỉ được trông đợi sẽ giảm đi một nửa khi bước tính được rút ngắn còn một nửa. thích ứng: Phương thức có thể điều chỉnh được bước thời gian để kiểm soát sai số. độ cứng: Đặc trưng của PVT khiến cho một số hàm giải PVT chạy rất chậm (hoặc tính ra kết quả không chính xác). Một số hàm giải PVT như ode23s, được thiết kế để giải các bài toán có độ cứng. tham số: Giá trị xuất hiện trong mô hình nhằm định lượng một khía cạnh vật lý nào đó của hệ được mô phỏng.

Bài tập

Giả sử rằng bạn có tách chứa 8 ounce cà-phê nóng 90 ∘ C và 1 ounce kem ở nhiệt độ phòng, 20 ∘ C. Bạn biết qua một kinh nghiệm cay đắng là chỉ uống được cà-phê nóng nhất ở mức 60 ∘ C.

Lưu ý rằng bạn cho kem vào tách cà-phê, và muốn uống càng sớm càng tốt, thì liệu có nên cho kem vào ngay không, hay đợi một lúc? Và nếu cần đợi, thì nên đợi bao lâu?

Để trả lời câu hỏi này, bạn phải mô phỏng quá trình nguội lạnh của nước nóng trong không khí. Cà-phê nóng truyền nhiệt ra ngoài môi trường bằng quá trình đối lưu, bức xạ, và bốc hơi. Việc định lượng từng hiệu ứng nói trên sẽ rất khó và cũng không cần thiết cho việc giải bài toán này.

Để đơn giản hóa, ta có thể dùng Định luật Newton về làm lạnh 2:

df / dt = –r(f – e)
trong đó f là nhiệt độ cà-phê, được viết theo hàm của thời gian còn df / dt là đạo hàm theo thời gian của nó; e là nhiệt độ môi trường, bằng hằng số trong trường hợp này, và r là một thông số (cũng không đổi) để đặc trưng cho tốc độ truyền nhiệt.

Sẽ dễ ước tính r cho một tách cà phê cụ thể bằng cách đo nhiệt độ vài lần theo thời gian. Ta hãy giả sử rằng điều đó đã được thực hiện và r đã tính được bằng 0. 001 với đơn vị là nghịch đảo của giây, 1/s.

• Hãy dùng kí hiệu toán học để biểu diễn hàm tốc độ, g, dưới dạng hàm của y, trong đó y là nhiệt độ của cà-phê tại thời điểm bất kì.
• Tạo ra tập tin M có tên coffee và viết một hàm có tên coffee nhận vào hai biến và không trả lại giá trị đầu ra. Hãy đặt lệnh đơn giản như x=5 vào trong phần thân hàm rồi gọi coffee() từ Command Window.
• Bổ sung một hàm có tên rate_func nhận vào t và y rồi tính g(t, y). Lưu ý rằng trong trường hợp này g thực ra không phụ thuộc vào t; song dù sao hàm cần viết vẫn phải nhận vào t làm đối số đầu tiên để có thể dùng được ode45.

  Hãy thử hàm vừa viết được bằng cách thêm một dòng lệnh như rate_func(0,90) vào trong coffee, rồi gọi coffee từ Command Window.

• Một khi rate_func(0,90) chạy được, hãy sửa coffee để dùng ode45 tính nhiệt độ của cà-phê không (chưa tính đến kem) trong vòng 60 phút. Hãy chắc rằng cà-phê ban đầu nguội nhanh, rồi chậm dần, cuối cùng thì tiến gần về và đạt (xấp xỉ) nhiệt độ phòng sau khoảng 1 giờ.
• Hãy viết một hàm có tên mix_func để tính nhiệt độ nhiệt độ cuối cùng của hỗn hợp hai chất lỏng. Hàm phải nhận vào các tham số là thể tích và nhiệt độ của từng loại chất lỏng.

  Nói chung nhiệt độ của một hỗn hợp thì phụ thuộc vào nhiệt dung riêng của hai chất3. Nhưng nếu làm đơn giản bằng cách giả sử rằng cà-phê và kem có cùng khối lượng riêng và nhiệt dung riêng thì nhiệt độ cuối cùng là (v1y1 + v2y2) / (v1 + v2), trong đó v1 và v2 là các thể tích của hai loại chất lỏng, và y1 và y2 là nhiệt độ của chúng.

  Hãy thêm mã lệnh vào trong coffee để thử nghiệm hàm mix_func vừa viết.

• Hãy dùng mix_func và ode45 để tính thời điểm bắt đầu lúc uống được cà-phê, trong trường hợp rót kem vào ngay.
• Hãy thay đổi coffee để nó nhận vào biến t (có vai trò quyết định sau bao nhiêu giây cà-phê được để nguội trước khi rót kem vào), và trả lại nhiệt độ của tách cà-phê sau khi trộn.
• Hãy dùng fzero tính thời gian t cần để nhiệt độ của tách cà-phê sau khi trộn hạ xuống còn 60∘C.
• Những kết quả trên sẽ cho bạn biết gì về đáp số của bài toán ban đầu? Có phải đó là đáp số như mong muốn không? Những giả thiết nhằm giản hóa nào đã chi phối lời giải này? Theo bạn thì giả thiết nào có ảnh hưởng đáng kể nhất? Bạn có nghĩ rằng ảnh hưởng này đủ mạnh để tác độ đến kết quả không? Xét toàn diện thì bạn tin tưởng đến mức nào vào việc mô hình này sẽ cho ra đáp số cuối cùng? Bạn có thể làm gì để cải thiện mô hình?