9/3/2019
1
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 1
Giảng viên: Ths. Trần Kim Bằng
Bộ môn Cơ Kỹ Thuật, P.106B4
Khoa Khoa Học Ứng Dụng
Đại học Bách Khoa TpHCM
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 2
Chương 8
Các phép toán trên đa thức
và biến biểu tượng
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 3
8.1. Đa thức trong MATLAB
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 4
Trong Matlab, có thể biểu diễn một hàm đa thức dưới dạng mảng một
chiều như sau
p=[1 -12 25 1 -116]
Trong đó, các phần tử trong mảng p là các hệ số nhân với x mũ n với n
giảm dần từ trái qua phải. Như vậy, ta có thể ngầm hiểu mảng p đang biểu
diễn một hàm đa thức
p = x4 - 12x3 + 25x2 + x - 116
Có thể biểu diễn đa thức p trên đồ thị bằng cách kết hợp lệnh polyval và
lệnh plot
p=[1 -12 25 1 -116]
x=linspace(-1,3,500)
y=polyval(p,x) %polyval là hàm để tính giá trị của đa thức p tại một giá trị x
plot(x,y)
8.1.1. Định nghĩa
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 5
Để tìm nghiệm của đa thức, dùng lệnh roots
p=[1 -12 25 1 -116]
r=roots(p)
8.1.2. Tìm nghiệm của đa thức
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 6
Nếu biết trước được tập nghiệm, cũng có thể tìm lại đa thức từ tập nghiệm
đó qua lệnh poly. Ta có thể truy ngược lại các hệ số của đa thức p thông
qua lệnh poly
Ví dụ: Dùng máy tính bỏ túi để tìm nghiệm của đa thức 7x^3+5x^2+3x+1
r=[-0.45881893;-0.127733391+0.543178251i;-0.127733391-0.543178251i]
p=poly(r)
p=p*7
8.1.3. Tìm hệ số đa thức từ nghiệm
9/3/2019
2
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 7
Cộng đa thức
Hai đa thức cùng bậc:
p1=[2 5 4 1];
p2=[1 -3 9 3];
p = p1 + p2
Tương tự cho trừ đa thức
Hai đa thức khác bậc:
p1=[2 5 4 1];
p2=[1 0 9 3];
p = p1– p2;
Thêm các hệ số 0 vào đa thức có bậc thấp hơn để
hai đa thức có cùng bậc
8.1.4. Cộng trừ hai đa thức
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 8
Hàm conv thực hiện nhân hai đa thức
Ví dụ:
Nhân 2 đa thức sau:
A(x)=x3+2x2+3x+4 và b(x)=x3+4x2+9x+16
Viết lại hai đa thức này dưới dạng mảng một chiều trong Matlab
a=[1 2 3 4]; b=[1 4 9 16];
c=conv(a,b) %dùng lệnh conv để nhân hai đa thức
8.1.5. Nhân hai đa thức
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 9
Chia đa thức:
Ví dụ: Cho 2 đa thức:
y = x3 - 5x2 + 4x - 9
z = x2 + 3
Viết lại dưới dạng mảng một chiều trong Matlab
y = [1 -5 4 -9];
z = [1 0 3];
p = deconv(y,z)
p = x – 5, là phần thương
[p,r] = deconv(y,z)
p = x – 5, là phần thương
r = x + 6, là phần dư
Vậy y = pz + r = (x – 5)(x2 + 3) + (x + 6)
8.1.6. Chia hai đa thức
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 10
y = x3 - 5x2 + 4x - 9
y = [1 -5 4 -9]
Dùng lệnh polyder để tính đạo hàm
đa thức
z = polyder(y)
8.1.7. Đạo hàm đa thức
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 11
8.2. Các phép toán trên
biến biểu tượng
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 12
syms a b c x
f=a*x^2+b*x+c;
solve(f) %giải phương trình
Nếu lấy đạo hàm f
diff(f)
Sẽ cho ra kết quả như sau
ans = b + 2*a*x
Tại sao Matlab chọn biểu tượng x là biến để lấy đạo hàm mà
không chọn a hay b? Matlab mặc định chọn biểu tượng gần nhất
với chữ cái x làm biến.
có thể dùng lệnh findsym(f,n) để tìm n biến trong hàm f gần nhất
với chữ cái x
8.2.1. Khai báo biến biểu tượng
9/3/2019
3
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 13
Thế hàm vào biến bằng lệnh compose
syms a b c x y
f = a*x^2+b*x+c
g = sin(y);
h=compose(f,g)
8.2.2. Thế một hàm vào biến biểu tượng
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 14
Lệnh subs dùng để thay thế một giá trị hay một biểu thức vào biến
symbolic
syms a b alpha
f=cos(a)+sin(b);
g=subs(f,{a,b},{alpha,2}) % thế alpha vào a và 2 vào b
g1=subs(f,{a,b},{[2 3],[12 8]})
Ta được: g = sin(2) + cos(alpha)
8.2.3. Thế một giá trị hay một biểu thức vào biến biểu tượng
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 15
Biểu diễn biểu thức toán học trên đồ thị bằng lệnh latex
syms x y t
F1 = diff(atan(x));
F2 =(x^2+1/x)/sqrt((sin(x)+cos(x))/x^3.5);
H1=latex(F1);
H2=latex(F2);
text(.2,.25,['$$' H1 '$$'],'interpreter','latex','fontsize',14)
text(.6,.2, ['$$' H2 '$$'],'interpreter','latex','fontsize',16)
8.2.4. Biểu diễn biểu thức toán học trên đồ thị
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 16
Tính giới hạn bằng lệnh limit
Tính giới hạn của hàm f khi h tiến đến 0
syms x h n
limit((cos(x+h)-cos(x))/h,h,0)
Tính giới hạn của hàm g khi n tiến về vô cùng
g = (1+x/n)^n
limg = limit(g,n,inf)
8.2.5. Tính giới hạn hàm số
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 17
Tính giới hạn bên phải hay bên trái (0+ hay 0-) của hàm 1/x
syms x
lim_plus=limit(1/x,x,0,'right')
lim_minus=limit(1/x,x,0,'left')
Trong trường hợp này, nếu ta không ghi rõ giới hạn trái hay phải thì kết
quả sẽ không xác định
syms x
limit(1/x,x,0)
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 18
Dùng lệnh finverse để tìm hàm ngược
syms x
f=tan(1/x)
finverse(f)
8.2.6. Tìm hàm ngược
9/3/2019
4
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 19
Giải phương trình theo biến symbolic bằng lệnh solve
syms a b c x
f=a*x^2+b*x+c;
y=solve(f)
Ta có thể giải theo biến ta mong muốn, giả sử là biến b
syms a b c x
f=a*x^2+b*x+c;
y=solve(f,b)
8.2.7. Giải phương trình theo biến biểu tượng
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 20
syms a b c d x y t
[x,y]=solve(a*x+b*y-t,c*x+d*y)
Hoặc ta có thể dùng cách khác
syms a b c d x y t
g=solve(a*x+b*y-t,c*x+d*y)
g.x
g.y
8.2.8. Giải hệ phương trình theo biến biểu tượng
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 21
Thu gọn biểu thức bằng lệnh collect
syms x t b
f=(1+x)*t*b^2+x*t*b^3+b
f1=collect(f)
f2=collect(f,x) % Đặt x ra ngoài
f3=collect(f,t) % Đặt t ra ngoài
f4=collect(f,b) % Đặt b ra ngoài
Khai triển biểu thức bằng lệnh expand
syms x y
g=cos(x+y);
g=expand(g)
8.2.9. Thu gọn và khai triển biểu thức
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 22
Lệnh factor biến đổi biểu thức từ dạng tổng thành tích
syms x
f=x^3-6*x^2+11*x-6;
f=factor(f)
Ta được: f = (x-3)*(x-1)*(x-2)
syms x
g=x^6+1;
g=factor(g)
Ta được: g = (x^2+1)*(x^4-x^2+1)
8.2.10. Biến đổi tổng thành tích
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 23
Tối giản biểu thức bằng lệnh simplify
syms x
f=(1-x^2)/(1-x);
f=simplify(f)
8.2.11. Tối giản biểu thức
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 24
Lọc ra các hệ số đi kèm biến x bằng lệnh coeffs
syms x
f=x^3-2*x^2+3*x-4;
c=coeffs(f,x)
Ta được: c =[ -4, 3, -2, 1]
8.2.12. Lấy các hệ số trong đa thức
9/3/2019
5
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 25
Lệnh lấy đạo hàm diff
syms a x
f=sin(a*x^2);
diff(f) %lấy đạo hàm của f theo biến gần nhất với ký tự x dựa trên quy
luật alphabet
diff(f,a) %lấy đạo hàm của f theo biến ta mong muốn là biến a
diff(f,2) %lấy đạo hàm cấp 2 của f theo biến gần nhất với ký tự x dựa trên
quy luật alphabet
diff(f,a,2) %lấy đạo hàm cấp 2 của f theo biến ta mong muốn là biến a
Lấy đạo hàm cho mảng hai chiều
syms a x
A=[cos(a*x),sin(a*x); -sin(a*x),cos(a*x)];
diff(A) %lấy đạo hàm theo biến x
8.2.13. Đạo hàm theo biến biểu tượng
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 26
Lệnh lấy tích phân int
syms x
f=sin(x);
F=int(f) %lấy tích phân của f theo x
F1=int(f,0,pi) %lấy tích phân của f theo x trên khoảng từ 0 đến pi
syms x y
f=sin(x) +y^2;
F=int(f,y) %lấy tích phân của f theo y
F1=int(f,y,0,2) %lấy tích phân của f theo y trên khoảng từ 0 đến pi
8.2.14. Tích phân theo biến biểu tượng
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 27
dsolve('Df = f + sin(t)')
Giải phương trình vi phân không có điều kiện ban đầu
sin( )
df
f t
dt
8.2.15. Tìm nghiệm tường minh của phương trình vi phân
Giải phương trình vi phân có điều kiện ban đầu y =1 lúc t = 0
2
1 , (0) 1
dy y y
dt
dsolve('Dy=1+y^2','y(0)=1')
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 28
Giải hệ phương trình vi phân có điều kiện ban đầu x=2 và y =1 lúc t = 0
s = dsolve('Dy=1+2*y','Dx=4*y+5*x','x(0)=2','y(0)=1')
1 2 , 4 5 , x 0 2 , (0) 1
dy dx
y y x y
dt dt
Để xem kết quả x(t) và y(t) gọi lệnh s.x và s.y
z = dsolve('Dx = y', 'Dy = -x')
Giải hệ phương trình vi phân không có điều kiện ban đầu
(t)
(t)
dx
y
dt
dy
x
dt
Matlab sẽ ngầm hiểu y và x là hàm theo t
8.2.16. Tìm nghiệm tường minh của hệ phương trình vi phân
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 29
s = dsolve('Dy=1+2*y','Dx=4*y+5*x','x(0)=2','y(0)=1')
• Để có thể dễ dàng thay đổi các tham số 2, 4 và 5, ta dùng cách ghép
chuỗi
• Bằng cách sử dụng lệnh strcat để ghép chuỗi và num2str để thay đổi
tham số
• Gán các phương trình vi phân vào một biến nào đó, giả sử đặt tên là
Eq1 và Eq2
• Gán điều kiện ban đầu vào một biến nào đó, giả sử đặt tên là Intcond
a = 2; b = 4; c = 5
Eq1 = strcat('Dy=1+',num2str(a),'*y')
Eq2 = strcat('Dx=',num2str(b),'*y+',num2str(c),'*x')
Intcond = 'x(0)=2,y(0)=1'
s = dsolve(Eq1, Eq2, Intcond)
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 30
Để tính x(t) và y(t) theo một dãy thời gian t=0:0.01:5, ta dùng các lệnh sau
syms t
t=0:0.05:5;
xt = eval(s.x);
plot(t,xt);
yt = eval(s.y);
figure
plot(t,yt);