1
BÁO CÁO TỐT NGHIỆP

Về một phương pháp không cổ điển giải số
phương trình vi phân bậc nhất và bậc hai 2

LỜI NÓI ĐẦU Phương trình vi phân là mô hình mô tả khá tốt các quá trình chuyển động
trong tự nhiên và kĩ thuật. Để nghiên cứu phương trình vi phân, người ta thường
tiếp cận theo hai hướng: nghiên cứu định tính và giải số.

tính toán trên máy các ví dụ của M. V. Bulatov và G. V. Berghe.
Luận văn được hoàn thành dưới sự hướng dẫn khoa học của PGS-TS Tạ Duy
Phượng (Viện Toán học). Xin được tỏ lòng cám ơn chân thành nhất tới Thầy.
Tác giả xin tỏ lòng cám ơn Ban Chủ nhiệm , các Thày Cô và các cán bộ khoa
Toán- Cơ – Tin học, trường Đại học Khoa học Tự nhiên, Đại học Quốc gia Hà Nội
đã nhiệt tình giảng dạy và giúp đỡ tác giả trong suốt quá trình học Cao học.
Tác giả xin chân thành cám ơn Ban Lãnh đạo và các cán bộ, giáo viên Học
viện Quân y đã tạo mọi điều kiện để tác giả hoàn thành tốt khóa học Cao học.
Và cuối cùng, xin được cám ơn Gia đình, bạn bè đã thông cảm, sẻ chia, hy
sinh và tạo mọi điều kiện cho tác giả trong suốt thời gian học Cao học và viết luận
văn.

Hà Nội, ngày 30 tháng 11 năm 2009
Tác giả
Vũ Thị Thanh Bình 4

CHƯƠNG 1
Kiến thức chuẩn bị

Trong Chương 1 chúng tôi nhắc lại những khái niệm cơ bản nhất của giải số
phương trình vi phân nhằm thuận tiện cho trình bày ở các mục sau.
1.1. Bài toán Cauchy giải hệ phương trình vi phân
Xét bài toán Cauchy tìm nghiệm của hệ phương trình


Ở đây ta hiểu nghiệm theo nghĩa cổ điển và địa phương, tức là nghiệm của (1.1)-
(1.2) là một hàm khả vi
( )
x t
trên


0,

,
1


sao cho
( ) ( ( ), )
x t f x t t


trên


0,

và
0
(0)
x x

.
Cùng với bài toán (1.1), ta cũng xét trường hợp hàm



,
f x t
,
( )
g t
là
đủ trơn (có đạo hàm đến cấp cần thiết trong tính toán). Khi ấy theo định lí Picard-
Lindelöf, hệ (1.1)-(1.2) có duy nhất nghiệm
( )
x t
trên toàn đoạn


0,1
(nghiệm có
thể kéo dài được trên toàn bộ khoảng xác định, hay tồn tại nghiệm toàn cục, xem
[8], trang 467). Lưu ý này là quan trọng trong giải số hệ phương trình (1.1)-(1.2). 5

1.2. Giải số bài toán Cauchy
Để chứng minh định lý về sự tồn tại và duy nhất nghiệm của hệ phương trình vi
phân (1.1)-(1.2), ta có thể xây dựng dãy nghiệm xấp xỉ hội tụ tới nghiệm của bài
toán (1.1)-(1.2) trên khoảng tồn tại nghiệm. Có hai phương pháp xây dựng dãy
nghiệm xấp xỉ: phương pháp giải tích và phương pháp số kết quả được cho dưới
dạng bảng, như phương pháp Euler, phương pháp Runge-Kutta, phương pháp đa
bước,

( )
f t
được
xấp xỉ bởi tích phân của hàm đa thức (tính được chính xác).
Giả sử ta có
s
điểm nội suy khác nhau
1 2
, , ,
s
c c c
trong khoảng


,
a b
. Đa thức nội
suy Lagrange bậc nhỏ hơn
s
có dạng (xem [1]):
1
( ) ( ) ( )
s
j j
j
t f c L t





f t dt f c





.
Các trọng số
j

được tính theo công thức
( ) .
b
j j
a
L t dt




Nếu
1
s

thì đa thức nội suy
1
( ) ( )
t f c

 và ta có:






Nếu
0( )
b a h
 
