Xem thêm: Chương trình mô phỏng chuyên nghiệp về con lắc lò xo viết bằng java.
Tại bài viết này, chúng ta phân tích hoạt động của con lắc lò xo trên nền Matlab. Với một con lắc nằm ngang, xét theo phương chuyển động, ta có phương trình động lực học đơn giản:
\(-kx=ma.\tag{1}\)
Tham số ban đầu của hệ do người dùng cài đặt:
1 2 3 4 5 6 7 | %% INPUT DATA l0 = 0.1; A = 0.02; m = 0.1; k = 1; dt = 0.01; tmax = 100; |
Trong vòng lặp tính toán sự tiến hoá của trạng thái theo thời gian, gia tốc \(a\) được gán theo quy luật của phương trình (1):
1 2 3 4 5 6 | t = t+dt; t_array = [t_array t]; a = -k*x/m; v = v+a*dt; x = x+v*dt+0.5*a*dt*dt; |
Ngoài ra code chương trình cung cấp dưới đây còn giúp vẽ các đồ thị phụ thuộc vào thời gian của li độ, vận tốc, gia tốc. Đặc biệt quỹ đạo pha của chuyển động cũng được diễn tả qua hàm vẽ:
1 |
Video minh hoạ
Code chương trì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 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 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 | function conlacloxo clc close all clear all %% INPUT DATA l0 = 0.1; A = 0.02; m = 0.1; k = 1; dt = 0.01; tmax = 100; %% Xu ly x = A; v = 0; t = 0; y = 0; a = -k*x/m; t_array = t; x_array = x; v_array = v; a_array = a; %% FIGURE fig1 = figure('name','Con lac lo xo','color','black','numbertitle','off'); hold on set(gcf,'Units','normalized'); set(gcf,'Position',[0 0.1 0.5 0.32]); fig_loxo = veloxo(-l0,x,(l0+A)/20,10,[0.5 0.5 0.5],'',false); fig_quanang = plot([x y],'ro','MarkerSize',30,'markerfacecolor','r'); l_vector_danhoi = l0/3; ve_vec_danhoi = drawArrow('',[x 0],[x-50*k*x*l_vector_danhoi 0],'m',false); l_vector_v = l0/3; ve_vec_v = drawArrow('',[x 0],[x-20*l_vector_v*v 0],'y',false); ht = title(sprintf('t = %0.2f s',t)); set(gca,'color','k','xcolor','w','ycolor','w') axis equal axis([-l0 3*A -0.15*l0 0.15*l0]); fig_x = figure('name','Toa do','color','white','numbertitle','off'); set(gcf,'Units','normalized'); set(gcf,'Position',[0.5 0.1 0.3 0.38]); graph_x = plot(t_array,x_array,'k','linewidth',2); xlabel('t [s]'); ylabel('x [m]'); fig_v = figure('name','Van toc','color','white','numbertitle','off'); set(gcf,'Units','normalized'); set(gcf,'Position',[0.5 0.52 0.3 0.38]); graph_v = plot(t_array,v_array,'y','linewidth',2); xlabel('t [s]'); ylabel('v [m/s]'); fig_a = figure('name','Gia toc','color','white','numbertitle','off'); set(gcf,'Units','normalized'); set(gcf,'Position',[0.2 0.52 0.3 0.38]); graph_a = plot(t_array,a_array,'m','linewidth',2); xlabel('t [s]'); ylabel('a [m/s^2]'); fig_x_v = figure('name','Gia toc','color','white','numbertitle','off'); set(gcf,'Units','normalized'); set(gcf,'Position',[0 0.52 0.2 0.38]); graph_x_v = plot(x_array,v_array,'k','linewidth',1); xlabel('x [m]'); ylabel('v [m/s]'); %% CALCULATION while t<tmax t = t+dt; t_array = [t_array t]; a = -k*x/m; v = v+a*dt; x = x+v*dt+0.5*a*dt*dt; x_array = [x_array x]; v_array = [v_array v]; a_array = [a_array a]; figure(fig1); veloxo(-l0,x,(l0+A)/20,10,[0.5 0.5 0.5],fig_loxo,true); set(fig_quanang,'xdata',x,'ydata',y); set(ht,'string',sprintf('t = %0.2f s',t)); drawArrow(ve_vec_danhoi,[x 0],[x-50*k*x*l_vector_danhoi 0],'m',true); drawArrow(ve_vec_v,[x 0],[x+20*l_vector_v*v 0],'y',true); figure(fig_x); set(graph_x,'xdata',t_array,'ydata',x_array); figure(fig_v); set(graph_v,'xdata',t_array,'ydata',v_array); figure(fig_a); set(graph_a,'xdata',t_array,'ydata',a_array); figure(fig_x_v); set(graph_x_v,'xdata',x_array,'ydata',v_array); pause(0.002); end end %% Ham ve lo xo function loxo = veloxo(x1,x2,A,q,color,fun,velai) dis = abs(x2-x1); lambda = dis/q; N = (q+1)*2; x_lx = linspace(x1-lambda/4,x2,N); x_lx(1) = x_lx(1)+lambda/4; x_lx(N) = x_lx(N)-lambda/4; y_lx = zeros(1,N); for i = 2:N-1 if mod(i,2)==0 y_lx(i) = A; else y_lx(i) = -A; end end if velai == false loxo = plot(x_lx,y_lx,'y','linewidth',2); loxo.Color = color; else set(fun,'xdata',x_lx,'ydata',y_lx,'color',color); end end %% Ham ve mui ten function hArrow = drawArrow(fun,p0,p1,color,draw) % drawArrow(p0,p1) % draw: boolean if nargin == 2 color = 'k'; end % Parameters: W1 = 0.08; % half width of the arrow head, normalized by length of arrow W2 = 0.014; % half width of the arrow shaft L1 = 0.18; % Length of the arrow head, normalized by length of arrow L2 = 0.13; % Length of the arrow inset % Unpack the tail and tip of the arrow x0 = p0(1); y0 = p0(2); x1 = p1(1); y1 = p1(2); % Start by drawing an arrow from 0 to 1 on the x-axis P = [... 0, (1-L2), (1-L1), 1, (1-L1), (1-L2), 0; W2, W2, W1, 0, -W1, -W2, -W2 ]; % Scale,rotate, shift and plot: dx = x1-x0; dy = y1-y0; Length = sqrt(dx*dx + dy*dy); Angle = atan2(-dy,dx); P = Length*P; %Scale P = [cos(Angle), sin(Angle); -sin(Angle), cos(Angle)]*P; %Rotate P = p0(:)*ones(1,7) + P; %Shift % Plot! if draw == 0 hArrow = patch(P(1,:), P(2,:),color); hArrow.EdgeColor = color; else set(fun,'xdata',P(1,:),'ydata',P(2,:),'EdgeColor',color); end end |