Đề cương chi tiết bài giảng Lập trình Matlab

 computing the motion of a nonlinear pendulum 
clear; 
Tspan=[0 25]; 
Theta0 = 0; 
Theta1 = input(' Give initial velocity : '); %%Van toc ban dau 
y0=[Theta0 Theta1]; 
[t,y]=ode45('Pendu', Tspan,y0); 
figure(1);plot(t,y(:,1),'r',t,y(:,2),'b'); 
xlabel('Time');ylabel('Angle, Angular Velocity '); 
pause; 
figure(2); plot(y(:,1),y(:,2)); 
xlabel('Angle');ylabel('Angular Velocity '); 
Kết quả chạy chương trình với Theta1=2 là hai đồ thị: 
Bài toán 3. Xét sự chuyển động của một con tàu vũ trụ dưới tác động của trường hấp 
dẫn của 2 thiên thể. Mỗi thiên thể đều chịu tác động của một lực sinh bởi lực hấp dẫn của 
con tàu. Nhưng khối lượng của con tàu quá nhỏ nên ảnh hưởng không đáng kể đến chuyển 
động của các thiên thể. Vì thế ta có thể bỏ qua ảnh hưởng này. Trong lĩnh vực cơ học vũ trụ, 
bài toán này còn có tên gọi là bài toán ba vật cô lập. Bài toán được áp dụng cho việc tính 
toán các quĩ đạo của con tàu, trái đất và mặt trăng của trong hệ toạ độ (quay) có gốc là tâm 
trường hấp dẫn của trấi đất và mặt trăng. Trong hệ toạ độ này, trái đất nằm ở vị trí (-M,0) và 
mặt trăng ở vị trí (E,0), tàu vũ trụ ở tọa độ (x,y). Các phương trình được cho như sau: 
 2
2 3 3
1 2
2
2 3 3
1 2
( )
2 ; 
2
E x Md x dy M x E
x
dtdt r r
d y dx Ey My
y
dtdt r r
 
   
    
, 
trong đó r1 =  
22 yMx  , r2 =  
22 yEx  và E=1-M. 
Chọn M=0.012277471 (cho hệ mặt trăng- trái đất) và điều kiện đầu: 
x(0) = 1,15; 0
dx
( )
dt
 = 0 ; y(0) = 0; 0
dy
( )
dt
 = 0,0086882909. 
 y
 Earth 
 x Moon 
 M E 
Hình 6.3 Vị trí tương đối của 3 vật thể 
Một lần nữa ta lại rút gọn hệ các phương trình vi phân cấp 2 về hệ phương trình vi 
phân cấp 1 gồm 4 phương trình: 
X1(t) = x(t), X2(t) = 
dx
dt
, Y1(t) = y(t) và Y2(t) =
dy
.
dt 
Hệ phương trình vi phân cần giải là: 
2
1 11
2 1 3 3
1 22
1 2
2 1 1
2 1 3 3
1 2
( ) ( )
2
,
2
X
E X M M X EX
Y X
R RXd
Y Ydt
Y EY MY
X Y
R R
 
 
         
    
   
   
      
 
 
trong đó R1 =  
2
1
2
1 YMX  và R2 =  
2 2
1 1 .X E Y  
Thủ tục giải bài toán như sau: 
- Cài đặt hàm vế phải của hệ phương trình 
% Function defining the Restricted Three-Body 
function yp =Apollo(t,y) 
M=0.012277471 ;E=1-M; 
xx = y(1); yy = y(3); 
r1 = sqrt((xx+M)^2 +yy^2); 
r2 = sqrt((xx-E)^2 +yy^2); 
yp(1) = y(2); 
yp(2) = 2*y(4) +y(1)- E*(y(1)+M)/r1^3-M*(y(1)-E)/r2^3; 
yp(3) = y(4); 
yp(4) = -2*y(2) +y(3)- E*y(3)/r1^3-M*y(3)/r2^3; 
yp=yp.'; 
- Cài đặt chương trình tính toán: 
 % MATLAB code for the Restricted Three-Body Problem 
