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à

Nguyễn Phước Lộc Võ Nhựt Tiến Cù Văn Tiến Huỳnh Trung Trực Nguyễn Kim Triển Phan Đức Trí

Các thành viên trong nhóm: 40901457: 40902767: 40902741: 40903057: 40902907: 40902935:

Đề 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.

ng trình vi phân c p 1 có th vi t d i d ng gi i

c

c hàm y t

M t ph (x,y) mà ta có th tìm

o hàm c a nó. T n t i vô s ng trình trên. M i nghi m ph thu c vào m t h ng s

c giá tr ban

u c a y là yo t i giá tr

u xo ta nh n

c ng trình. Bài toán Cauchy (hay bài toán có i u a, tìm y(x) tho mãn i u

x

u) tóm l i nh sau: cho x sao cho b

yf nghi m tho mãn ph tu ý. Khi cho tr m t nghi m riêng c a ph ki n ki n:

)x(y )y,x(f (1) )a(y

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

y

)y,x(f 2

)y,x(f 1

1

2

tho mãn i u ki n Lipschitz: yL

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.

i ta nh ngh a h ph ng trình b c 1:

,...,

1

1

,...,

y

y,y,x(f 1 2 y,y,x(f 2

1

2

)y n )y n

2

y

,...,

y,y,x(f n

1

)y n

n

2 Ta ph i tìm nghi m y1, y2,..., yn sao cho:

)X,x(f

)x(Y )a(Y

v i:

M t cách t ng quát h n, ng y

1

1

2

2

y

F Y Y

y .. .. f 1 f 2 .. .. y y .. ..

n

n

y y f n

N u ph

ng trình vi phân có b c cao h n (n), nghi m s ph thu c vào u.

c m t nghi m riêng, ta ph i cho n i u ki n

nh n

u n u v i giá tr xo ã cho ta cho y(xo), y (xo), y (xo),....

ng trình vi phân b c n có th

a v thành m t h ph

ng

n h ng s tu ý. Bài toán s có giá tr M t ph

ng trình vi phân c p 2:

trình vi phân c p 1. Ví d n u ta có ph

y

)y,y,x(f

)a(y

)a(y,

c h ph

ng trình vi phân c p 1:

Khi

t u = y và v = y ta nh n u

v

v

)v,u,x(g

và v(a) =

v i

u: u(a) = ng pháp gi i ph

ng trình vi phân

i u ki n Các ph ng này là các ph

ng pháp r i r c: o n [a, b]

c trình bày trong c chia thành n o n

c g i là các b

c tích phân h = ( b a) / n.

ch nh b ng nhau

§2. PH

NG PHÁP EULER

ch a chính xác i v i các bài toán th c t .

c giá tr g n úng yi c a

2

m

1m

)m(

)1m(

Xét bài toán Cauchy (1). Gi s y(xi) và mu n tính yi+1 c a y(xi+1). Tr ta ã tìm c h t ta vi t công th c Taylor:

1i

y y )c( (11) )x(yh)x(y)x(y i i )x(y i )x( i h 2 h !m h !m

v i c

(xi, xi+1) và: )x(y i

)x(y,xf i 1k

)k(

i

1k

i d dx

2

1m

m

)m(

)1m(

)x(y,xf i )x(y i

1i

i

y y y y )c( (12) )x(yh i )x(y i )x( i Ta vi t l i (11) d i d ng: h 2 h !m

h !m k t qu chính xác h n. tính y i, y i v.v. ta

t:

y

kr 4

kr 2

kr 3

kr 1

)i( 1

)i( 4

)i( 2

i

)y,x(hf i i

1i trong ó: )i( k 1

x(hf

y,ah

k

i

i

)i( )k 1

)i( 2

(14)

x(hf

y,bh

k

k

i

i

)i( 1

)i( )k 2

)i( 3 .......

nh các h s a, b,..;

,

,

,...; r1, r2,.. sao cho v ph i c a (13)

i v i h.

và ta c n xác khác v i v ph i c a (12) m t vô cùng bé c p cao nh t có th có Khi dùng công th c Runge Kutta b c hai ta có:

k

)i( 1

)y,x(hf i i

(15)

k

x(hf

)i( 2

)i( )k 1

y

y

(16)

y,ah )i( 1

i kr 2

i kr 1

)i( 2

1i

i

và Ta có:

y (x) = f[x,y(x)] )x(y

)x(y,xf y

)x(y,xf x ................

Do ó v ph i c a (12) là:

2

(17)

)y,x(hf i i

)y,x(f x i i

)x(y)y,x(f y i

i

h 2

M t khác theo (15) và theo công th c Taylor ta có:

Ta ã kéo dài khai tri n Taylor có th dùng ph y (13) ng pháp Runge Kutta b ng cách )i( 3

i

]

)i( 1 )i( 2

)y,x(fah i i

x

)y,x(fk i

y

)i( 1

i

(18)

Do ó v ph i c a (16) là: 2 )y,x(far[h)y,x(f)r 2 i

k k yh)y,x(hf i i )y,x(f[h i i

r(h 1

i

x2

yi

i

)]y,x(fyr 2

i

i

c các

i Bây gi cho (17) và (18) khác nhau m t vô cùng bé c p O(h3) ta tìm 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

= a, r1 = (2a 1)/ 2a, r2 = 1/ 2a v i a

c ch n b t kì.

c công th c Euler. N u

c công th c Euler c i ti n.

Nh v y: N u a = 1 / 2 thì r1 = 0 và r2 = 1. Lúc này ta nh n a=1 thì r1 = 1 / 2 và r2 = 1/2. Lúc này ta nh n ng t

chúng ta nh n

c công th c Runge Kutta b c 4.

M t cách t

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

0≤ x ≤1

y’ = y – x2 +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<4, error('Vui long nhap du doi so !! '), end; if nargin<5, m = input('Nhap so doan chia n = ');,h=(x1-x0)/m; end; x=[]; x(1)=[x0]; n=(x1-x0)/h; for i=1:n, x(i+1)=x(i)+h; end; y=[]; y(1)=[y0]; for i=1:n K1=h*f(x(i),y(i)); K2=h*f(x(i)+h/2,y(i)+K1/2); K3=h*f(x(i)+h/2,y(i)+K2/2); K4=h*f(x(i)+h,y(i)+K3); y(i+1)=y(i)+(K1+2*K2+2*K3+K4)/6; end;

Hàm rk4 sẽ nhận vào 4 đến 5 đối số. Nếu ta nhập số đối số bé hơn 4 thì chương trình sẽ báo lỗi để nhắc nhỡ. Nếu ta nhập số đối số bằng 4 thì chương trình sẽ yêu cầu ta nhập số đoạn chia n. Để giải phương trình vi phân dùng hàm trên trước hết ta phải định nghĩa hàm f.

Ví dụ: Ta giải lại phương trình: y’ = y – x2 +1, 0≤ x ≤1, y(0) = 0.5, với n = 5



>> [X,Y]=rk4(inline('y-x^2+1','x','y'),0,1,0.5)

Định nghĩa trực tiếp hàm f từ của sổ Command Windows:

 Đị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

>>