際際滷

際際滷Share a Scribd company logo

Chuong13

H畛 L畛i
H畛 L畛i
1 of 7
Download to read offline
211 Ch測ng 13 : Gi其i ph測ng trnh vi ph息n 則1.B袖i to存n Cauchy M辿t ph測ng trnh vi ph息n cp 1 c達 th vit d鱈i d孫ng gi其i 速樽c y = f(x,y) m袖 ta c達 th tm 速樽c h袖m y t探 速孫o h袖m c単a n達.T奪n t孫i v束 s竪 nghim tho其 m揃n ph測ng trnh tr捉n.M巽i nghim 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 nhn 速樽c m辿t nghim ri捉ng c単a ph測ng trnh.B袖i to存n Cauchy ( hay b袖i to存n c達 速iu kin 速u) t達m l孫i nh sau : cho x sao cho b x a,tm y(x) tho其 m揃n 速iu kin : ｩ ｨ ｧ 留= = )a(y )y,x(f)x(y (1) Ng棚i ta ch淡ng minh r損ng b袖i to存n n袖y c達 m辿t nghim duy nht nu f tho其 m揃n 速iu kin Lipschitz : 1 2 1 2 f x y f x y L y y( , ) ( , ) v鱈i L l袖 m辿t h損ng s竪 d測ng. Ng棚i ta c嘆ng ch淡ng minh r損ng nu fy ( 速孫o h袖m c単a f theo y ) l袖 li捉n t担c v袖 b chn th f tho其 m揃n 速iu kin Lipschitz. M辿t c存ch t脱ng qu存t h測n,ng棚i ta 速nh ngha h ph測ng trnh bc 1 : 1 1 1 2 , ( ), , ,...,y f x y y yn = 2 2 1 2 , ( ), , ,...,y f x y y yn = ................ n n n y f x y y y , ( ), , ,...,= 1 2 Ta ph其i tm nghim y1,y2,...,yn sao cho : = = ｧ ｨ ｩ Y x f x Y Y a ( ) ( , ) ( ) 留 v鱈i : = Y y y yn 1 2 , , , . . F f f fn = 1 2 . . Y y y yn = 1 2 . . Nu ph測ng trnh vi ph息n c達 bc cao h測n (n),nghim s ph担 thu辿c v袖o n h損ng s竪 tu炭 箪.則 nhn 速樽c m辿t nghim ri捉ng,ta ph其i cho n 速iu kin 速u.B袖i to存n s c達 gi存 tr 速u nu v鱈i gi存 tr xo 速揃 cho ta cho y(xo),y(xo),y(xo),.... M辿t ph測ng trnh vi ph息n bc n c達 th 速a v th袖nh m辿t h ph測ng trnh vi ph息n cp 1.V d担 nu ta c達 ph測ng trnh vi ph息n cp 2 : 霞 = = = ｧ ｨ ｪ ｩ ｪ y f x y y y a y a ( , , ) ( ) ( ),留 硫 Khi 速t u = y v袖 v = y ta nhn 速樽c h ph測ng trnh vi ph息n cp 1 : = = ｧ ｨ ｩ u v v g x u v( , , ) t鱈i 速iu kin 速u : u(a) = 留 v袖 v(a) = 硫 C存c ph測ng ph存p gi其i ph測ng trnh vi ph息n 速樽c trnh b袖y trong ch測ng n袖y l袖
212 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" tch ph息n h = ( b - a) / n. 則2.Ph測ng ph存p Euler v袖 Euler c其i tin Gi其 s旦 ta c達 ph測ng trnh vi ph息n : = = ｧ ｨ ｩ y x f x y y a ( ) ( , ) ( ) 留 (1) v袖 cn tm nghim c単a n達.Ta chia 速o孫n [xo,x ] th袖nh n phn b谷i c存c 速im chia : xo < x1 < x2 <...< xn = x Theo c束ng th淡c khai trin Taylor m辿t h袖m l息n cn xi ta c達 : i i i i i i i i i i i y x y x x x y x x x y x x x y x + + + + = + + + + 霞 霞霞 1 1 1 2 1 3 2 6 ( ) ( ( ) ( ) ( ) ( ) ( ) ( ) ) Nu (xi+1 - xi) kh存 b th ta c達 th b叩 qua c存c s竪 h孫ng (xi+1 - xi)2 v袖 c存c s竪 h孫ng bc cao y(xi+1) = y(xi) + (xi+1- xi) y(xi) Tr棚ng h樽p c存c m竪c c存ch 速u : (xi-1 - xi) = h = (x - xo)/ n th ta nhn 速樽c c束ng th淡c Euler 速測n gi其n : yi+1 = yi + hf(xi,yi) (2) V mt hnh h辰c ta thy (1) cho kt qu其 c袖ng chnh x存c nu b鱈c h c袖ng nh叩.則 t即ng 速辿 chnh x存c ta c達 th d誰ng c束ng th淡c Euler c其i tin.Tr鱈c ht ta nh他c l孫i 速nh l Lagrange: Gi其 s旦 f(x) l袖 h袖m li捉n t担c trong[a,b] v袖 kh其 vi trong (a,b)th c達 t nht m辿t 速im c(a,b) 速 cho : ab )a(f)b(f )c(f = Theo 速nh l Lagrange ta c達 : ))c(y,c(hf)x(y)x(y iii1i +=+ Nh vy v鱈i c(xi,xi+1) ta c達 th thay : [ ])y,x(f)y,x(f 2 1 ))c(y,c(f 1i1iiiii +++= T探 速達 ta c達 c束ng th淡c Euler c其i tin : i i i i i i y y h f x y f x y+ + + = + + 1 1 12 [ ( ) ( )], , (3) Trong c束ng th淡c n袖y gi存 tr yi+1 cha bit.Do 速達 khi 速揃 bit yi ta ph其i tm yi+1 b損ng c存ch gi其i ph測ng trnh 速孫i s竪 tuyn tnh (3).Ta th棚ng gi其i (3) b損ng c存ch lp nh sau:tr鱈c ht ch辰n xp x 速u ti捉n c単a php lp )0( 1iy + chnh l袖 gi存 tr yi+1 tnh 速樽c theo ph測ng ph存p Euler sau 速達 d誰ng (3) 速 tnh c存c )s( 1iy + ,c担 th l袖 : [ ])y,x(f)y,x(f 2 h yy )y,x(hfyy )1s( 1i1iiii )s( 1i iii )0( 1i +++ + ++= += y b a yi yi+1 h xi xi+1 x
213 Qu存 trnh tnh kt th坦c khi )s( iy 速単 gn )1s( iy Ch測ng trnh gi其i ph測ng trnh vi ph息n theo ph測ng ph存p Euler nh sau : Ch測ng trnh 13-1 //pp_Euler; #include <conio.h> #include <stdio.h> #include <math.h> 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]; 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 quan"); printf("n"); printf("Phuong phap Eulern"); h=(b-a)/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"); }
214 printf("n"); getch(); printf("Phuong phap Euler cai tienn"); 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(); } V鱈i ph測ng trnh cho trong function v袖 速iu kin 速u xo = 0,yo= 0, nghim trong 速o孫n [0,1] v鱈i 10 速im chia l袖 : x y(Euler) y(Euler c其i tin) 0.0 0.00 0.00 0.1 0.00 0.01 0.2 0.01 0.02 0.3 0.03 0.05 0.4 0.06 0.09 0.5 0.11 0.15 0.6 0.17 0.22 0.7 0.25 0.31 0.8 0.34 0.42 0.9 0.46 0.56 1.0 0.59 0.71 則3.Ph測ng ph存p Runge-Kutta Xt b袖i to存n Cauchy (1).Gi其 s旦 ta 速揃 tm 速樽c gi存 tr gn 速坦ng yi c単a y(xi) v袖 mu竪n tnh yi+1 c単a y(xi+1).Tr鱈c ht ta vit c束ng th淡c Taylor : i i i i m i m y x y x hy x h y x h m y x h m ym m + + = + + + + + + 霞 +1 2 1 2 1 1 ( ) ( ) ( ) ( ) ! ( ) ( )! (c)... ( ) ( ) ( ) (11) v鱈i c (xi,xi+1) v袖 : i i i y x f x y x =( ) [ ( )], ( ) ( ) [ ( ( )], k i k ky x d dx f x y x ix x = = 1 1 Ta vit l孫i (11) d鱈i d孫ng : i i i i m i m m m y y hy h y h m y h m y+ + + = + + + + +1 2 1 1 2 1 , ,, ( ) ( ) ( ) ... ! ( )! (c) (12) Ta 速揃 ko d袖i khai trin Taylor 速 kt qu其 chnh x存c h測n.則 tnh yi,yi v.v.ta c達 th d誰ng ph測ng ph存p Runge-Kutta b損ng c存ch 速t :
215 i i i i i s s i y y r k r k r k r k+ = + + + + 1 1 1 2 2 3 3 ( ) ( ) ( ) ( ) ... (13) trong 速達 : 1 2 1 3 1 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) . . . . . . . . . . . . . . . , , , i i i i i i i i i i i i k hf x y k hf x ah y k k hf x bh y k k = = + + = + + + ｧ ｨ ｪ ｪ ｪ ｪ ｩ ｪ ｪ ｪ ｪ 留 硫 粒 (14) v袖 ta cn 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 cp cao nht c達 th c達 速竪i v鱈i h. Khi d誰ng c束ng th淡c Runge-Kutta bc hai ta c達 : 1 2 1 ( ) ( ) ( ) ( ) ( ) , , i i i i i i i k hf x y k hf x ah y k = = + + ｧ ｨ ｪ ｩ ｪ 留 (15) v袖 i i i i y y r k r k+ = + 1 1 1 2 2 ( ) ( ) (16) Ta c達 : y(x) = f[x,y(x)] 霞 = +y x f x y x f x y x y xx y( ) [ , ( )] [ , ( )] ( ), , ................ Do 速達 v ph其i c単a (12) l袖 : i i x i i y i i hf x y h f x y f x y y x( ) [ ( ) ( ) ( )], , , ... , , + + + 2 2 (17) Mt kh存c theo (15) v袖 theo c束ng th淡c Taylor ta c達 : 1 ( ) , ( ), i i i ik hf x y hy= = 2 1 ( ) , ( ) , [ ( ) ( ) ( ) ], , , ..... i i i x i i i y i ik h f x y ahf x y k f x y= + + +留 Do 速達 v ph其i c単a (16) l袖 : 1 2 2 2 2h r r x y h f x y r y f x yi i x i i i x i i ( )f( ) [ar ( ) ( )], , , .... , , , + + + +留 (18) B息y gi棚 cho (17) v袖 (18) kh存c nhau m辿t v束 c誰ng b cp O(h3 ) ta tm 速樽c c存c h s竪 cha bit 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 vy : 留 = a,r1 = (2a - 1)/ 2a,r2 = 1/ 2a v鱈i a 速樽c ch辰n bt k. Nu a = 1 / 2 th r1 = 0 v袖 r2 = 1.L坦c n袖y ta nhn 速樽c c束ng th淡c Euler.Nu a = 1 th r1 = 1 / 2 v袖 r2 = 1/2.L坦c n袖y ta nhn 速樽c c束ng th淡c Euler c其i tin. M辿t c存ch t測ng t湛 ch坦ng ta nhn 速樽c c束ng th淡c Runge - Kutta bc 4.C束ng th淡c n袖y hay 速樽c d誰ng trong tnh 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 Ch測ng trnh gi其i ph測ng trnh vi ph息n b損ng c束ng th淡c Runge - Kutta bc 4 nh sau : Ch測ng trnh 11-2 //Phuong phap Runge_Kutta;
216 #include <conio.h> #include <stdio.h> #include <math.h> #define k 10 float f(float x,float y) { float a=x+y; return(a); } 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 - Kuttan"); 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)((b-a)/h); printf(" x yn"); 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.4fn",x[i],y[i]); } getch(); } Kt qu其 tnh to存n v鱈i f = x + y,h = 0.1,a = 0,b =1,yo = 1 l袖 : x y 0.0 1.0000 0.1 1.1103 0.2 1.2427 0.3 1.3996 0.4 1.5834
217 0.5 1.7971 0.6 2.0440 0.7 2.3273 0.8 2.6508 0.9 3.0190 1.0 3.4362