clear; 
Tspan=[0 23.7]; 
M=0.012277471 ; E=1-M; 
Options= odeset('reltol',1e-5, 'abstol',[ 1e-5 1e-5 1e-5 1e-5]); 
x0=1.15; xdot0= 0; y0=0; 
ydot0 = 0.0086882909; 
V0=[x0 xdot0 y0 ydot0]; 
 [t,y]= ode45('Apollo',Tspan,V0,Options); 
figure(1); plot(y(:,1),y(:,3)); 
axis([-.8 1.2 -.8 .8]); 
grid on; hold on; 
plot(-M,0,'o'); plot(E,0,'o'); 
hold off; 
xlabel('X-axis'); ylabel('Y-axis'); 
text(0, 0.05,'Earth'); text(0.9, -0.15,'Moon'); 
figure(2); 
for i=1:length(t); 
 tt=t(i); 
 xn(i) = cos(tt)*y(i,1) + sin(tt)*y(i,3); 
 yn(i) = -sin(tt)*y(i,1) + cos(tt)*y(i,3); 
 xmoon(i)= E*cos(tt); ymoon(i)= -E*sin(tt); 
 xearth(i) = -M*cos(tt); yearth(i)=M*sin(tt); 
end; 
plot(xn,yn,'r',xmoon,ymoon,'g:',xearth,yearth,'bo'); 
axis([-2 2 -1.5 1.5]);grid on; hold on; 
legend('SpaceCraft','Moon', 'Earth'); 
xlabel('X-axis'); ylabel('Y-axis'); 
plot(xn(1),yn(1),'rx'); hold off; 
 % Animation 
figure(3); 
 for i=1:length(t) 
 plot(xn(i),yn(i),'r+',xmoon(i),ymoon(i),'go',xearth(i),yearth(i),'bo'); 
 axis([-2 2 -1.5 1.5]); 
 pause(0.01); 
 end; 
A. CHỮA BÀI TẬP 
 Sử dụng các hàm nội trú của MATLAB. 
1. Hãy dùng phương pháp Euler để giải gần đúng bài toán Cauchy sau đây với bước đi h = 
0,02: 
' , 1 5
(1) 1
x y
y x
y
y

  

 
2. Dùng phương pháp Hiện ẩn hình thang tính gần đúng hàm y= f(x) cho bởi bài toán 
Cauchy sau đây với h = 0,05: 
' , 0 2
(0) 5
x
y x x
y
y

   

 
3. Dùng phương pháp Hiện ẩn trung điểm giải gần đúng bài toán Cauchy sau đây với h = 
0,1: 
' , 0 10
1
(0) 2
x
y x
y
y

  

 
4. Dùng phương pháp Runge-Kutta giải gần đúng bài toán Cauchy sau đây với bước đi h = 
0,1: 
' , 0 10
2
(0) 0
xy
y x
y

  

 
B. THỰC HÀNH 
 Cài đặt chƣơng trình và lập hàm 
1. Cài đặt chương trình tìm y=f(x) trong đoạn [0,10] thoả mãn phương trình vi phân: 
‟‟‟‟ 2. ‟‟‟. ‟. – 4 ‟‟ 3. . 6 0   y y y y y x y 
với điều kiện đầu        ‟‟‟ 0 ‟‟ 0 1; ‟ 0 0 0,5.   y y y y Sai số tương đối là 10-4 và sai số 
tuyệt đối là 10- 6. 
 Vẽ đồ thị phụ thuộc của y‟‟ đối với y, y‟dưới dạng đường cong Contour 3D. 
2. Cài đặt chương trình giải hệ phương trình vi phân: 
' 2 3
' 3 2
' 7 5
y x g hy
g x hy
h gy xh
  

 
  
