69
Dưới đây là ba m-file được soạn.
% file_name = ptvp_vantoc.m
function dydt=ptvp_vantoc(t,y)
global m b F0 % khai bao cac bien tong the
dydt=1/m*(F0-b*y);
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
% file_name = Euler_method.m
function [t, z]=Euler_method(ptvp, t0, tf, y0, N)
% Cac bien:
% ptvp - ten file chua ve phai cua phuong trinh vi phan can
giai
% t0 - thoi diem dau
% tf - thoi diem cuoi (t0 < tf)
% y0 - gia tri ham tai thoi diem dau
% N - so doan ma khoang thoi gian (tf-t0) duoc chia ra
global m b F0
% khai bao cac bien tong
h = (tf - t0)/N; % buoc thoi gia
t = t0+[0:N]'*h;
% vector chua cac diem c
z(1,:) = y0;
for k = 1:N
z(k + 1,:) = z(k,:)+h*ptvp(t(k),z(k,:));
end
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
% file_name = tuy_chon__.m
% chuong trinh chinh
global m b F0 % khai bao cac bien tong the
m=1; b=0.2; F0=4; % Gan gia tri cho cac tham so
t0=0; tf=10; N=20; y0=0;
% nghiem gan dung z
[t, z]=Euler_method(@ptvp_vantoc, t0, tf, y0, N);
%nghiem chinh xac
y=F0/b*(1-exp(-b/m.*t));
figure(1)
subplot(211)
plot(t,y,'k',t,z,'-.k','Linewidth',1);
% so sanh ket qua tinh voi nghiem chinh xac
legend('y','z')
ylabel('y, z [m/s]');
subplot(212)
plot(t,abs(z-y),'k','Linewidth',1);
% sai so cua phuong phap
xlabel('t [s]'); ylabel('|y-z|');