More Related Content

Chuong13

  • 1. 211 Ch測ng 13 : Gi其i ph測ng trnh vi ph息n 則1.B袖i to存n Cauchy M辿t ph測ng trnh vi ph息n cp 1 c達 th vit d鱈i d孫ng gi其i 速樽c y = f(x,y) m袖 ta c達 th tm 速樽c h袖m y t探 速孫o h袖m c単a n達.T奪n t孫i v束 s竪 nghim tho其 m揃n ph測ng trnh tr捉n.M巽i nghim 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 nhn 速樽c m辿t nghim ri捉ng c単a ph測ng trnh.B袖i to存n Cauchy ( hay b袖i to存n c達 速iu kin 速u) t達m l孫i nh sau : cho x sao cho b x a,tm y(x) tho其 m揃n 速iu kin : ｩ ｨ ｧ 留= = )a(y )y,x(f)x(y (1) Ng棚i ta ch淡ng minh r損ng b袖i to存n n袖y c達 m辿t nghim duy nht nu f tho其 m揃n 速iu kin Lipschitz : 1 2 1 2 f x y f x y L y y( , ) ( , ) v鱈i L l袖 m辿t h損ng s竪 d測ng. Ng棚i ta c嘆ng ch淡ng minh r損ng nu fy ( 速孫o h袖m c単a f theo y ) l袖 li捉n t担c v袖 b chn th f tho其 m揃n 速iu kin Lipschitz. M辿t c存ch t脱ng qu存t h測n,ng棚i ta 速nh ngha h ph測ng trnh bc 1 : 1 1 1 2 , ( ), , ,...,y f x y y yn = 2 2 1 2 , ( ), , ,...,y f x y y yn = ................ n n n y f x y y y , ( ), , ,...,= 1 2 Ta ph其i tm nghim y1,y2,...,yn sao cho : = = ｧ ｨ ｩ Y x f x Y Y a ( ) ( , ) ( ) 留 v鱈i : = Y y y yn 1 2 , , , . . F f f fn = 1 2 . . Y y y yn = 1 2 . . Nu ph測ng trnh vi ph息n c達 bc cao h測n (n),nghim s ph担 thu辿c v袖o n h損ng s竪 tu炭 箪.則 nhn 速樽c m辿t nghim ri捉ng,ta ph其i cho n 速iu kin 速u.B袖i to存n s c達 gi存 tr 速u nu v鱈i gi存 tr xo 速揃 cho ta cho y(xo),y(xo),y(xo),.... M辿t ph測ng trnh vi ph息n bc n c達 th 速a v th袖nh m辿t h ph測ng trnh vi ph息n cp 1.V d担 nu ta c達 ph測ng trnh vi ph息n cp 2 : 霞 = = = ｧ ｨ ｪ ｩ ｪ y f x y y y a y a ( , , ) ( ) ( ),留 硫 Khi 速t u = y v袖 v = y ta nhn 速樽c h ph測ng trnh vi ph息n cp 1 : = = ｧ ｨ ｩ u v v g x u v( , , ) t鱈i 速iu kin 速u : u(a) = 留 v袖 v(a) = 硫 C存c ph測ng ph存p gi其i ph測ng trnh vi ph息n 速樽c trnh b袖y trong ch測ng n袖y l袖
  • 2. 212 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" tch ph息n h = ( b - a) / n. 則2.Ph測ng ph存p Euler v袖 Euler c其i tin Gi其 s旦 ta c達 ph測ng trnh vi ph息n : = = ｧ ｨ ｩ y x f x y y a ( ) ( , ) ( ) 留 (1) v袖 cn tm nghim c単a n達.Ta chia 速o孫n [xo,x ] th袖nh n phn b谷i c存c 速im chia : xo < x1 < x2 <...< xn = x Theo c束ng th淡c khai trin Taylor m辿t h袖m l息n cn xi ta c達 : i i i i i i i i i i i y x y x x x y x x x y x x x y x + + + + = + + + + 霞 霞霞 1 1 1 2 1 3 2 6 ( ) ( ( ) ( ) ( ) ( ) ( ) ( ) ) Nu (xi+1 - xi) kh存 b th ta c達 th b叩 qua c存c s竪 h孫ng (xi+1 - xi)2 v袖 c存c s竪 h孫ng bc cao y(xi+1) = y(xi) + (xi+1- xi) y(xi) Tr棚ng h樽p c存c m竪c c存ch 速u : (xi-1 - xi) = h = (x - xo)/ n th ta nhn 速樽c c束ng th淡c Euler 速測n gi其n : yi+1 = yi + hf(xi,yi) (2) V mt hnh h辰c ta thy (1) cho kt qu其 c袖ng chnh x存c nu b鱈c h c袖ng nh叩.則 t即ng 速辿 chnh x存c ta c達 th d誰ng c束ng th淡c Euler c其i tin.Tr鱈c ht ta nh他c l孫i 速nh l Lagrange: Gi其 s旦 f(x) l袖 h袖m li捉n t担c trong[a,b] v袖 kh其 vi trong (a,b)th c達 t nht m辿t 速im c(a,b) 速 cho : ab )a(f)b(f )c(f = Theo 速nh l Lagrange ta c達 : ))c(y,c(hf)x(y)x(y iii1i +=+ Nh vy v鱈i c(xi,xi+1) ta c達 th thay : [ ])y,x(f)y,x(f 2 1 ))c(y,c(f 1i1iiiii +++= T探 速達 ta c達 c束ng th淡c Euler c其i tin : i i i i i i y y h f x y f x y+ + + = + + 1 1 12 [ ( ) ( )], , (3) Trong c束ng th淡c n袖y gi存 tr yi+1 cha bit.Do 速達 khi 速揃 bit yi ta ph其i tm yi+1 b損ng c存ch gi其i ph測ng trnh 速孫i s竪 tuyn tnh (3).Ta th棚ng gi其i (3) b損ng c存ch lp nh sau:tr鱈c ht ch辰n xp x 速u ti捉n c単a php lp )0( 1iy + chnh l袖 gi存 tr yi+1 tnh 速樽c theo ph測ng ph存p Euler sau 速達 d誰ng (3) 速 tnh c存c )s( 1iy + ,c担 th l袖 : [ ])y,x(f)y,x(f 2 h yy )y,x(hfyy )1s( 1i1iiii )s( 1i iii )0( 1i +++ + ++= += y b a yi yi+1 h xi xi+1 x
  • 3. 213 Qu存 trnh tnh kt th坦c khi )s( iy 速単 gn )1s( iy Ch測ng trnh gi其i ph測ng trnh vi ph息n theo ph測ng ph存p Euler nh sau : Ch測ng trnh 13-1 //pp_Euler; #include <conio.h> #include <stdio.h> #include <math.h> 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]; 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 quan"); printf("n"); printf("Phuong phap Eulern"); h=(b-a)/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"); }
  • 4. 214 printf("n"); getch(); printf("Phuong phap Euler cai tienn"); 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(); } V鱈i ph測ng trnh cho trong function v袖 速iu kin 速u xo = 0,yo= 0, nghim trong 速o孫n [0,1] v鱈i 10 速im chia l袖 : x y(Euler) y(Euler c其i tin) 0.0 0.00 0.00 0.1 0.00 0.01 0.2 0.01 0.02 0.3 0.03 0.05 0.4 0.06 0.09 0.5 0.11 0.15 0.6 0.17 0.22 0.7 0.25 0.31 0.8 0.34 0.42 0.9 0.46 0.56 1.0 0.59 0.71 則3.Ph測ng ph存p Runge-Kutta Xt b袖i to存n Cauchy (1).Gi其 s旦 ta 速揃 tm 速樽c gi存 tr gn 速坦ng yi c単a y(xi) v袖 mu竪n tnh yi+1 c単a y(xi+1).Tr鱈c ht ta vit c束ng th淡c Taylor : i i i i m i m y x y x hy x h y x h m y x h m ym m + + = + + + + + + 霞 +1 2 1 2 1 1 ( ) ( ) ( ) ( ) ! ( ) ( )! (c)... ( ) ( ) ( ) (11) v鱈i c (xi,xi+1) v袖 : i i i y x f x y x =( ) [ ( )], ( ) ( ) [ ( ( )], k i k ky x d dx f x y x ix x = = 1 1 Ta vit l孫i (11) d鱈i d孫ng : i i i i m i m m m y y hy h y h m y h m y+ + + = + + + + +1 2 1 1 2 1 , ,, ( ) ( ) ( ) ... ! ( )! (c) (12) Ta 速揃 ko d袖i khai trin Taylor 速 kt qu其 chnh x存c h測n.則 tnh yi,yi v.v.ta c達 th d誰ng ph測ng ph存p Runge-Kutta b損ng c存ch 速t :
  • 5. 215 i i i i i s s i y y r k r k r k r k+ = + + + + 1 1 1 2 2 3 3 ( ) ( ) ( ) ( ) ... (13) trong 速達 : 1 2 1 3 1 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) . . . . . . . . . . . . . . . , , , i i i i i i i i i i i i k hf x y k hf x ah y k k hf x bh y k k = = + + = + + + ｧ ｨ ｪ ｪ ｪ ｪ ｩ ｪ ｪ ｪ ｪ 留 硫 粒 (14) v袖 ta cn 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 cp cao nht c達 th c達 速竪i v鱈i h. Khi d誰ng c束ng th淡c Runge-Kutta bc hai ta c達 : 1 2 1 ( ) ( ) ( ) ( ) ( ) , , i i i i i i i k hf x y k hf x ah y k = = + + ｧ ｨ ｪ ｩ ｪ 留 (15) v袖 i i i i y y r k r k+ = + 1 1 1 2 2 ( ) ( ) (16) Ta c達 : y(x) = f[x,y(x)] 霞 = +y x f x y x f x y x y xx y( ) [ , ( )] [ , ( )] ( ), , ................ Do 速達 v ph其i c単a (12) l袖 : i i x i i y i i hf x y h f x y f x y y x( ) [ ( ) ( ) ( )], , , ... , , + + + 2 2 (17) Mt kh存c theo (15) v袖 theo c束ng th淡c Taylor ta c達 : 1 ( ) , ( ), i i i ik hf x y hy= = 2 1 ( ) , ( ) , [ ( ) ( ) ( ) ], , , ..... i i i x i i i y i ik h f x y ahf x y k f x y= + + +留 Do 速達 v ph其i c単a (16) l袖 : 1 2 2 2 2h r r x y h f x y r y f x yi i x i i i x i i ( )f( ) [ar ( ) ( )], , , .... , , , + + + +留 (18) B息y gi棚 cho (17) v袖 (18) kh存c nhau m辿t v束 c誰ng b cp O(h3 ) ta tm 速樽c c存c h s竪 cha bit 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 vy : 留 = a,r1 = (2a - 1)/ 2a,r2 = 1/ 2a v鱈i a 速樽c ch辰n bt k. Nu a = 1 / 2 th r1 = 0 v袖 r2 = 1.L坦c n袖y ta nhn 速樽c c束ng th淡c Euler.Nu a = 1 th r1 = 1 / 2 v袖 r2 = 1/2.L坦c n袖y ta nhn 速樽c c束ng th淡c Euler c其i tin. M辿t c存ch t測ng t湛 ch坦ng ta nhn 速樽c c束ng th淡c Runge - Kutta bc 4.C束ng th淡c n袖y hay 速樽c d誰ng trong tnh 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 Ch測ng trnh gi其i ph測ng trnh vi ph息n b損ng c束ng th淡c Runge - Kutta bc 4 nh sau : Ch測ng trnh 11-2 //Phuong phap Runge_Kutta;
  • 6. 216 #include <conio.h> #include <stdio.h> #include <math.h> #define k 10 float f(float x,float y) { float a=x+y; return(a); } 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 - Kuttan"); 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)((b-a)/h); printf(" x yn"); 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.4fn",x[i],y[i]); } getch(); } Kt qu其 tnh to存n v鱈i f = x + y,h = 0.1,a = 0,b =1,yo = 1 l袖 : x y 0.0 1.0000 0.1 1.1103 0.2 1.2427 0.3 1.3996 0.4 1.5834
  • 7. 217 0.5 1.7971 0.6 2.0440 0.7 2.3273 0.8 2.6508 0.9 3.0190 1.0 3.4362