Chöông 4 : Ñaïi Soá Tuyeán Tính Vaø Phöông Phaùp Tính Trong Matlab

2. Chöông 4 : Ñaïi Soá Tuyeán Tính Vaø Phöông Phaùp Tính Trong Matlab Chöông goà caùc ñeà muïc sau : Ñaïi soá tuyeán tính trong Matlab. Caùch taïo bieåu thöùc kyù hieäu toùan hoïc trong Matlab. Tính tích phaân xaùc ñònh. Giaûi heä phöông trình vi phaân. Phöông phaùp noäi suy.

3. 4.1) Ma Traän Vaø Ñaïi Soá Tuyeán Tính Trong Matlab : • Heä phöông trình tuyeán tính : Cho heä phöông trình tuyeán tính laø Ax = b Trong ñoù, • A : laø ma traän vuoâng vôùi kích thöôùc laø n×n, • b laø vector heä soá , • x laø vector nghieäm soá cuûa heä phöông trình. • Vector nghieäm soá x cuûa heä phöông trình ñöôïc xaùc ñònh baèng coâng thöùc laø x = A-1b Hay x = inv(A)*b Trong ñoù, inv(A) laø haøm traû veà ma traän ñaûo cuûa ma traän A.

4. Trò rieâng vaø vector rieâng : Cho A laø ma traän thoûa maõn phöông trình laø Laø trò rieâng vaø xi laø vector rieâng töông öùng cuûa A. Phöông trình treân coù theå ñöôïc vieát laïi laø Ñeå tính trò rieâng vaø vector rieâng cuûa A duøng haøm eig cuûa matlab vôùi cuù phaùp laø [V, D] = eig(A) • Moãi vector coät cuûa ma traän V laø moät vector rieâng töông öùng. • D laø ma traän ñöôøng cheùo, moãi phaàn töû töông öùng treân ñöôøng cheùo laø moät giaù trò rieâng.

5. Phöông trình ñaëc tính cuûa ma traän A : Phöông trình ñaëc tính cuûa A ñoù laø moät ñònh thöùc coù daïng laø • Ñeå tìm vector haøng chöùa caùc heä soá an, an-1, . . .,a1 vaø a0 cuûa phöông trình ñaëc tính cuûa ma traän A duøng haøm poly vôùi cuù phaùp laø • P = poly(A) • Ñeå tìm nghieäm cuûa phöông trình ñaëc tính cuûa A duøng haøm roots vôùi cuù phaùp laø roots(P) Trong ñoù, P laø vector haøng chöùa caùc heä soá cuûa ña thöùc phöông trình ñaëc tính cuûa A.

6. 4.2) Caùch Taïo Bieåu Thöùc Kyù Hieäu Toùan Hoïc Trong Matlab : • Ñeå taïo caùc bieåu thöùc kyù hieäu toùan hoïc cho caùc haøm nhö sin, cos hoaëc ma traän söû duïng haøm sym hoaëc leänh syms. • Haøm sym vôùi cuù phaùp laø sym(‘Bieåu thöùc toùan hoïc’). • Leänh syms vôùi cuù phaùp laø syms <Caùc bieán trong bieåu thöùc toùan hoïc>. Ví duï : Caùc leänh sau laø caùch söû duïng sym vaø syms. >>sym(‘sin(2*x)’) sin(2*x) >>diff(sin(2*x)) 2*cos(2*x)

7. >>syms x; >>diff(sin(2*x)) 2*cos(2*x) >>A = sym(‘[a b;c d]’) A = [ a, b] [ c, d] >>syms a b c d; >>A = [a b;c d] A = [ a, b] [ c, d] >>inv(A) ans = [ d/(a*d-b*c), -b/(a*d-b*c)] [ -c/(a*d-b*c), a/(a*d-b*c)]

8. >> det(A) ans = a*d-b*c >> [V,D] = eig(A) V = [ -(1/2*d-1/2*a-1/2*(d^2-2*a*d+a^2+4*b*c)^(1/2))/c, -(1/2*d- 1/2*a+1/2*(d^2-2*a*d+a^2+4*b*c)^(1/2))/c] [ 1, 1] D = [ 1/2*d+1/2*a+1/2*(d^2-2*a*d+a^2+4*b*c)^(1/2), 0] [ 0, 1/2*d+1/2*a-1/2*(d^2-2*a*d+a^2+4*b*c)^(1/2)] >>poly(A) ans = x^2-x*d-a*x+a*d-b*c • Laáy 4 laàn ñaïo haøm cuûa haøm sin(2x) vôùi leänh laø >> diff(sin(2*x),4) ans = 16*sin(2*x) • Taïo ña thöùc vôùi haøm poly2sym >> poly2sym([3 2 4],'y') 3*y^2+2*y+4

