Giáo trình Turbo C nâng cao và C++ - Chương 13: Giải phương trình vi phân

211 Ch−ơng 13 : Giải ph−ơng trình vi phân Đ1.Bài toá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 y′ = f(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 : ⎩⎨ ⎧ α== ′ )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 nghiệm duy nhất nếu f thoả mãn điều kiện Lipschitz : 1 2 1 2f 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 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 : 1 1 1 2 , ( ), , , ...,y f x y y yn= 2 2 1 2 , ( ), , , ...,y f x y y yn= ................ n n ny f x y y y , ( ), , , ...,= 1 2 Ta phải tìm nghiệm y1,y2,...,yn sao cho : ′ == ⎧⎨⎩ Y x f x Y Y a ( ) ( , ) ( ) α với : ′ = ⎛ ⎝ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜ ⎞ ⎠ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟ Y y y y n 1 2 , , , . . F f f f n = ⎛ ⎝ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜ ⎞ ⎠ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟ 1 2 . . Y y y y n = ⎛ ⎝ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜ ⎞ ⎠ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟ 1 2 . . 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 ( , , ) ( ) ( ),α β Khi đặt u = y và v = y′ ta nhận đ−ợc hệ ph−ơng trình vi phân cấp 1 : ′ = ′ = ⎧⎨⎩ u v v g x u v( , , ) tớ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à 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" tích phân h = ( b - a) / n. Đ2.Ph−ơng pháp Euler và Euler cải tiến Giả sử ta có ph−ơng trình vi phân : ′ = = ⎧⎨⎩ y x f x y y a ( ) ( , ) ( ) α (1) và cần tìm nghiệm của nó.Ta chia đoạn [xo,x ] thành n phần bởi các điểm chia : xo < x1 < x2 <...< xn = x Theo công thức khai triển Taylor một hàm lân cận xi ta có : i i i i i i i i i i iy x y x x x y x x x y x x x y x + + + += + − + − + − + ⋅ ⋅ ⋅′ ′′ ′′′ 1 1 1 2 1 3 2 6 ( ) ( ( ) ( ) ( ) ( ) ( ) ( ) ) Nếu (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 bậc 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 nhận đ−ợc công thức Euler đơn giản : yi+1 = yi + hf(xi,yi) (2) Về mặt hình học ta thấy (1) cho kết quả càng chính xác nếu b−ớc h càng nhỏ.Để tăng độ chính xác ta có thể dùng công thức Euler cải tiến.Tr−ớc hết 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 nhất một điểm 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− vậy 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 tiến : 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 ch−a biết.Do đó khi đã biết yi ta phải tìm yi+1 bằng cách giải ph−ơng trình đại số tuyến tính (3).Ta th−ờng giải (3) bằng cách lặp nh− sau:tr−ớc hết chọn xấp xỉ đầu tiên của phép lặp )0( 1iy + chính là giá trị yi+1 tính đ−ợc theo ph−ơng pháp Euler sau đó dùng (3) để tính 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á trình tính kết thúc khi )s(iy đủ gần )1s( iy − 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 13-1 //pp_Euler; #include #include #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]; 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=(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 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(); } Với ph−ơng trình cho trong function và điều kiện đầu xo = 0,yo= 0, nghiệm trong đoạn [0,1] với 10 điểm chia là : x y(Euler) y(Euler cải tiến) 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 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 : i i i i m i m y x y x hy x h y x h m y x h m y m 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 viết lại (11) d−ới dạng : i i i i m i m m my y hy h y h m y h m y+ + +− = + + + + +1 2 1 1 2 1 , ,, ( ) ( ) ( ) ... ! ( )! (c) (12) 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 : 215 i i i i i s s iy 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 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ó : 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 iy 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) Mặt 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é 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 Ch−ơng trình giải ph−ơng trình vi phân bằng công thức Runge - Kutta bậc 4 nh− sau : Ch−ơng trình 11-2 //Phuong phap Runge_Kutta; 216 #include #include #include #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 - 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)((b-a)/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(); } Kết quả tính 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