thì sai số trong quy tắc cầu phương của độ chính xác
p
là
1
0( ).
p
h


Ta xét một số trường hợp đặc biệt.
 Nếu chọn
1
s

và
1
c a

thì ta có công thức xấp xỉ tích phân bởi diện tích hình
chữ nhật ABCD (Hình 1.1):


(
h

có thể dương hoặc âm, khi
h

dương thì nghiệm được xây dựng về bên phải của điểm
0
t

và ngược lại, khi
h

âm

7

thì nghiệm được xây dựng về bên trái của
0
t
). Dưới đây ta coi
0
h

, trường hợp
0
h

có thể được xét tương tự.
Từ công thức (1.7) ta có

a
f t dt b a f b
 


Từ đây ta có:
( ) ( ) . ( ( ), )
x t h x t h f x t h t h
    

Suy ra công thức Euler lùi:
1 1 1
. ( , )
n n n n
x x h f x t
  
 
. Hình 1.2
Hai phương pháp Euler tiến và Euler lùi là những phương pháp Runge-Kutta bậc
nhất (có độ xấp xỉ bậc nhất).
 Nếu chọn
1
s

và
1
2
a b
c



b

B

a

A

E

b

B

a

A

F

O

f

x

2
ba


   
.
Đây chính là phương pháp trung điểm (midpoind method).
 Nếu chọn
2
s

và
1 2
,
c a c b
 
thì
1
( )
( )
t b
L t
a b



và
2
( ) .
( )
t a
L t
b a


a a
b
t a t a b a
L t dt dt
b a b a

  
   
 
 

Chứng tỏ
( ) [ ( ) ( )]
2
b
a
b a
f t dt f a f b

 

.
Như vậy nếu xấp xỉ tích phân
( ( ), )
t h
t
f x s s ds