9. 4.3) Tính Tích Phaân Xaùc Ñònh : • Tích phaân cuûa haøm moät bieán : Cho tích phaân xaùc ñònh moät bieán laø • Ñeå tính tích phaân xaùc ñònh moät bieán duøng moät trong hai haøm sau : • q = quad(f, a, b, tol, trace) • q = quadl(f, a, b, tol, trace). • Ví duï : Tính ñöôøng cuûa ñöôøng cong vôùi caùc phöông trình cho laø x(t) = sin(2t); y(t) = cos(t); z(t) = t;

10. Chöông trình veõ ñoà thò ñöôøng cong cuûa caùc phöông trình cho treân laø clear t = 0; x = sin(2*t); y = cos(t); z = t p = plot3(x,y,z,'.','EraseMode','none','MarkerSize',5); axis([-1 1 -1 1 0 10]),grid hold on for t = 0:0.01:3*pi x = sin(2*t); y = cos(t); z = t; set(p,'XData',x,'YData',y,'Zdata',z) drawnow pause(0.01) end %end of program

12. Coâng thöùc tính chieàu daøi cuûa ñöôøng cong noùi treân vôùi phöông trình laø • Ñeå tính chieàu daøi cuûa ñöôøng cong vôùi tích phaân cho treân duøng caùc leänh laø • >> f = inline('sqrt(4*cos(2*t).^2 + sin(t).^2 + 1)') • f = • Inline function: • f(t) = sqrt(4*cos(2*t).^2 + sin(t).^2 + 1) • >> q = quad(f,0,3*pi) • q = • 17.2220

13. Tích phaân cuûa haøm hai bieán : Cho tích phaân xaùc ñònh cuûa haøm hai bieán laø Tính tích phaân xaùc ñònh cuûa haøm duøng haøm q = dblquad(f,xmax, xmin, ymax, ymin) Ví duï : Cho tích phaân xaùc ñònh cuûa haøm hai bieán laø

14. >> q = dblquad(inline('y*sin(x)+x*cos(y)'),0,pi,pi,2*pi) q = 29.6088 • Tích phaân cuûa haøm ba bieán : Cho tích phaân cuûa haøm ba bieán laø Ñeå tính tích phaân xaùc ñònh cuûa haøm ba bieán duøng haøm q = triplequad(f,xmin, xmax, ymin, ymax,zmin,zmax) Ví duï : Cho tích phaân xaùc ñònh cuûa haøm ba bieán laø • >> q = triplequad(inline('y*sin(x)+z*cos(x)'),0,pi,0,1,-1,1) • q = • 2.0000

15. 4.4) Nghieäm Cuûa Heä Phöông Trình Vi Phaân Baäc 1 : • Cho heä cuûa caùc phöông trình vi phaân daïng toång quaùt laø Ñeå tính nghieäm cuûa heä phöông trình vi phaân naøy theo thôøi gian duøng moät trong hai haøm sau : [t,x] = ode23(‘odefun_name’,[t0 tfinal],x0) [t,x] = ode45(‘Odefun_name’,[t0 tfinal],x0) Ví duï 1 : Cho heä thoáng nhuùn ñöôïc moâ taû baèng heä phöông trình vi phaân baäc 1 laø

16. Vôùi ñieàu kieän ban ñaàu laø y(0) = 0 vaø v(0) = 0. function der_yv = smd(t,x) global m omega_n zeta; der_yv = zeros(size(x)); der_yv(1) = x(2); der_yv(2) = -(2*omega_n*zeta)*x(2) -(omega_n^2)*x(1); %Beginning of program global m omega_n zeta; m = 1; omega_n = 1; zeta = 0.1; tspan = [0 10*pi]; y0 = 1; v0 = 0; yv0 = [y0 v0]'; [t,yv] = ode45('smd',tspan,yv0); y = yv(:,1); v = yv(:,2); [m,n] = size(yv); [t(m),y(m) v(m)] subplot(2,1,1); plot(t,y,'-',t,y,'o'),grid title('Position vs Time') subplot(2,1,2) plot(y,v,'-',y,v,'o'),grid title('Velocity vs Postion') %End of program

