Giải phương trình vi phân bằng phương pháp Runge-Kuta 4

Chia sẻ: Nguyễn Phước Lộc | Ngày: | Loại File: PDF | Số trang:9

Thêm vào BST

Báo xấu

1.794
lượt xem 174
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Đây là báo cáo bài tập lớn bộ môn phương pháp tính khoa điện - điện tử nghiên cứu về đề tài " giải phương trình vi phân bằng phương pháp Runge-Kuta 4 gửi đến các bạn độc giả tham khảo.

Chủ đề:

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

Lưu

Nội dung Text: Giải phương trình vi phân bằng phương pháp Runge-Kuta 4

ĐẠI HỌC QUỐC GIA THÀNH PHỐ HỒ CHÍ MINH TRƯỜNG ĐẠI HỌC BÁCH KHOA *** KHOA ĐIỆN – ĐIỆN TỬ Giáo viên hướng dẫn: Lê Thị Quỳnh Hà Các thành viên trong nhóm: 40901457: Nguyễn Phước Lộc 40902767: Võ Nhựt Tiến 40902741: Cù Văn Tiến 40903057: Huỳnh Trung Trực 40902907: Nguyễn Kim Triển 40902935: Phan Đức Trí
Đề tài: Viết chương trình giải phương trình bằng phương pháp Runge-Kutta bậc 4. Vẽ đồ thị hàm nhận được.
CÔ SÔÛ LYÙ THUYEÁT BAØI TOAÙN CAUCHY M t ph ng trình vi phân c p 1 có th vi t d i d ng gi i c yf (x,y) mà ta có th tìm c hàm y t o hàm c a nó. T n t i vô s nghi m tho mãn ph ng trình trên. M i nghi m ph thu c vào m t h ng s tu ý. Khi cho tr c giá tr ban u c a y là yo t i giá tr u xo ta nh n c m t nghi m riêng c a ph ng trình. Bài toán Cauchy (hay bài toán có i u ki n u) tóm l i nh sau: cho x sao cho b x a, tìm y(x) tho mãn i u ki n: y ( x) f ( x , y ) (1) y(a ) Ng i ta ch ng minh r ng bài toán này có m t nghi m duy nh t n u f tho mãn i u ki n Lipschitz: f( x , y 1 ) f( x , y 2 ) L y 1 y 2 v i L là m t h ng s d ng. Ng i ta c ng ch ng minh r ng n u f y ( o hàm c a f theo y ) là liên t c và b ch n thì f tho mãn i u ki n Lipschitz. M t cách t ng quát h n, ng i ta nh ngh a h ph ng trình b c 1: y 1 f1 ( x , y 1 , y 2 ,..., y n ) y2 f2 ( x , y 1 , y 2 ,..., y n ) yn fn ( x , y 1 , y 2 ,..., y n ) Ta ph i tìm nghi m y1, y2,..., yn sao cho: Y ( x) f ( x , X ) Y(a ) v i: y1 f1 y1 y2 f2 y2 F Y Y .. .. .. .. .. .. fn yn yn
N u ph ng trình vi phân có b c cao h n (n), nghi m s ph thu c vào n h ng s tu ý. nh n c m t nghi m riêng, ta ph i cho n i u ki n u. Bài toán s có giá tr u n u v i giá tr xo ã cho ta cho y(xo), y (xo), y (xo),.... M t ph ng trình vi phân b c n có th a v thành m t h ph ng trình vi phân c p 1. Ví d n u ta có ph ng trình vi phân c p 2: y f( x , y , y ) y(a ) , y (a ) t u = y và v = y ta nh n c h ph ng trình vi phân c p 1: Khi uv v g( x , u , v ) v i i u ki n u: u(a) = và v(a) = Các ph ng pháp gi i ph ng trình vi phân c trình bày trong ch ng này là các ph ng pháp r i r c: o n [a, b] c chia thành n o n nh b ng nhau c g i là các b c tích phân h = ( b a) / n. §2. PH NG PHÁP EULER
PHÖÔNG PHAÙP RUNGE KUTTA ch a chính xác i v i các bài toán th c t . Xét bài toán Cauchy (1). Gi s ta ã tìm c giá tr g n úng yi c a y(xi) và mu n tính yi+1 c a y(xi+1). Tr c h t ta vi t công th c Taylor: h2 h m (m) h m 1 ( m 1) y( x i 1 ) y( x i ) hy ( x i ) y (x i ) y (x i ) y (c) (11) 2 m! m! v i c (xi, xi+1) và: y ( x i ) f x i , y( x i ) dk 1 ( k) y (xi ) f x i , y(x i ) dx k 1 Ta vi t l i (11) d i d ng: h2 h m (m) h m 1 ( m 1) y i 1 y i hy ( x i ) y (x i ) y (x i ) y ( c) (12) 2 m! m! Ta ã kéo dài khai tri n Taylor k t qu chính xác h n. tính y i, y i v.v. ta có th dùng ph ng pháp Runge Kutta b ng cách t: y i 1 y i r1 k (1i ) r2 k (2i ) r3 k (3i ) r4 k (4i ) (13) trong ó: k (1i ) hf( x i , y i ) k (2i ) k (1i ) ) hf( x i ah , y i (14) k (3i ) k (1i ) k (2i ) ) hf( x i bh , y i ....... và ta c n xác nh các h s a, b,..; , , ,...; r1, r2,.. sao cho v ph i c a (13) khác v i v ph i c a (12) m t vô cùng bé c p cao nh t có th có i v i h. Khi dùng công th c Runge Kutta b c hai ta có: k (1i ) hf( x i , y i ) (15) k (2i ) hf( x i ah , y i k (1i ) ) y i 1 y i r1 k 1i ) r2 k (2i ) ( và (16) Ta có: y (x) = f[x,y(x)] y ( x) fx x , y( x) fy x , y( x) ................ Do ó v ph i c a (12) là: h2 hf( x i , y i ) fx (x i , y i ) fy ( x i , y i ) y ( x) (17) 2 M t khác theo (15) và theo công th c Taylor ta có:
k 1i ) ( hf ( x i , y i ) hy i k (2i ) k (1i ) fy ( x i , y i ) h[f( x i , y i ) ahfx ( x i , y i ) ] Do ó v ph i c a (16) là: h(r1 r2 )f( x i , y i ) h 2 [ar2 fx ( x i , y i ) r2 y i fy ( x i , y i )] (18) Bây gi cho (17) và (18) khác nhau m t vô cùng bé c p O(h3) ta tìm c các h s ch a bi t khi cân b ng các s h ng ch a h và ch a h2: r1 + r2 = 1 a.r1 = 1/ 2 .r2 = 1 Nh v y: = a, r1 = (2a 1)/ 2a, r2 = 1/ 2a v i a c ch n b t kì. N u a = 1 / 2 thì r1 = 0 và r2 = 1. Lúc này ta nh n c công th c Euler. N u a=1 thì r1 = 1 / 2 và r2 = 1/2. Lúc này ta nh n c công th c Euler c i ti n. M t cách t ng t chúng ta nh n c công th c Runge Kutta b c 4. Công th c này hay c dùng trong tính toán th c t : k1 = h.f(xi, yi) k2 = h.f(xi+h/ 2, yi + k1/ 2) k3 = h.f(xi+h/ 2, yi + k2/ 2) k4 = h.f(xi+h, yi + k3) yi+1 = yi + (k1 + 2k2 + 2k3 + k4) / 6 Ta xây d ng hàm rungekutta() th c hi n công th c Runge Kutta b c 4:
Ví dụ : Dùng công thức Runge-Kutta tìm nghiệm gần đúng của bài toán Cauchy y’ = y – x2 +1, 0≤ x ≤1 y(0) = 0.5 với n = 5 Tính sai số biết nghiệm chính xác là : y(x) = (x+1)2 – 0.5ex Giải:  Ta có h = 0.2 x0 = 0, x1 = 0.2, x2 = 0.4, x3 = 0.6, x4 = 0.8, x5 = 1
Xây dựng hàm rk4 trong matlab để giải phương trình vi phân theo phương pháp trên. function[x,y]=rk4(f,x0,x1,y0,h) if nargin [X,Y]=rk4(inline('y-x^2+1','x','y'),0,1,0.5) Định nghĩa hàm f trong của sổ Script rồi lưu thành file f.m  function [dy] = f(a,b) dy = b-a^2 +1; end >> [X,Y]=rk4(@f,0,1,0.5) Vẽ đồ thị từ các giá trị xk,yk (giá trị trả về của hàm rk4):  plot(X,Y,'rd-', 'LineWidth',3) xlabel('Truc X') ylabel('Truc Y') title('DO THI CUA HAM DA CHO') grid on
Kết quả cửa sổ Command Windows: >> [X,Y]=rk4(inline('y-x^2+1','x','y'),0,1,0.5) Nhap so doan chia n = 5 n= 5 X= 0 0.2000 0.4000 0.6000 0.8000 1.0000 Y= 0.5000 0.8293 1.2141 1.6489 2.1272 2.6408 >> ve >>