bởi công thức trên (bởi diện tích hình

 Nếu chọn
3
s

và
1 2 3
, ,
2
a b
c a c c b

  
thì, đặt
h b a
 
, ta có:
E

D

b

B

a

A

O

2
2
( )( ) 4
( ) ( )( ),
( )( )
2 2
t a t b
L t t a t b
a b a b
h
a b
  
   
 
 3
2
( )( )
2
2
( ) ( )( ).
2
( )( )
2
a b
t a t

a a a
b
b
a
a
a b a b
L t dt t t b dt t b t b dt
h h
a b t b a b t b
t b t b dt
h h
b a h
h

 
       
   
     

 
  


và
2 2
2 2
3 2
2
4 4
( ) ( )( ) ( )( ( ))

a
h a b
f t dt f a f f b

  

.

Suy ra công thức xấp xỉ nghiệm của phương trình vi phân

10
( ) ( ) [ ( ( ), ) 4 ( ( ), ) ( ( ), )]
6 2 2
h h h
x t h x t f x t t f x t t f x t h t h
        

và công thức sai phân
1 1 1
[ ( , ) 4 ( ( ), ) ( , )]
6 2 2
n n n n n n n n
h h h
x x f x t f x t t f x t
  
     
.
Đây là công thức ẩn của phương pháp Runge-Kutta kinh điển cấp bốn (classical
fourth-order Runge-Kutta method).
1.2.2. Phương pháp Runge-Kutta

[ ( , ) ( , )]
2
n n n n n n
h
x x f x t f x t
  
  
Công thức này được gọi là phương pháp hình thang hiển (explicit trapezoidal
method).
 Bằng cách sử dụng xấp xỉ bậc nhất của
( )
2
n
h
x t

theo phương pháp Euler tiến:
1
2
ˆ
( , )
2
n n n
n
h
x x f x t

 

và thay vào công thức của phương pháp trung điểm ẩn

n n n n n n n n
h h h
x x f x t f x t t f x t
  
     
,
ta có công thức Runge-Kutta hiển bậc bốn kinh điển sau:
1 1 2 3 4
( 2 2 ), 0,1,2,
6
n n
h
x x k k k k n

     
trong đó:
1
1
2
2
3
4 3
( , );
( , );
2 2
( , );
2 2
( , ).
n n
n n

1
h
N

,
và kí hiệu
i
x
là giá trị xấp xỉ


i
x t
,
( )
i i
B B t
 ,


i i
g g t
 .
Phương pháp Runge-Kutta cho bài toán (1.1)-(1.2) có dạng (xem, [2], [4]-[7])
1
1
s
i i i i
i
x x h b X

ij
a
,
i
c
,
i
b
xác định bậc xấp xỉ của phương pháp, còn
s
được gọi là số

12
nấc. Nếu
0
ij
a

với mọi
j i

thì ta có phương pháp Runge-Kutta hiển. Khi ấy tính
toán khá đơn giản (
i
X
được tính theo công thức truy hồi). Nếu
0
ij
a


b

Hai phương pháp Runge-Kutta quan trọng thường hay được sử dụng là phương
pháp Runge-Kutta bậc hai và phương pháp Runge-Kutta bậc bốn.
1.2.2.3. Công thức lặp của phương pháp Runge-Kutta bậc hai
Giả thiết rằng ta đã biết giá trị của
x
tại
n
t
là
n
x
. Phương pháp Runge-Kutta hiển
hai nấc cấp hai sử dụng điểm
( , )
n n
x t

để xấp xỉ giá trị của
x
tại điểm tiếp theo bằng
công thức
1 1 1 2 2
( ),
n n
x x h b k b k

  
(1.10)

   
.

13
Đây là một hệ ba phương trình (phi tuyến) bốn ẩn. Ta có thể chọn một hệ số, thí
dụ,
2
0
b

tự do. Khi ấy các hệ số còn lại biểu diễn qua
2
0
b

bởi các công thức:
1 2
1
b b
 
,
2
2
1
2
c
b
 ,
21
2

1 1
( , ) ( ( , ) )
2 2
n n n n n n n n
x x hf x t hf x hf x t t h

    
.
Công thức này được gọi là phương pháp Runge-Kutta đơn giản (Simple Runge-
Kutta Method) hoặc phương pháp tiếp tuyến cải tiến (Impoved Tangent Method), vì
nó trùng với phương pháp Euler cải tiến.
Nếu chọn
2
1
b

thì
21
1
2
a

,
1
0
b

và
2
1

. (1.11)
Phương pháp một tựa tương ứng với nó là
1 1 1
0 0 0
( , )
k k k
j i j j i j j i j
j j j
x hf x t
  
     
  

  
. (1.12)
Đối với phương pháp này ta giả thiết rằng các giá trị xuất phát
1 2 1
, , ,
k
x x x

đã được
tính tương đối chính xác.
Nếu
0
0


và
0


   
, (1.13)
trong đó

là một hằng số (thực hoặc phức). Nghiệm của phương trình này là
0
)( xetx
t

 .
Ta thường viết
IR
i


, trong đó
R

và
I

tương ứng là các phần thực và phần
ảo của

.
Tương ứng với phương trình (1.13), xét phương trình sai phân tuyến tính bậc nhất
nn
xx


xxexex
nnh
t
n
n




, trong đó
h
e



.
Nếu nghiệm chính xác bị chặn thì 1
h
e


. Điều này chỉ có thể xảy ra nếu
0Re  hh
R

. Điều này có nghĩa là, trên mặt phẳng với trục hoành )Re( h

và trục
tung )Im( h



, mọi phương pháp Runge-Kutta (1.8)-(1.9) đều có thể viết dưới dạng
1
( )
i i
x R z x

 ,
trong đó
( )
R z
được gọi là hàm ổn định.
Định nghĩa 3.1
Tập tất cả các điểm của mặt phẳng phức
M
mà
( ) 1
R z

được gọi là miền ổn định
của phương pháp (1.8)-(1.9). Nếu tập hợp đó chứa toàn bộ nửa mặt phẳng trái thì
phương pháp được gọi là ổn định-A, còn nếu ngoài ra
lim ( ) 0
z
R z


thì phương pháp
được gọi là ổn định-L(hay còn gọi là ổn định tiệm cận).
1.3.2. Sự ổn định của phương pháp Euler

.
Phương pháp số là ổn định nếu
1

.
Xét các trường hợp sau
1)

là số thực. Khi ấy 11  h

, hay 02



h

.
2)

là thuần ảo (


i

, trong đó

là số thực khác 0). Khi ấy
111
22
 hhi
Hình 1.5