trong đoạn x[ 1,10] , với điều kiện đầu là: y(1) =3, g(1)=2 và h(1)=1 với sai số tương đối 
là 10
-5
 và sai số tuyệt đối là 10- 8. 
 Vẽ đồ thị của 3 hàm y, g và h trên cùng một hệ trục toạ độ bằng 3 màu đỏ, xanh 
dương và xanh lá cây. 
3. Cài đặt chương trình giải hệ phương trình vi phân: 
3
2
2 3 5
2
3
dx
t y yz
dt
dy
ty xzt
dt
dz
x txz
dt

  


 


 

trong đoạn t[2,15],với điều kiện đầu là: x(2) =1,5 y(2)=2,34 và z(2)=3,33 
 Vẽ đồ thị đường cong Contour 3D với 3 trục toạ độ là x, y và z. 
4. Tìm 3 hàm y=y(x), g=g(x) và h=h(x) thoả hệ phương trình vi phân: 
' 4
' 5 2
' 2
y xg gyh
g xy hy
h xgy xh
 

 
   
trong đoạn x[ 0,20] , với điều kiện đầu là: y(0) =1, g(0)=1,2 và h(0)=5 với sai số tương 
đối là 10-4 và sai số tuyệt đối là 10- 6. 
 Vẽ đồ thị của 3 hàm y, g và h trên cùng một hệ trục toạ độ bằng 3 màu đỏ, xanh 
dương và xanh lá cây, cùng với các chú thích trên đồ thị. 
5. Cài đặt chương trình giải hệ phương trình vi phân sau: 
3
2
3 2
2 5
3
dx
t y tyz
dt
dy
txy txz
dt
dz
tx xz
dt

  


 


  

trong đoạn t[1,10],với điều kiện đầu là: x(1)=-1,10 y(1)=2,33 và z(1)=5,33. 
 Vẽ đồ thị đường cong Contour 3D với 3 trục toạ độ là x, y và z. 
- Yêu cầu SV chuẩn bị: 
 Cài đặt các chương trình các phương giải bài toán Cauchy. 
- Ghi chú: Đọc các tài liệu tham khảo 1,2,3,4,5. 
 BỘ MÔN DUYỆT ĐỀ CƢƠNG CHI TIẾT BÀI GIẢNG Thay mặt nhóm 
 Chủ nhiệm bộ môn ( Dùng cho 4 tiết giảng) môn học 
 Học phần: Lập trình MATLAB 
 Nhóm môn học: 
 Bộ môn: Toán 
 Tô Văn Ban Khoa : Công nghệ thông tin Nguyễn Trọng Toàn 
Thông tin về nhóm môn học: 
TT Họ tên giáo viên Học hàm Học vị 
1 Nguyễn Trọng Toàn Giảng viên chính TS 
2 Vũ Thanh Hà Giảng viên chính TS 
3 Vũ Anh Mỹ Giảng viên ThS 
Địa điểm làm việc: Bộ môn Toán (1301 nhà S4) 
Điện thoại, email: 069 515 330, bomontoan_hvktqs@yahoo.com 
Bài giảng 15: Phƣơng trình vi phân thƣờng 
Chương, mục: Chương 6, mục 6.5 
Tiết thứ: 57-60 Tuần thứ: 15 
- Mục đích, yêu cầu: 
Hướng dẫn các phương pháp giải phương trình vi phân với điều kiện biên. 
- Hình thức tổ chức dạy học: 
 Lý thuyết, thảo luận, tự học, tự nghiên cứu 
- Thời gian: 
 Lý thuyết trên lớp: 2 tiết Bài tập + Thực hành: 2 tiết 
Tự học, tự nghiên cứu : 4 tiết 
- Địa điểm: Giảng đường do P2 phân công. 
- Nội dung chính: 
6.5 FTVP với điều kiện biên. 
Thực hành chương trình về Hệ phương trình vi phân. 
Ôn tập và giải đáp thắc mắc. 
Kiểm tra. 
 - Yêu cầu SV chuẩn bị: 
 Ôn tập toàn bộ chương đã học. 
