Phương trình tiến hoá

Dẫn nhập

Những phương trình cơ bản trong mỗi lĩnh vực của vật lý học thường biểu diễn quy luật biến đổi của trạng thái theo không gian và thời gian. Trong cơ học cổ điển, phương trình cơ bản của động lực học chất điểm biểu diễn qua định luật Newton thứ hai:

\frac{d^2x}{dt^2}=f(x,\dot{x})

trong đó f(x,\dot{x}) - hàm số của lực tác dụng lên một đơn vị khối lượng. Sự lan truyền xung được diễn tả qua phương trình sóng:

\frac{\partial^2u}{\partial x^2}=\frac{1}{v^2}\frac{\partial^2u}{\partial t^2},

trong đó v - vận tốc truyền sóng. Quá trình phân bố nhiệt độ diễn ra theo phương trình truyền nhiệt:

\frac{\partial T}{\partial t}=\chi\frac{\partial^2T}{\partial x^2},

với \chi - độ dẫn nhiệt của môi trường. Trong điện động lực học, hệ phương trình Maxwell đóng vai trò nền tảng:

\begin{eqnarray}
\nabla\cdot\vec{D}&=&\rho,\nonumber\\
\nabla\cdot\vec{B}&=&0,\nonumber\\
\nabla\times\vec{E}&=&-\frac{\partial\vec{B}}{\partial t},\nonumber\\
\nabla\times\vec{H}&=&\vec{j}+\frac{\partial\vec{D}}{\partial t}.\nonumber
\end{eqnarray}

Ở đây dấu \nabla=\dfrac{\partial}{\partial x}\vec{i}+\dfrac{\partial}{\partial y}\vec{j}+\dfrac{\partial}{\partial z}\vec{k} là toán tử Hamilton. Trong tĩnh điện học, phương trình Poisson diễn tả quy luật phân bố điện trường:

\nabla^2\varphi=\frac{\rho}{\varepsilon_0}.

Trong cơ học lượng tử cổ điển, sự tiến hoá của hàm sóng diễn ra theo quy luật của phương trình Schrodinger:

i\hbar\frac{\partial}{\partial t}\psi=\left[-\frac{\hbar^2}{2m}\nabla^2+U(\vec{r},t)\right]\psi,

trong đó U(\vec{r},t) - hàm thế năng.

Như vậy, việc giải các bài toán vật lý thường quy về việc giải phương trình vi phân thông thường và phương trình vi phân đạo hàm riêng. Tuy nhiên trong nhiều trường hợp ngoài thực thế khoa học - kĩ thuật, việc đi tìm nghiệm giải tích hết sức khó khăn và đôi khi không cần thiết. Ngày nay với công nghệ máy tính phát triển, các phương pháp tính toán số được ứng dụng rộng rãi. Matlab là một công cụ rất tốt để giải các bài toán đó.

Bài toán

Bài toán sau đây sẽ là một ví dụ cho thuật toán giải phương trình vi phân thông thường (ODE - Ordinary Differential Equation) thực thi trên Matlab. Vấn đề được trình bày theo tuần tự: Đặt vấn đề - Mô hình hoá - Thuật toán - Giải bằng hàm chuyên dụng.

Xét một hệ sinh thái đơn giản gồm thỏ và sói. Thỏ ăn cỏ để sống, còn sói chỉ ăn thịt thỏ. Số lượng bầy thỏ và bầy sói sẽ thay đổi như thế nào theo thời gian?

Bài toán sinh thái học này đã được nghiên cứu từ lâu và trở thành kinh điển trong các sách giáo khoa về các hệ động lực và cân bằng sinh thái. Số liệu thực nghiệm về quần thể thỏ - sói được ghi lại trong biểu đồ sau:

Đồ thị diễn biến số lượng cá thể thỏ và sói

Dưới đây chúng ta sẽ lặp lại phân tích và giải thích diễn biến thực nghiệm ấy.

Mô hình hoá

Gọi y_1 là số thỏ, y_2 là số sói. Mục đích của chúng ta là tính số lượng thỏ và sói tại một thời điểm nào đó trong tương lai. Để làm được điều đó, ta cần đến hàm số miêu tả tốc độ biến đổi \dfrac{dy_1}{dt}\dfrac{dy_2}{dt} theo thời gian.