18. Ví duï 2 : Baøi toùan con laéc ngöôïc ñöôïc moâ taû baèng heä phöông trình vi phaân baäc 1 laø • function zdot = invpnl(t,z) • global u M m g len • zdot = zeros(size(z)); • c1 = (M+m); c2 = m*len; c3 = m*g; c4 = (M+m)*len; c5 = (M+m)*g; • zdot(1) = z(2); • top2 = u*cos(z(1)) - c5*sin(z(1)) + c2*cos(z(1))*sin(z(1))*z(2)^2; • zdot(2) = top2/(c2*cos(z(1))^2 - c4); • zdot(3) = z(4); • top4 = u + c2*sin(z(1))*z(2)^2 - c3*cos(z(1))*sin(z(1)); • zdot(4) = top4/(c1-m*cos(z(1))^2);

19. global u M m g len M = 2.0; m = 0.1; len = 0.5; g = 9.81; u = 1; t0 = 0; tff = 1.0; z0 = [0 0 0 0]'; tol = 1.0e-6; options = odeset('RelTol',tol); [tn,znl] = ode23('invpnl',[t0 tff],z0,options); subplot(2,1,1) plot(tn,znl(:,1),'r'),grid title('Step response inverted pendulum') xlabel('Time(Sec)'),ylabel('Tod angle(Radians') subplot(2,1,2) plot(tn,znl(:,2),'b'),grid xlabel('Time(Sec)'),ylabel('car position') %End of program

21. 4.5) Noäi Suy : • Noäi suy moät chieàu : Cho taäp döõ lieäu vaøo ra {xk, yk} cho k = 0, . . ., n vôùi x0 < x1 < x2 < . . . < xn . Noäi suy laø tìm haøm f(x) maø ñoà thò cuûa noù noäi suy töø taäp caùc ñieåm döõ lieäu f(xk) = yk . • Haøm noäi suy moät chieàu vôùi cuù phaùp laø yi = interp1(x, y, xi, method) Trong ñoù, method nhaän moät trong caùc giaù trò sau : • ‘nearest’ : phöông phaùp noäi suy caän gaàn nhaát. • ‘linear’ : phöông phaùp noäi suy tuyeán tính (maëc ñònh). • ‘spline’ hay ‘cubic’ : phöông phaùp noäi suy phi tuyeán. • ‘extra’ : noäi suy ngoøai taàm cuûa x.

22. Ví duï : Cho taäp döõ lieäu veà maät ñoä daân soá cuûa Hoa Kyø tính töø naêm 1900 ñeán 1990 caùch nhau 10 laø >>t = 1900:10:1990; >>p = [75.995 91.972 105.711 123.203 131.669… 150.697 179.323 203.212 226.505 249.633]; • Noäi suy töø taäp döõ lieäu veà maät ñoä daân soá naêm 1975 vôùi leänh laø >> interp1(t,p0,1975) ans = 214.8585 • Bieåu dieãn maät ñoä daân soá töø naêm 1900 ñeán naêm 2000. >> x = 1900:1:2000; >> y = interp1(t,p,x,'spline'); >> plot(t,p,'o',x,y) Cho keát qua ñoà thò noäi suy cuûa y theo x nhö hình

23. Ví duï : Ñoà thò noäi suy töø döõ lieäu sin vôùi caùc leänh laø • >> x = 0:0.2:2*pi; • >> y = sin(x); • >> xi = 0:0.1:2*pi; • >> yi = interp1(x,y,xi,'cubic'); • >> plot(x,y,'ko') • >> plot(xi,yi,'r') • >> plot(x,y,'ko'),hold on,plot(xi,yi,'r:'),hold off,grid

24. Noäi suy hai chieàu : Cho taäp döõ lieäu {xk, yk, zk}, haøm noäi suy töø taäp döõ lieäu naøy laø zi = f(xi, yi) = interp2(x, y, z, xi, yi, method) Ví duï : Ñoà thò noäi suy hai chieàu vôùi caùc leänh laø

25. >> [x,y] = meshgrid(-3:0.25:3); >> z = peaks(x,y); >> [xi,yi] = meshgrid(-3:0.125:3); >> zi = interp2(x,y,z,xi,yi); >> mesh(x,y,z),hold,mesh(xi,yi,zi+15),hold off

26. Ví duï : Cho taäp döõ lieäu tieàn löông lao ñoäng trung bình töø naêm 1950 ñeán naêm 1990 caùch nhau 10 vaø thôøi gian laøm vieäc töø 10 naêm ñeán 30 naêm caùch nhau 10 naêm vôùi caùc leänh laø >>service = 10:10:30; >>years = 1950:10:1990; >>wage = [150.679 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243]; Noäi suy thôøi gian laøm vieäc 15 naêm taïi thjo72i ñieåm naêm 1975 vôùi leänh laø >>w = interp2(service,years,wage,15,1975) w = 190.6287