- Ghi chú: Đọc các tài liệu tham khảo 1,2,3,4,5. 
6.5 GIẢI PHƢƠNG TRÌNH VI PHÂN VỚI ĐIỀU KIỆN BIÊN 
Mục này sẽ giới thiệu phương pháp bắn (Shooting) để giải một dạng phương trình vi 
phân với điều kiện biên. Đây là một bài toán rất khó nếu như ta không sử dụng các phương 
pháp số và chương trình hóa quá trình giải. 
Giả sử ta muốn giải bài toán Blasius: Tìm hàm y=f(x) xác định trong khoảng 
[0, ) thỏa mãn phương trình vi phân: 
1
''' . '' 0
2
f f f  , với các điều kiện: 
   0 ' 0 0, '( ) 1f f f    . 
 Bài toán này có thể giải được bằng các phương pháp đã giới thiệu ở mục trước nếu ta 
biết điều kiện đầu đối với đạo hàm bậc hai: ''f (0)=?. Thay vào đó ta lại biết điều kiên biên 
vô cùng đối với đạo hàm bậc nhất: '( ) 1f   . 
Để giải bài toán này, người ta thay điều kiện biên vô cùng bằng điều kiện biên tại 
x=10 để thành lập phương trình tìm điều kiện ban đầu phù hợp cho đạo hàm bậc hai. Sau đó 
giải bài toán Cauchy với điều kiện đầu tìm được. Kỹ thuật này được gọi là phương pháp 
bắn. Thủ tục giải bài toán này như sau: 
- Lập hàm về phải cho phương trình của bài toán Blasius: 
% Function for computing solutions to the Blasius equation 
function fder =Blasius(x,f) 
fder(1)= f(2); 
fder(2)= f(3); 
fder(3)= -f(1)*f(3)/2; 
fder=fder.'; 
- Lập hàm Shooting để tìm điều kiện đầu phù hợp cho đạo hàm bậc hai: 
% Function for computing solutions to the Blasius equation 
function dev =Shooting(z) 
global Xinf; 
f0=([ 0 0 z]); 
Xspan = [ 0 Xinf]; 
[x, f] = ode45('Blasius', Xspan,f0); 
n =length(x); 
dev=f(n,2)-1; %% f‟(Xinf)-1 
- Cài đặt chương trình tính toán: 
% MATLAB code computing solutions to the Blasius equation 
clear; 
global Xinf; 
Xinf=10; 
z0=0.5; 
z=fzero('Shooting',z0); %% Giải phương trình f‟(Xinf)=1 
f0=([ 0 0 z]); 
Xspan = [ 0 Xinf]; 
[x, f] = ode45('Blasius', Xspan,f0); 
figure(1); 
plot(f(:,2),x); 
axis([ 0 1.1 0 10]); axis('square'); 
xlabel(' Velocity U/Uinf'); 
ylabel(' Distance from the wall'); 
Chú ý chương trình đã sử dụng khai báo biến toàn cục trong MATLAB: 
global X Y Z ...: Khai báo các biến X, Y và Z là các biến toàn cục. Giá trị ban đầu của 
mỗi biến là một ma trận rỗng. Nếu các hàm có sử dụng các biến này thì cũng phải có khai 
báo GLOBAL. 
 ÔN TẬP VÀ KIỂM TRA 
- Rà soát tất cả các nội dung đã học 
- Giải đáp thắc mắc 
- Hướng dẫn phương pháp thi hết môn học 
- Cho SV làm bài kiểm tra 60 phút 
- Yêu cầu SV chuẩn bị: 
 Ôn tập toàn bộ chương đã học. 
- Ghi chú: Đọc các tài liệu tham khảo 1,2,3,4,5.