Có thể xem rằng lượng cỏ đủ nhiều thỏ ăn mãi không hết. Thỏ mẹ sẽ sinh thỏ con, thỏ con lớn lên sẽ sinh thêm thỏ mới. Dễ hình dung rằng số lượng thỏ con sinh ra trong một đơn vị thời gian tỉ lệ với chính dân số bầy thỏ:

\begin{equation}
\frac{dy_1}{dt}=\alpha y_1,
\label{eq:tho1}
\end{equation}

trong đó \alpha - hệ số tỉ lệ phụ thuộc vào môi trường. Phương trình này có nghiệm giải tích đơn giản:

y_1=y_1(0)e^{\alpha t}.

Theo đó số thỏ không ngừng tăng lên theo hàm mũ! Thực tế tất nhiên không diễn ra như thế, thỏ sẽ bị sói ăn thịt bớt. Càng nhiều sói, thỏ càng bị ăn thịt nhiều. Ngoài ra số lượng thỏ bị ăn thịt cũng phải tỉ lệ với chính dân số của chúng. Phương trình \eqref{eq:tho1} cần viết lại dưới dạng:

\begin{equation}
\frac{dy_1}{dt}=\alpha y_1-\beta y_1y_2.
\label{eq:tho2}
\end{equation}

Bây giờ xem qua đàn sói. Nếu không có thỏ sói sẽ bị thiếu thức ăn, lượng sói sẽ suy giảm với tốc độ tỉ lệ với chính dân số của chúng. Dân số càng cao thì số bị chết đi trong một đơn vị thời gian càng lớn:

\begin{equation}
\frac{dy_2}{dt}=-\gamma y_2.
\label{eq:soi1}
\end{equation}

Chúng sẽ chết dần theo quy luật:

y_2=y_2(0)e^{-\gamma t}.

Tất nhiên sói không chết hết mà vẫn duy trì quần thể nhờ ăn thịt thỏ. Thỏ càng nhiều, sói càng nhiều thức ăn và càng sinh sôi. Ngoài ra càng nhiều sói mẹ thì sói con cũng càng thêm được sinh ra. Phương trình \eqref{eq:soi1} cần viết thành

\begin{equation}
\frac{dy_2}{dt}=-\gamma y_2+\delta y_1y_2.
\label{eq:soi2}
\end{equation}

Từ \eqref{eq:tho2} và \eqref{eq:soi2} ta có được hệ phương trình cần thiết:

\begin{eqnarray}
\frac{dy_1}{dt}&=&\alpha y_1-\beta y_1y_2,\label{eq:tho3}\\
\frac{dy_2}{dt}&=&-\gamma y_2+\delta y_1y_2.\label{eq:soi3}
\end{eqnarray}

Hệ phương trình này do Lotka và Volterra đưa ra vào năm 1925.

Thuật toán

Chương trình dưới đây cho phép miêu tả sinh động số lượng thỏ và sói theo thời gian:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
function tho_soi
clc;
clear all
close all

%% INPUT DATA
alpha = 0.07; beta = 0.0025; gamma = 0.2; delta = 0.003 ;

y1 = 30;
y2 = 10;

t = 0;
tmax = 200;
dt = 0.05;

%% FIGURE
fig1 = figure('name','Tho va soi','color','white','numbertitle','off');
hold on
h_tho = plot([t y1],'bo','MarkerSize',10,'markerfacecolor','b');
h_soi = plot([t y2],'ro','MarkerSize',10,'markerfacecolor','r');
ht = title(sprintf('t = %0.2f days',t));

axis([0 tmax -1 200]);

%% CALCULATION
while t<tmax
t = t+dt;
dy1_dt = (alpha-beta*y2)*y1;
dy2_dt = (-gamma+delta*y1)*y2;
y1 = y1+dy1_dt*dt;
y2 = y2+dy2_dt*dt;