Như vậy, chỉ có một phần rất nhỏ (hình tròn bán kính bằng 1) của nửa mặt phẳng
trái là miền ổn định của phương pháp Euler.
Với mọi giá trị khác của
h

trong nửa mặt phẳng trái và bên ngoài hình tròn này,
nghiệm số sẽ bị phóng đại (blow-up) khi nghiệm chính xác triệt tiêu (decays).
Phương pháp số này được gọi là ổn định có điều kiện.
Để nhận được nghiệm số ổn định, bước
h
phải được chọn sao cho
h

nằm trong
hình tròn. Nếu

là số thực âm thì từ điều kiện 02



h

suy ra 0

1
1
z

,
1
2
1
2
z
z


và
1
z

.
Im
h

Re
h

-

O

h
hxhhhxkkxx










2
11
2
1
2
1
22
211


.
Để phương pháp ổn định thì
1

, trong đó
2
1
22


thuần ảo,
0


.
Khi ấy 1
4
1
1
44
 h

. Phương pháp không ổn định.
Trường hợp 3.
IR
i

 là số phức. Khi ấy
 
IRI
IR
R
ihh
h
h





và tìm nghiệm phức
h

của phương trình bậc hai
theo các giá trị của

. Nhận xét rằng
1

với mọi giá trị của

.
Miền ổn định được chỉ ra trên Hình 1.6.
1.3.3.2. Sự ổn định của phương pháp Runge-Kutta bậc bốn
Xét phương pháp Runge-Kutta bậc bốn cho phương trình thử (1.13). Ta có
nnn
hxtxhfk

 ),(
1
;
 
nnnn
x
h
hhxxh
k
xhk





















































42
1
42
1
332222
34






!
3
2
1
443322
hhh
h


 .
Trường hợp 1.

là số thực. Khi ấy 0785.2



h

.
Trường hơp 2.


i

thuần ảo, 0


. Khi ấy 220  h

.

.
Miền ổn định được chỉ ra trên Hình 1.6. Hình 1.6 19
1.3.4. Sự ổn định của phương pháp đa bước
Áp dụng các phương pháp (1.11) và (1.12) cho bài toán (1.13) ta được
 
1
0
0
k
j j i j
j
z x
 
 

 

. (1.14)
Phương trình đặc trưng của phương trình sai phân tuyến tính trên có dạng
 
0
0
k
k j

s
-nấc cho phương trình (1.13)-(1.15)
không vượt quá
2
s
(chắn Butcher).
2)
p
là cấp chính xác,
k
là số bước của phương pháp (1.11).Nếu phương pháp ổn
định thì
p
không vượt quá
2
k

khi
k
chẵn và
1
k

khi
k
lẻ (chắn Dahlquist thứ
nhất).
3) Phương pháp (1.11) ổn định – A không thể có cấp chính xác vượt quá 2 (chắn
Dahlquist thứ hai).
Trong Chương sau ta sẽ trình bày phương pháp do Bulatov đề nghị cải tiến được



x
thì nghiệm bị chặn. Đại lượng
kx
e

là nghịch biến khi
0

k
và
đồng biến khi
0

k
.
Xấp xỉ sai phân trung tâm cho hệ là
   
0
2
2
1
1111
2

 iiiii
xx
h
k

kh
nên
i
kh
kh








2
2
là nghịch biến. Khi 0

k thì
1
2
2



kh
kh
nếu
k
h
2

(0)
x x

, (2.2)
trong đó




, ,
f x t x t
là các hàm vectơ
n
- chiều, hàm
f
xác định trên hình hộp
chữ nhật vô tận


1,0: 
n
RD .
Để giải bài toán (2.1)-(2.2), ta bắt đầu đi từ phương pháp- và phương pháp một tựa
của nó (hai phương pháp này đều có cấp chính xác bằng một):

22


1 1 1 1 1
( , ) (1 ) ( ,

ta được:
      


2
1 1 1 1 1 1
, , ,
i i i i i i i i i i
f x t f x t J x t x x O x x
     
    

và




   
 
 
 
 
    
 
2 1 2 2
2
2 2 2 1 2 1
2
2 2 2 1 1
1 ,


1 1 1 1 1 1 1
( , ) (1 ) ( , ) ( , )( )
i i i i i i i i i i
x x h c f x t c f x t c J x t x x
   
     
, (2.3’)

1 2 2 2 1
( , ) ( , )( )]
i i i i i i i i
x x h f x t c h c J x t c h x x
 
     
. (2.4’)
Viết lại hai công thức trên thành một hệ phương trình đại số tuyến tính có ẩn là
1
i
x


