intTypePromotion=1
zunia.vn Tuyển sinh 2024 dành cho Gen-Z zunia.vn zunia.vn
ADSENSE
 » 
 » 

Bài giảng Matlab ứng dụng: Phần II - TS. Nguyễn Hoài Sơn

Chia sẻ: Nguyễn Thị Ngọc Lựu | Ngày: | Loại File: PDF | Số trang:45

Thêm vào BST
Báo xấu
395
lượt xem 93
download
 
  Download Vui lòng tải xuống để xem tài liệu đầy đủ

Bài giảng Matlab ứng dụng - Phần II: Phương trình vi phân thường, trình bày các nội dung: bài toán giá trị đầu, ví dụ định luật 2 Newton, phương pháp Euler, phương pháp điểm giữa, phương pháp Runger - Kutta, bài toán giá trị biên, phương pháp vi phân cấp 2, phương trình vi phân cấp 4.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Bài giảng Matlab ứng dụng: Phần II - TS. Nguyễn Hoài Sơn

  1. PHÖÔNG TRÌNH VI PHAÂN THÖÔØNG NG TS Nguyễn Hoài Sơn
  2. NOÄI DUNG: Baøi toaùn giaù trò ñaàu : Ví duï ñònh luaät 2 Newton Phöông phaùp Euler Phöông phaùp ñieåm giöõa Phöông phaùp Runge-Kutta Baøi toaùn giaù trò bieân : Phöông trình vi phaân caáp 2 : TS Nguyễn Hoài Sơn Phöông trình vi phaân caáp 4
  3. Ví duï ñònh luaät 2 Newton 1.1 Ví duï ñònh luaät 2 Newton : r r F = ma Gia toác laø ñaïo haøm baäc 1 cuûa vaän toác theo thôøi gian, do ñoù : r dv r =a vaø dt r r dv F = dt m Minh hoïa: Ñònh luaät 2 Newton cho moät vaät noùng boû vaøo trong moâi tröôøng chaát loûng. Söï thay ñoåi nhieät TS Nguyễn Hoài Sơn ñoä theo thôøi gian cuûa vaät ñöôïc moâ taû bôûi phöông Ts trình vi phaân caân baèng naêng löôïng. dT mc = −Q dt
  4. Ví duï ñònh luaät 2 Newton Vôùi nhieät naêng do laøm laïnh: Q = hA(Ts − T∞ ) Giaû söû vaät lieäu coù tính caùch nhieät cao : => Ts = T dT dT hA mc = − hA(T − T∞ ) hoaëc =− (T − T∞ ) dt dt mc Ví duï 1: dy = −y y (0) = y 0 dt Phöông trình naøy coù theå tích phaân tröïc tieáp : TS Nguyễn Hoài Sơn y dy = − dt ln = −t C2 y ln y = -t + C y = C2e-t ln y –lnC2 = -t y = y0e−-t
  5. Ví duï ñònh luaät 2 Newton Tích phaân soá cuûa caùc phöông trình vi phaân Cho : dy y (0) = y0 = f (t , y ); dt Tìm keát quaû chính xaùc taïi giaù trò t baát kì : tj = t0 + jh Vôùi h laø böôùc thôøi gian. y Goïi: Keát quaû soá taïi t3 y( t ) = keát quaûchính xaùc f ( t0,y0) = ñoä doác ñoà thò taïi (t0,y0) TS Nguyễn Hoài Sơn y( tj )= keát quaû chính xaùc taïi tj yj = keát quaû gaàn ñuùng taïi tj Keát quaû chính xaùc y(t) y0 f(tj , yj ) = keát quaû gaàn ñuùng cuûa haøm veà phía phaûi taïi t t
  6. Phöông phaùp Euler Cho h = t1 – t0 vaø ñieàu kieän ban ñaàu, y = y(t0), tính : y1 = y 0 + hf (t 0 , y 0 ) y 2 = y1 + hf (t1 , y1 ) ... y j +1 = y j + hf (t j , y j ) Hoaëc y j = y j −1 + hf (t j −1 , y j −1 ) Ví duï 2: Söû duïng phöông phaùp Euler ñeå tính TS Nguyễn Hoài Sơn dy = t − 2y y(0) = 1 Keát quaû chính xaùc dt laø : 1 y j = y j −1 + hf (t j −1 , y j −1 ) y = [2t − 1 + 5e − 2t ] 4
  7. Phöông phaùp Euler f (t j −1 , y j −1 ) Euler C.xaùc Sai soá j tj yj=yj+hf(tj-1,yj-1) y(tj) yj-y(tj) 0 0.0 NaN (Ñk ban ñaàu) 1.000 1.000 0 1 0.2 0-(2)(1) = - 2.000 1.0 + (0.2)(-2.0) = 0.60 0.688 -0.0879 2 0.4 0.2 – (2)(0.6) = -1.000 0.6 + (0.2)(-1.0) = 0.40 0.512 -0.1117 3 0.6 0.4 – (2)(0.4) = -0.4 0.4 + (0.2)(- 0.4) = 0.32 0.427 -0.1065 So saùnh vôùi ñoà thò : Ñoái vôùi h ñaõ bieát, sai soá lôùn nhaát trong keát quaû soá ñoù laø sai soá rôøi raïc toaøn cuïc max ( ∑ ( y j − y (t j )) ) j TS Nguyễn Hoài Sơn
  8. Phöông phaùp Euler Ñaùnh giaù sai soá : Sai soá ñòa phöông taïi moãi böôùc laø: ej = yj – y(tj) vôùi y(tj) laø keát quaû chính xaùc taïi tj GDE = max( ej ) j = 1, … Giaûi baèng Matlab: h max(ej ) function [t,y] = odeEuler(diffeq,tn,h,y0)] 0.200 0.1117 0.100 0.0502 t = (0:h:tn)’; 0.050 0.0240 n = length(t); TS Nguyễn Hoài Sơn 0.025 0.0117 y = y0 + ones(n , 1); for j = 2 : n y(j) = y(j – 1) + h* feval(diffeq,t(j -1),y(j-1)); end >> rhs = inline(‘cos(t)’,’t’,’y’) ; >> [t,Y] = odeEuler(rhs,2*pi,0.01, 0) ; >> plot(t,Y,’o’) ;
  9. Phöông phaùp ñieåm giöõa Taêng möùc ñoä chính xaùc baèng caùch tính ñoä nghieâng 2 laàn trong moãi böôùc cuûa h: k1 = f (t j , y j ) Tính moät giaù trò cuûa y taïi ñieåm giöõa : f (t j , y j ) h y j +1 / 2 = y j + 2 Ñaùnh giaù laïi ñoä nghieâng h h k 2 = f (t j + , y j + k1 ) 2 2 Tính giaù trò cuoái cuøng cuûa y y j +1 = y j + hk 2 TS Nguyễn Hoài Sơn
  10. Phöông phaùp ñieåm giöõa yj töø phöông phaùp Euler ñaùnh giaù ñoä doác taïi t +0.5h j -1 keát quaû chính xaùc taïi y j -1 y +0.5hk 1 j- 1 y töø phöông phaùp ñieåm giöõa y j j -1 0.5h 0.5h t t j -1 j Giaûi baèng Matlab: function [t,y] = odeMidpt(diffeq,tn,h,y0)] t = (0:h:tn)’; >> rhs = inline(‘cos(t)’,’t’,’y’) ; n = length(t) ; TS Nguyễn Hoài Sơn >> [t,Y] = odeEuler(rhs,2*pi,0.01, 0) ; y = y0 + ones(n , 1) ; >> plot(t,Y,’o’) ; h2= h /2 ; for j = 2 : n k1 = feval(diffeq,t(j -1),y(j-1)) ; k2 = feval(diffeq,t(j -1)+h2,y(j-1)+h2*k1) ; y(j) = y(j – 1) + h* k2 ; end
  11. Phöông phaùp ñieåm giöõa So saùnh phöông phaùp Midpoint vôùi phöông phaùp Euler Giaûi: dy = −y y(0) = 1 0 ≤ t ≤ 1 dt Keát quaû chính xaùc laø : y = e-t h flopeE errE flopeH errH 0.20000 31 4.02e-02 57 2.86e-03 0.10000 61 1.92e-02 112 6.62e-04 0.05000 121 9.39e-03 222 1.59e-04 0.02500 241 4.66e-03 442 3.90e-05 0.01250 481 2.31e-03 882 9.67e-06 TS Nguyễn Hoài Sơn 0.00625 961 1.15e-03 1762 2.41e-06
  12. Phöông phaùp Runge-Kutta Tính ñoä doác ôû 4 vò trí öùng vôùi moãi böôùc laëp: k1 = f (t j , y j ) h h k 2 = f (t j + , y j + k1 ) 2 2 h h k 3 = f (t j + , y j + k2 ) 2 2 k 4 = f (t j + h, y j + hk 3 ) Ta tính ñöôïc yj+1 ⎛k k k k ⎞ y j +1 = y j + h⎜ 1 + 2 + 3 + 4 ⎟ TS Nguyễn Hoài Sơn ⎝6 3 3 6⎠
  13. Phöông phaùp Runge-Kutta Giaûi baèng Matlab >> rhs = inline(‘cos(t)’,’t’,’y’) ; function [t,y] = odeRK4(diffeq,tn,h,y0)] >> [t,Y] = odeRK4(rhs,2*pi,0.01, 0) ; t = (0:h:tn)’; >> plot(t,Y,’o’) ; n = length(t) ; Haøm thö vieän Matlab y = y0 + ones(n , 1) ; h2= h /2 ; h3= h /3 ; h6= h /6 ; [t,Y] = ode45(diffep,tn,y0) for j = 2 : n >> rhs = inline(‘cos(t)’,’t’,’y’) ; k1 = feval(diffeq, t(j -1), y(j-1)) ; >> [t,Y] = ode45(rhs,[0 2*pi], 0) ; k2 = feval(diffeq , t(j -1)+h2, y(j-1)+h2*k1 ) ; >> plot(t,Y,’r’,’linewidth’,2) ; k3 = feval(diffeq , t(j -1)+h2, y(j-1)+h2*k2 ) ; k4 = feval(diffeq , t(j -1)+h , y(j-1)+h*k3) ; y(j) = y(j – 1) + h6* (k1+k4) + h3*(k2+k3); TS Nguyễn Hoài Sơn end
  14. Phöông phaùp Runge-Kutta So saùnh Euler, Midpoint vaø RK4: Giaûi: dy 0 ≤ t ≤1 = −y y(0) = 1 dt h flope errE flopeM errM flope4 err4 E 0.20000 31 4.02e-02 57 2.86e-3 129 5.80e-6 0.10000 61 1.92e-02 112 6.62e-4 254 3.33e-7 0.05000 121 9.39e-03 222 1.59e-4 504 2.00e-8 0.02500 241 4.66e-03 442 3.90e-5 1004 1.22e-9 0.01250 481 2.31e-03 882 9.67e-6 2004 7.56e-11 0.00625 961 1.15e-03 1762 2.41e-6 4004 4.70e-12 Söû duïng haøm cuûa Matlab: TS Nguyễn Hoài Sơn Söû duïng ode45 Cuù phaùp : [t,Y] = ode45(diffep,tn,y0) [t,Y] = ode45(diffep,[t0 tn],y0) [t,Y] = ode45(diffep,[t0 tn],y0,options) [t,Y] = ode45(diffep,[t0 tn],y0,options,arg1,arg2,…)
  15. Phöông phaùp Runge-Kutta Ví duï dy y(0) = 0 = cos( t ) dt >> rhs = inline(‘cos(t)’,’t’,’y’) ; >> [t,Y] = ode45(rhs,[0 2*pi], 0) ; >> plot(t,Y,’o’) ; TS Nguyễn Hoài Sơn
  16. Baøi toaùn giaù trò bieân : Phöông trình vi phaân caáp 2 : ÖÙng duïng cho caùc baøi toaùn veà thanh , truyeàn nhieät ,vv… Daïng : ay’’(x)+by’(x)+cy(x)=f(x) 0 < x < 1 (7.10) Ñieàu kieän bieân : a/ y(x=0) = y0 * y(x= L) = y L * b/ y’(x=0) = yo' * y(x= L) = yL * c/ y(x=0) = y0 * y’(x=L) = y 'L * TS Nguyễn Hoài Sơn Xaáp xæ (7.10) baèng löôùi ñeàu sai phaân trung taâm :ho= Δx yi +1 − yi −1 1 yi’ = + 0(h2) vôùi O(h2) = - h2fi’’’ (7.11) 2h 6 y i +1 − 2 f i + y i −1 1 yi’’ = 2 + 0(h2) vôùi 0(h2) = - h2fi’’’ (7.12) h 12
  17. (7.10), (7.11) vaø(7.12) cho ta phöông trình sai phaân ⎡ y i +1 − 2 y i + y i −1 ⎤ ⎡ y i +1 − y i −1 ⎤ a⎢ h2 ⎥ + b⎢ 2h ⎥ + cy i = f ( x) (7.13) ⎣ ⎦ ⎣ ⎦ (2a+ bh)yi+1+(2ch2 - 4a)yi + (2a - bh)yi-1 = 2h2f(x) (7.14) i=1 => (2a + bh)y2 + (2ch2 - 4a)y1 + (2a - bh)yo = 2h2f(x) i=2 => (2a + bh)y3 + (2ch2 - 4a)y2 + (2a - bh)y1 = 2h2f(x) i=3 => (2a + bh)y4 + (2ch2 - 4a)y3 + (2a - bh)y2 = 2h2f(x) TS Nguyễn Hoài Sơn i=4 => (2a + bh)y5 + (2ch2 - 4a)y4 + (2a - bh)y3 = 2h2f(x) i=5 => (2a + bh)y6 + (2ch2 - 4a)y5 + (2a - bh)y4 = 2h2f(x)
  18. Ñaët : A=2a + bh B=2ch2 - 4a C=2a – bh Ñöa heä 5 phöông trình treân veà daïng ma traän : a/ By1 +Ay2 = 2h2f(x)-Cy0* Cy1+By2+Ay3 = 2h2f(x) Cy2+By3+Ay4 = 2h2f(x) Cy3+By4+Ay5 = 2h2f(x) Cy4+By5 = 2h2f(x)-AyL* Hay dạng ma trận : ⎡B A 0 0 0 ⎤ ⎛ y1 ⎞ ⎡2h 2 f(x) - Cy 0 * ⎤ TS Nguyễn Hoài Sơn ⎢C B A ⎥ ⎜ y2 ⎟ 0 0⎥ ⎜ ⎟ ⎢ 2 ⎥ ⎢ ⎢2h f(x) ⎥ ⎢0 C B A 0⎥ ⎜ y3 ⎟ = ⎢2h 2 f(x) ⎥ ⎢ ⎥ ⎜ ⎟ ⎢ ⎢2h 2 f(x) ⎥ ⎥ ⎢0 0 C B A⎥ ⎜ y4 ⎟ ⎢ 2 ⎥ ⎢0 ⎣ ⎥ ⎜ y5 ⎟ 0 0 C B⎦ ⎝ ⎠ ⎢2h f(x) - Ay L *⎥ ⎣ ⎦
  19. y2 = y0 b/ y'1 = 2h = y'0 * (Biết) ⇒ y0=y2 -2hy0’* ⎡B A + C 0 0 0 ⎤ ⎛ y1 ⎞ ⎡2hCy 0 '* + 2h 2 f(x)⎤ ⎢C B ⎜ ⎟ ⎢ 2 ⎥ ⎢ A 0 0 ⎥ ⎜ y2 ⎟ ⎥ ⎢2h f(x) ⎥ (5) ⇒ ⎢0 ⎢ C ⎥⎜ ⎟ = B A 0 ⎥ ⎜ y3 ⎟ ⎢2h 2 f(x) ⎢ ⎥ ⎥ ⎢0 0 C B A⎥ ⎜ y4 ⎟ ⎢2h f(x) 2 ⎥ ⎢0 0 0 C B⎥ ⎜ y5 ⎟ ⎢ 2 ⎥ ⎣ ⎦⎝ ⎠ ⎢2h f(x) - Ay L * ⎥ ⎣ ⎦ y6 − y4 c/ y'5 = 2h = y L '* ⇒ y 6 = y 4 + 2hy L '* TS Nguyễn Hoài Sơn ⎡B A 0 0 0⎤ ⎛ y1 ⎞ ⎡2h 2 f(x) - Cy 0 * ⎤ ⎢C ⎜ ⎟ ⎢ 2 ⎥ ⎢ B A 0 0⎥ ⎥ ⎜ y2 ⎟ ⎢2h f(x) ⎥ ⎢0 ⎢ C B A 0⎥ ⎥ ⎜ y3 ⎟ ⎜ ⎟ = ⎢2h 2 f(x) ⎢ ⎥ ⎥ ⎢0 0 C B A⎥ ⎜ y4 ⎟ ⎢2h 2 f(x) ⎥ ⎢0 0 0 A + C B⎥ ⎜ y5 ⎟ ⎢ 2 ⎥ ⎣ ⎦ ⎝ ⎠ ⎢2h f(x) - 2hAy L '*⎥ ⎣ ⎦
  20. • I. Phöông phaùp soâ: Chia ñoâi khoaûng, Newton-Raphson, Daây cung 1. Chia ñoâi khoaûng (Bisection Method) (1) : a, b, tol Matlab code. (2) : k = 1,2,.. clear all a+b clc (3) : xm = a=3;b=4;tol=0.0001 2 for k=1:10 ( 4) : f ( x m ) ≠ 0 x=(a+b)/2; if sign(f(x))==sign(f(a)) (5) : f (a). f ( xm ) < 0 a=x; else ( 6) : b = x m b=x; end (7) : a = x m if abs(f(x)>tol) TS Nguyễn Hoài Sơn break (8) : f ( xm ) ≤ tol end a+b end (9) : xm = function gg=f(x) 2 gg=x-x.^(1/3)-2; (10) : in x
ADSENSE

CÓ THỂ BẠN MUỐN DOWNLOAD

LV.15: Bộ Đồ Án Tốt Nghiệp Chuyên Ngành Cơ Khí
65 tài liệu
2431 lượt tải
THÔNG TIN
TRỢ GIÚP
HỖ TRỢ KHÁCH HÀNG
Theo dõi chúng tôi

Chịu trách nhiệm nội dung:

Nguyễn Công Hà - Giám đốc Công ty TNHH TÀI LIỆU TRỰC TUYẾN VI NA

LIÊN HỆ

Địa chỉ: P402, 54A Nơ Trang Long, Phường 14, Q.Bình Thạnh, TP.HCM

Hotline: 093 303 0098

Email: support@tailieu.vn

Giấy phép Mạng Xã Hội số: 670/GP-BTTTT cấp ngày 30/11/2015 Copyright © 2022-2032 TaiLieu.VN. All rights reserved.

 

Đồng bộ tài khoản
3=>0