Bài viết “Phương trình Newton” đã đưa ra phương pháp để giải quyết các bài toán động lực học trong vật lý. Áp dụng thuật toán ấy, ta giải bài toán chuyển động rơi tự do dưới tác dụng của trọng lực đều. Sự rơi tự do cung cấp cho chúng ta một mô hình rất đơn giản của chuyển động: gia tốc luôn là một vec-tơ không đổi:
\[g=9.81\,\mathrm{m/s^2}=\mathrm{const}.\]
Xét một vật bắt đầu rơi từ độ cao \(h_0=50\) m, lấy mặt đất làm gốc toạ độ, chiều dương hướng lên. Chương trình Matlab sau đây cho phép mô phỏng sự rơi ấy ngay trên màn hình máy tí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 | function su_roi_cua_vat clc close all clear all %% CONSTANTS g = 9.81; %% INPUT DATA h0 = 50; v = 0; t = 0; dt = 0.01; %% FIGURE x = 0; y = h0; figure('name','Vat roi','color','white','numbertitle','off'); hold on fig_quanang = plot(x,y,'ro','MarkerSize',20,'markerfacecolor','r'); axis equal axis([-10 10 -1 60]); %% CALCULATION while y>0 t = t+dt; a = -g; v = v+a*dt; y = y+v*dt+0.5*a*dt.^2; set(fig_quanang,'xdata',x,'ydata',y); pause(0.002); end end |
Cần lưu ý về điều kiện trong vòng lặp “while – end” đang sử dụng ở đây. Điều kiện \(y>0\) giúp cho phép tính được lặp đi lặp lại chừng nào vật còn chưa chạm đất. Sau khi chạm đất chương trình tính toán dừng lại. Kết quả mô phỏng là một biểu diễn sinh động theo thời gian:
Khảo sát đầy đủ
Ta có thể mở rộng việc khảo sát bài toán trên nhiều phương diện: quãng đường đi được, diễn biến thay đổi của vận tốc và gia tốc.
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 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 | function roi_tu_do clc close all clear all %% CONSTANTS g = 9.81; %% INPUT DATA y = 50; % Do cao ban dau %% x = 0; v = 0; t = 0; dt = 0.01; s = 0; t_array = t; s_array = s; v_array = abs(v); a_array = g; %% FIGURE fig_vatroi = figure('name','Vat roi','color','white','numbertitle','off'); set(gcf,'Units','normalized'); set(gcf,'Position',[0 0.1 0.2 0.8]); graph_quanang = plot(x,y,'ro','MarkerSize',20,'markerfacecolor','r'); xlabel('X [m]'); ylabel('Y [m]'); axis equal axis([-10 10 -1 60]); fig_quangduong = figure('name','Quang duong','color','white','numbertitle','off'); set(gcf,'Units','normalized'); set(gcf,'Position',[0.2 0.1 0.3 0.38]); graph_quangduong = plot(t_array,s_array,'linewidth',2); xlabel('Thoi gian [s]'); ylabel('Quang duong [m]'); fig_vantoc = figure('name','Van toc','color','white','numbertitle','off'); set(gcf,'Units','normalized'); set(gcf,'Position',[0.5 0.1 0.3 0.38]); graph_vantoc = plot(t_array,v_array,'linewidth',2); xlabel('Thoi gian [s]'); ylabel('Van toc [m/s]'); fig_giatoc = figure('name','Gia toc','color','white','numbertitle','off'); set(gcf,'Units','normalized'); set(gcf,'Position',[0.2 0.5 0.3 0.38]); graph_giatoc = plot(t_array,a_array,'linewidth',2); xlabel('Thoi gian [s]'); ylabel('Gia toc [m/s^2]'); %% CALCULATION while y>0 t = t+dt; t_array = [t_array t]; a = -g; v = v+a*dt; y = y+v*dt+0.5*a*dt.^2; s_array = [s_array 50-y]; v_array = [v_array abs(v)]; a_array = [a_array abs(a)]; figure(fig_vatroi); set(graph_quanang,'xdata',x,'ydata',y); figure(fig_quangduong); set(graph_quangduong,'xdata',t_array,'ydata',s_array); figure(fig_vantoc); set(graph_vantoc,'xdata',t_array,'ydata',v_array); figure(fig_giatoc); set(graph_giatoc,'xdata',t_array,'ydata',a_array); pause(0.002); end end set(graph_giatoc,'xdata',t_massiv,'ydata',a_array); pause(0.002); end end |