(kí hiệu
E
là ma trận đơn vị cấp
n
):
1 1 1 1 1 1 1
1
2 2 2 2 2 2

   cấp
2
n n

, trong đó
3
c
là hằng số tùy
ý. Khi ấy (2.5) trở thành:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
2 2
1 1 3 2 2 1
2 2
1 1 3 2 2

, ,
c c c
và có bậc hội
tụ tối thiểu là bậc một. Ta có ba tham số tự do, vì vậy có thể chọn được bộ ba số
1 2 3
, ,
c c c
sao cho họ lược đồ sai phân (2.6) hội tụ cấp hai.
Khai triển Taylor theo
t
tại điểm


,
i i
x t
ta được:








     
 
     
 
     



 

  

   

  

   
 
    

Thay vào (2.6) ta được





   
 
   
 

 
3 1 1 2 3 2 3
2
2 2
2 2 2

   
 

     
 
 
     
 
2 2
2 2 2
1 3 2 2
1 1 1 1 1
2 2 3
3 2 2 2
,
, , , ,
( , ) , , , 1 ,
( , ) , , ,
t i i
i i t i i i i t i i i
i i t i i i i t i i i i
i i t i i i i t i i
J x t
h c J x t hJ x t c c J x t c hJ x t x
h E hc J x t hc J x t c f x t c hf x t c f x t
hc E hc J x t h c J x t f x t hc f x t O h
 

 
 

 
2
2 2 2 2
3 1 2 3 1 2 3 1 2 3
2 2 3
1 1 3 2 3 2
(1 ) 2 , 2 , ,
2
, , , , , ,
i i t i i i i i i
i i i i t i i i i i i t i i
h
c E h c c c J x t h c c c J x t c c c J x t hx x
hE h c J x t f x t hc f x t hc E h c c J x t f x t hc f x t O h
 
 
  
       
 
 
 
 
      
hay
       
 
 
 
             
       

 
3
, .
t i i
x t O h

Đẳng thức trên đúng với mọi
h
. So sánh hệ số hai vế ta được:
 Hệ số của
1
h
:
),(),()1(
33 iiiii
txfctxfxc 

 .
Do
( ) ( ( ), )
x t f x t t


(phương trình (1.1)) nên hệ thức trên đúng với mọi
3
c
.
 Hệ số của
2
h :

 
  
   
 
.
Do đó hệ số của
2
h
viết lại thành:
     
3
1 2 3
1
( , , ( , )) ( )( ( , ) ( , ) , )
2
i i i i t i i i i i i t i i
c
J x t f x t f x t c c c J x t f x t f x t

 
    
hay
   
 
3
1 2 3
1
, , ( , ) 0
2
i i i i t i i

Nếu










1,0,,,
2
CttxJttxf  thì ta có đánh giá


2
( )
i i
x x t O h
  , trong đó
i
x

tìm được theo công thức (2.6).
Định lý này được chứng minh nhờ nhận xét là sai số địa phương có bậc ba, và ma
trận chuyển sang bước tiếp theo có chuẩn là 1
hK

với 0 K

( , ) , ( , )
i i i i i i
f x t c h x f x t x
 
  
   .
Khi đó (1.6) được viết lại thành:
       
2 2 2 2
1 3 2 1 1 3 2
1 1 1 3 2
( )( (1 ) ) ( )
i i
i i i
E hc c E hc x E hc c E hc x
h E hc c x c x hc E hc x
   
    

   
       
   
     

hay


2 2 2
3 1 2 3 1 2 3 1
2 2 2

trong đó
 
2 2 2
3 1 3 2 3 1 1 3 2 2 3
2 2 2
3 1 2 3 1 2 3
(1 ) (1 2 2 ) ( )
(1 ) 2( ) ( )
c c c c c z c c c c c c z
R z
c c c c z c c c z
        

    
.
Để (2.6) là ổn định-L thì


lim 0
z
R z


, suy ra trong


R z
bậc của tử thức phải nhỏ
hơn bậc của mẫu thức, tức là ta phải chọn
1 2 3