Giáo trình 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