figure(fig1)
plot(t,y1,'bo','markersize',2,'markerfacecolor','b');
plot(t,y2,'ro','markersize',2,'markerfacecolor','r');
set(h_tho,'xdata',t,'ydata',y1);
set(h_soi,'xdata',t,'ydata',y2);
set(ht,'string',sprintf('t = %0.2f days',t));
pause(0.002);
end
end

Đầu chương trình ta cần khai báo các tham số \alpha, \beta, \gamma, \delta của môi trường cũng như số lượng ban đầu của thỏ và sói:

1
2
3
4
5
6
7
8
9
%% INPUT DATA
alpha = 0.07; beta = 0.0025; gamma = 0.2; delta = 0.003 ;

y1 = 30;
y2 = 10;

t = 0;
tmax = 200;
dt = 0.05;

Trong vòng lặp while - end, số lượng thỏ và sói được tính toán theo đúng quy luật của hệ phương trình \eqref{eq:tho3} - \eqref{eq:soi3}

1
2
3
4
5
6
7
while t<tmax
t = t+dt;
dy1_dt = (alpha-beta*y2)*y1;
dy2_dt = (-gamma+delta*y1)*y2;
y1 = y1+dy1_dt*dt;
y2 = y2+dy2_dt*dt;
end

Hàm chuyên dụng của Matlab

Thuật toán nói trên giải phương trình vi phân theo phương pháp Euler, có sai số lớn, nhưng rất tiện để minh hoạ cho người mới tập giải hiểu được nguyên lý hoạt động. Có những phương pháp chính xác hơn, kéo dài và phức tạp hoá dòng lệnh hơn. Ở đây chúng ta sẽ dùng luôn các hàm giải phương trình vi phân chuyên dụng sẵn có trong Matlab.

Dưới đây là sơ đồ sử dụng lệnh ode45 trong Matlab để giải phương trình vi phân bằng phương pháp Runge-Kutta bậc 4-5. Hàm ode45 nằm trong chương trình chính sẽ gọi đến hệ phương trình vi phân cần giải nằm trong chương trình con.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function tho_soi_ode_system

global alpha beta gamma delta

alpha = 0.07; beta = 0.0025; gamma = 0.2; delta = 0.003;

......
[T Y] = ode45(@hethosoi,linspace(0,tmax,500),[x0 y0]);

......

end

function dy = hethosoi(t,y)
global alpha beta gamma delta
dy = zeros(2,1);
dy(1) = (alpha-beta*y(2))*y(1);
dy(2) = (-gamma+delta*y(1))*y(2);
end

Hàm ode45 có 3 tham số chính: tên của hệ phương trình cần giải đặt sau dấu @, tập hợp các thời điểm cần tính toán trạng thái, giá trị ban đầu của trạng thái.

Code chương trình hoàn chỉnh bên dưới đưa vào những lệnh giúp biểu diễn kết quả được thuận tiện.

Video minh hoạ

Code chương trình hoàn chỉnh

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
function tho_soi_ode_system
clc
clear all
close all

global alpha beta gamma delta

%% INPUT DATA
alpha = 0.07;
beta = 0.0025;
gamma = 0.2;
delta = 0.003;

x0 = 30;
y0 = 10;

tmax = 200;

x = x0;
y = y0;
t = 0;

%% CALCULATION
[T Y] = ode45(@hethosoi,linspace(0,tmax,500),[x0 y0]);

%% FIGURE
fig1 = figure('name','Tho va soi','color','white','numbertitle','off');
plot1 = animatedline('linewidth',2,'color','b');
plot2 = animatedline('linewidth',2,'color','r');
axis([0 tmax -inf inf]);
for i = 1:length(T)
addpoints(plot1,T(i),Y(i,1))
addpoints(plot2,T(i),Y(i,2))
drawnow limitrate
pause(0.05)
end

fig2= figure('name','Tho va soi','color','white','numbertitle','off');
plot1 = animatedline('linewidth',2,'color','b');
for i = 1:length(T)
addpoints(plot1,Y(i,1),Y(i,2))
drawnow limitrate
pause(0.05)
end
end

function dy = hethosoi(t,y)
global alpha beta gamma delta
dy = zeros(2,1);
dy(1) = (alpha-beta*y(2))*y(1);
dy(2) = (-gamma+delta*y(1))*y(2);
end