ƯƠ
Ả
ƯƠ
CH
NG 7: GI
I PH
NG TRÌNH VI PHÂN
§1. BÀI TOÁN CAUCHY
ể ế ạ ng trình vi phân c p 1 có th vi ộ M t ph
ượ ừ ạ
ả ệ
ấ c hàm y t ỗ ủ ướ ạ ị ầ c giá tr ban đ u c a y là y
ươ ệ c m t nghi m riêng c a ph (cid:0) ủ ư ầ ạ ả ượ ả ướ c i d ng gi t d i đ ố ạ ồ ủ i vô s đ o hàm c a nó. T n t ộ ằ ộ ụ ng trình trên. M i nghi m ph thu c vào m t h ng o ta nh nậ ị ầ i giá tr đ u x o t ng trình. Bài toán Cauchy (hay bài toán có a, tìm y(x) tho mãn i nh sau: cho x sao cho b x (cid:0)
ươ y(cid:0) =f(x,y) mà ta có th tìm đ ể ươ ệ nghi m tho mãn ph ố ỳ s tu ý. Khi cho tr ượ đ ề đi u ki n đ u) tóm l ề đi u ki n: (cid:0) (cid:0) (cid:0) ộ ệ ệ )x(y )y,x(f (cid:0) (1) (cid:0) (cid:0) (cid:0)
ấ ế ệ ộ )a(y ườ Ng
ả ằ tho mãn đi u ki n Lipschitz:
2
1
(cid:0) (cid:0) (cid:0) yL y ứ ệ )y,x(f 2
ớ 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,x(f 1 ộ ằ v i L là m t h ng s d ườ ủ ế ằ ạ (cid:0) y ( đ o hàm c a f theo y ) là liên
i ta cũng ch ng minh r ng n u f ặ ố ươ ứ ả ụ t c và b ch n thì f tho mãn đi u ki n Lipschitz.
1
2
1
ệ ị ườ ươ ậ Ng ị ộ ổ ơ ệ i ta đ nh nghĩa h ph ng trình b c 1: (cid:0) (cid:0) ề M t cách t ng quát h n, ng y ,...,
1
2
2
(cid:0) (cid:0) ,..., y y,y,x(f 1 y,y,x(f 2 )y n )y n
(cid:0) (cid:0) (cid:0) (cid:0)
n
(cid:0) (cid:0) )y n
1 ệ )X,x(f
y ả y,y,x(f n 2 Ta ph i tìm nghi m y ,..., 1, y2,..., yn sao cho: (cid:0) (cid:0) (cid:0) )x(Y (cid:0) (cid:0) (cid:0) )a(Y (cid:0)
1
1
2
2
v i:ớ (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) y y (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) y y (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) .. f 1 f 2 .. Y F Y .. (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) .. .. .. (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) f n
ụ ệ ậ
ẽ ả ể ậ
ượ ị ầ ộ ế ẽ ệ ớ ị y n ơ ộ ng trình vi phân có b c cao h n (n), nghi m s ph thu c ề ệ c m t nghi m riêng, ta ph i cho n đi u o đã cho ta cho y(x o),
168
y n ế ươ N u ph ằ ố ỳ vào n h ng s tu ý. Đ nh n đ ầ ki n đ u. Bài toán s có giá tr đ u n u v i giá tr x y(cid:0) (xo), y(cid:0) (xo),....
ệ ậ ươ ể ư ộ M t ph ng trình vi phân b c n có th đ a v thành m t h ph ng
ấ ụ ế ươ ộ ề ng trình vi phân c p 2: trình vi phân c p 1. Ví d n u ta có ph (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) ươ ấ )y,y,x(f y (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) )a(y )a(y, (cid:0)
ượ ươ ấ (cid:0) ta nh n đ ậ ệ c h ph ng trình vi phân c p 1: (cid:0) (cid:0) (cid:0) Khi đ t u = y và v = y v ặ u (cid:0) (cid:0) (cid:0) v )v,u,x(g (cid:0)
(cid:0) ớ v i đi u ki n đ u: u(a) =
và v(a) = (cid:0) ươ ả i ph
ờ ạ ng trình vi phân đ ượ ạ ầ ề ệ ươ ng pháp gi 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
ươ ỏ ằ ươ ọ ượ ướ ch nh b ng nhau đ c g i là các "b c" tích phân h = ( b a) / n.
Ả Ế
ƯƠ ươ NG PHÁP EULER VÀ EULER C I TI N ng trình vi phân: Gi (cid:0) (cid:0) (cid:0) §2. PH s ta có ph )y,x(f (cid:0) (1) (cid:0) (cid:0) ả ử )x(y )a(y (cid:0)
ở ệ ạ ể ầ o,x ] thành n ph n b i các đi m
ủ ầ và c n tìm nghi m c a nó. Ta chia đo n [x chia:
xo < x1 < x2 <...< xn = x
3
i ta có: x(
1i
1i
1i
1i
i
ứ ể ộ Theo công th c khai tri n Taylor m t hàm lân c n x ậ 2 (cid:0) (cid:0) (cid:0) (cid:0) x( )x i )x i (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) x(y ) x( )x(y i )x(y i )x(y i )x(y)x i 2 6
i+1 xi) khá bé thì ta có th b qua các s h ng (x
i+1 xi)2 và các s h ng ố ạ
ể ỏ ố ạ
ườ ợ ế N u (x ậ b c cao y(xi+1) = y(xi) + (xi+1 xi)y(cid:0) (xi) Tr y
ề ố ng h p các m c cách đ u: (xi1 xi) = h = (x xo)/ n
i+1
i
ậ ượ ứ ơ thì ta nh n đ y yi+1 = y i + hf(xi, yi) y ả c công th c Euler đ n gi n: (2) ả ế ề ặ
x ế ộ ả ế ứ x ướ ả ụ i đ nh lí Lagrange:
ể ỏ c h càng nh . ể ế c h t ta i x ả ử s f(x) là hàm liên t c trong [a, b] và kh vi trong i+1 (cid:0) (a ,b) đ choể : (cid:0) ọ ấ V m t hình h c ta th y (1) cho k t qu càng ướ chính xác n u b ể Đ tăng đ chính xác ta có th dùng công th c Euler c i ti n. Tr ắ ạ ị gi nh c l ấ (a,b) thì có ít nh t m t đi m c )b(f ộ )a(f (cid:0) (cid:0) )c(f (cid:0) ab
i
i
169
Theo đ nh lí Lagrange ta có : (cid:0) (cid:0) (cid:0) )x(y i (cid:0) ị )x(y 1i ớ ư ậ Nh v y v i c ))c(y,c(hf ể (xi, xi+1) ta có th thay :
(cid:0) (cid:0) (cid:0) (cid:0) (cid:0) ))c(y,c(f x(f y,
(cid:0))
1i
i
i
1i
)y,x(f i i
ả ế ừ 1 2 ứ T đó ta có công th c Euler c i ti n :
(cid:0))
1i
i
1i
1i
(cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) y y x(f y, (3) )y,x(f i i h 2 ị
ằ
ế t. Do đó khi đã bi ế ng trình đ i s tuy n tính (3). Ta th ỉ ầ ướ ủ
)s( 1iy (cid:0)
ượ ng gi ặ c h t ch n x p x đ u tiên c a phép l p ể ươ ế c theo ph ư i+1 ch a bi ạ ố ọ ấ ng pháp Euler sau đó dùng (3) đ tính các ế i ta ph i tìm y ả t y i+1 ằ ả ườ i (3) b ng )0( 1iy (cid:0) chính là ,
i
(cid:0) (cid:0) (cid:0) y ứ Trong công th c này giá tr y ươ ả i ph b ng cách gi ư ặ cách l p nh sau: tr giá tr yị i+1 tính đ ụ ể c th là: )0( y 1i )y,x(hf i i
(cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) y y x(f y,
(cid:0))
)s( 1i
i
1i
)1s( 1i
)s(
)1s(
iy
)y,x(f i i (cid:0) ủ ầ iy đ g n
ươ ươ ươ ư h 2 ế ả Quá trình tính k t thúc khi Ch ng trình gi i ph ng trình vi phân theo ph ng pháp Euler nh sau:
ươ Ch ng trình 71
/ / p p_Euler;
#include
float f(float x,float y) {
float a=x+y; return(a);
}
void main() { int i,n;
float a,b,t,z,h,x0,y0,c1,c2; float x[100],y[100];
170
clrscr(); printf("Cho can duoi a = "); scanf("%f",&a); printf("Cho can tren b = "); scanf("%f",&b); printf("Cho so buoc tinh n = "); scanf("%d",&n);
printf("Cho so kien x0 = "); scanf("%f",&x0); printf("Cho so kien y0 = "); scanf("%f",&y0); printf("\n"); printf("Bang ket qua \n"); printf("\n"); printf("Phuong phap Euler \n"); h=(ba)/n; x[1]=x0; y[1]=y0; printf(" x y"); printf("\n"); for (i=1;i<=n+1;i++) {
x[i+1]=x[i]+h; y[i+1]=y[i]+h*f(x[i],y[i]); printf("%3.2f%16.3f",x[i],y[i]); printf("\n");
} printf("\n"); getch(); printf("Phuong phap Euler cai tien \ n"); printf(" x y"); printf("\n"); for (i=1;i<=n+1;i++) {
x[i+1]=x[i]+h; c1=h*f(x[i],y[i]); c2=h*f(x[i]+h,y[i]+c1); y[i+1]=y[i]+(c1+c2)/2; printf("%3.2f%15.5f",x[i],y[i]); printf("\n");
} getch();
}
o = 0, y o = 0,
ươ ề ệ ầ ớ V i ph ng trình cho trong function và đi u ki n đ u x
ể ệ ạ ớ nghi m trong đo n [0, 1] v i 10 đi m chia là:
ả ế
171
x 0.0 0.1 y(Euler) 0.00 0.00 y(Euler c i ti n) 0.00 0.01
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.01 0.03 0.06 0.11 0.17 0.25 0.34 0.46 0.59 0.02 0.05 0.09 0.15 0.22 0.31 0.42 0.56 0.71
ƯƠ §3. PH NG PHÁP RUNGE KUTTA
i c aủ
Xét bài toán Cauchy (1). Gi ị ầ c giá tr g n đúng y
1m
)1m(
)m(
1i
i
i+1). Tr 2 h 2
i
1k
)k(
(cid:0))x(y,xf
i
i
1k
ố ả ử s ta đã tìm đ ế ế ướ y(xi) và mu n tính y ủ i+1 c a y(x (cid:0) c h t ta vi m (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (11) x(y ) y y )c( )x(yh)x(y i )x(y i )x( i ượ ứ t công th c Taylor: h !m h !m (cid:0) v i c ớ (cid:0) (cid:0) (cid:0) (xi, xi+1) và: )x(y i (cid:0) (cid:0) (cid:0) y (cid:0) )x( i
(cid:0))x(y,xf i d dx ướ ạ
m
1m
2
)m(
)1m(
i
1i
Ta vi ế ạ t l i (11) d (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (12) (cid:0) y y y )c( y )x(yh i )x(y i )x( i i d ng: h 2
i v.v.
h !m ả ể ế ể (cid:0) i, y(cid:0)
)i( 4
)i( 3
)i( 2
i
ằ h !m ơ ể ặ ng pháp RungeKutta b ng cách đ t: (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) Ta đã kéo dài khai tri n Taylor đ k t qu chính xác h n. Đ tính y ể ta có th dùng ph y y (13) ươ )i( 1 kr 1 kr 4 kr 2 kr 3
1i trong đó: )i( k 1
)i( 2
i
i
)i( )k 1
i
i
)i( 1
)i( )k 2
)i( 3 .......
(cid:0) (cid:0) )y,x(hf i i (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) k x(hf y,ah (cid:0) (14) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) k x(hf y,bh k (cid:0) (cid:0)
(cid:0) ế , (cid:0) ,...; r 1, r 2,.. sao cho v ph i c a (13)
ầ ớ ế ả ủ ố ớ ệ ố ộ ể ấ ấ
ậ
)i( 1
)i( 2
)i( )k 1
(cid:0) (cid:0) (cid:0) k , (cid:0) ị và ta c n xác đ nh các h s a, b,..; ả ủ 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 RungeKutta b c hai ta có: )y,x(hf i i (cid:0) (15) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) k x(hf (cid:0)
1i
i kr 2
i kr 1
)i( 2
i
(cid:0) (cid:0) (cid:0) (cid:0) y (16) y,ah )i( 1
y và Ta có:
(cid:0))x(y,xf
y
172
(cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) y(cid:0) (x) = f[x,y(x)] )x(y )x(y,xf x
................
2
ế Do đó v ph i c a (12) là:
i
i
(cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (17) )y,x(hf i i )y,x(f x i i )x(y)y,x(f y i ả ủ h 2 ặ ứ (cid:0) (cid:0) (cid:0)
)i( 1
y
x
i
(cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) k k ] M t khác theo (15) và theo công th c Taylor ta có: yh )y,x(fah i i )y,x(fk i
x2
yi
i
i
i
i
(cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0) (cid:0)
3) ta tìm đ
)]y,x(fyr 2 ấ (18) ượ c các
i cho (17) và (18) khác nhau m t vô cùng bé c p O(h 2:
)i( 1 )i( 2 ế r(h 1 ờ Bây gi h s ch a bi
ệ ố ư ứ ứ ằ )y,x(hf i i )y,x(f[h i i ả ủ Do đó v ph i c a (16) là: 2 )y,x(far[h)y,x(f)r 2 i ộ ố ạ t khi cân b ng các s h ng ch a h và ch a h
(cid:0) ế r 1 + r 2 = 1 a.r 1 = 1/ 2 .r2 = 1 (cid:0) ọ
ế ượ ượ ấ ứ c ch n b t kì. ế c công th c Euler. N u = a, r 1 = (2a 1)/ 2a, r 2 = 1/ 2a v i a đ 1 = 0 và r 2 = 1. Lúc này ta nh n đ
ớ ậ ượ ứ
ậ ượ ươ ả ế c công th c Euler c i ti n. ậ c công th c Runge Kutta b c 4. chúng ta nh n đ
ộ ứ ự ế ư ậ Nh v y: N u a = 1 / 2 thì r a=1 thì r 1 = 1 / 2 và r 2 = 1/2. Lúc này ta nh n đ ậ ự c dùng trong tính toán th c t M t cách t Công th c này hay đ ứ :
ng t ượ k1 = h.f(xi, yi) k2 = h.f(xi+h/ 2, yi + k 1/ 2) k3 = h.f(xi+h/ 2, yi + k 2/ 2) k4 = h.f(xi+h, yi + k 3) yi+1 = yi + (k1 + 2k2 + 2k3 + k 4) / 6 ằ ươ ả ứ ậ i ph ng trình vi phân b ng công th c Runge Kutta b c 4 ng trình gi
ươ Ch ư nh sau:
ươ Ch ng trình 72
/ / P h uong phap Runge_Kutta;
#include
float f(float x,float y) {
float a=x+y; return(a);
173
}
void main() {
float a,b,k1,k2,k3,k4; int i,n; float x0,y0,h,e; float x[k],y[k];
clrscr(); printf("Phuong phap Runge Kutta \n"); printf("Cho can duoi a = "); scanf("%f",&a); printf("Cho can tren b = "); scanf("%f",&b); printf("Cho so kien y0 = "); scanf("%f",&y[0]); printf("Cho buoc tinh h = "); scanf("%f",&h); n=(int)((ba)/h); printf(" x y\n"); for (i=0;i<=n+1;i++) {
x[i]=a+i*h; 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; printf("%12.1f%16.4f\n",x[i],y[i]);
} getch();
}
o = 1 là :
ế ả ớ K t qu tính toán v i f = x + y, h = 0.1, a = 0, b =1, y
174
x 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 y 1.0000 1.1103 1.2427 1.3996 1.5834 1.7971 2.0440 2.3273 2.6508
175
0.9 1.0 3.0190 3.4362
