1BỘ GIÁO DỤC VÀ ĐÀO TẠO
TRƯỜNG……………………

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


Chương 1 trình bày một số khái niệm và phương pháp cơ bản giải số phương
trình vi phân. Trong mục 1.2 của Chương, chúng tôi trình bày các phương pháp số
cổ điển theo một quan điểm nhất quán là xuất phát từ Quy tắc cầu phương cơ bản.
Chương 2 trình bày phương pháp không cổ điển (do Bulatov đề xuất vào
những năm 2003-2008) giải số hệ phương trình vi phân bậc nhất, phi tuyến và tuyến
tính, theo các tài liệu [9]-[11].
Chương 3 trình bày phương pháp không cổ điển giải số hệ phương trình vi
phân bậc hai, tuyến tính và phi tuyến, theo bài báo của M. V. Bulatov và G. V.
Berghe ([4], 2009).

3

Thông qua việc tính toán đạo hàm, phân tích các hàm nhiều biến vào chuỗi
Taylor và các phép biến đổi chi tiết, chúng tôi cố gắng trình bày các kết quả của M.
V. Bulatov và G. V. Berghe một cách rõ ràng và chi tiết nhất.
Để minh họa và kiểm chứng lý thuyết, chúng tôi đã lập trình trên MATLAB và
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ả




, ,
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 vô
hạn


n
RD  1,0:
.
Ở đâ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

g t
là vectơ
n
-chiều,
tức là hệ tuyến tính


( ) ( ) ( ), 0,1
x t B t x g t t

  
. (1.3)
Ta luôn giả thiết rằng các phần tử của ma trận
( )
B t
, của các vectơ


,
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


(1.4)
nên ta cũng có thể sử dụng quy tắc cầu phương cơ bản trong việc giải số phương
trình vi phân. Trong mục này ta sẽ chỉ ra rằng, nhiều công thức sai phân cổ điển giải
số phương trình vi phân có thể suy ra từ quy tắc cầu phương cơ bản. Trước tiên ta
nhắc lại quy tắc cầu phương cơ bản (xem, thí dụ, [1]).
Nội dung cơ bản của quy tắc cầu phương là: để tính tích phân
( )
b
a
f t dt

ta thay
( )
f t

bởi một đa thức nội suy (interpolating polynomial). Tích phân của hàm
( )
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



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

( )
k
P t
bậc nhỏ hơn
p

ta có:
1
( ) ( ).
b
s
k j j
j
a
P t dt f 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( ).

t h
t
x t h x t f x s s ds

  

(1.6)
Kết hợp với công thức (1.5) ta đi đến công thức:
( ) ( ) . ( ( ), )
x t h x t h f x t t
  
(1.7)
Gọi
h
là độ dài bước (stepsize) của biến độc lập
t
(
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


 
.
Đây chính là công thức Euler tiến quen thuộc. Hình 1.1
 Nếu chọn
1
s

và
1
c b

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 ABEF (Hình 1.2):
( ) ( ) ( ).
b
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

( ) ( ) . ( ( ), )
2 2
h h
x t h x t h f x t t
    
.

O

f

x

C

D

b

B

a

A

E

b

B


f

x8

Từ công thức trên ta có
1
[ ( ( ), ]
2 2
n n n n
h h
x x h f x t t

   
.
Đâ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

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

và
2
2 2
1 ( )
( ) .
( ) ( ) 2 2
a
b b
a a
b
t a t a b a
L t dt dt
b a b a

  
   
 
 

Chứng tỏ
( ) [ ( ) ( )]
2

Từ đây ta có công thức hình thang:
1 1 1
[ ( , ) ( , )]
2
n n n n n n
h
x x f x t f x t
  
   .
Phương pháp điểm giữa và phương pháp
hình thang là hai phương pháp ẩn, Hình 1.4
chúng có độ chính xác
2
p

.
 Nếu chọn
3
s

và
1 2 3
, ,
2
a b
c a c c b

  
thì, đặt
h b a

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

2 2
3 2
2
2 2
3
2
2 2
( ) ( )( ) ( )( )
2 2
2 2 ( ) ( )
[( ) ( )( )] [ ( ) ]
2 3 2 2
2 ( )
.
12 6
b b b
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


        
  
   
  

Do tính chất đối xứng (hoặc tính trực tiếp), ta có
3 1
6
h
 
 
.
Từ các tính toán trên ta đi đến công thức Simpson:
( ) [ ( ) 4 ( ) ( )]
6 2
b
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

ta thay giá trị
1
n
x

ở vế phải bằng công thức Euler tiến:
1
ˆ
( , )
n n n n
x x hf x t

  .
Khi ấy ta được công thức:
1 1 1
ˆ
[ ( , ) ( , )]
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

ˆ
( , )
2
n n n
n
h
x x hf x t


  

 Từ phương pháp Runge-Kutta ẩn cấp bốn kinh điể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
  
     
,
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



1.2.2.2. Phương pháp Runge-Kutta tổng quát
Nội dung cơ bản của phương pháp Runge-Kutta tổng quát như sau.
Chia đoạn


0,1
thành một lưới đều
i
t ih

,
0,1,2, ,
i N

,
1
h
N

,
và kí hiệu
i
x
là giá trị xấp xỉ


i
x t
,









s
j
iijjiii
hctXahxfX
1
,
. (1.9)
Các tham số
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

X

(mỗi vectơ
i
X
có
n
tọa độ).
Thường phương pháp Runge-Kutta được viết dưới dạng bảng Butcher (Butcher
table)
c

AT
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

(tại các
điểm khác nhau trong công thức Runge-Kutta) là
s
.
Để tìm các phương pháp Runge-Kutta bậc hai, ta làm như sau (xem [2]).
Khai triển Taylor hàm
( , )
f x t
theo phương trình (1.1) và theo công thức (1.10) rồi
so sánh, ta đi đến kết luận:
Các hệ số trong phương pháp Runge-Kutta cấp hai phải thoả mãn hệ phương trình
1 2 2 2 21 2
1 1
1, ,
2 2
b b c b a b
   
.

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

2
b

và
21 2
1
a c
 
. Khi ấy ta có một phương pháp Runge-Kutta
cấp hai cho phép tính
1
n
x

dựa trên công thức:
1
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

1.2.3. Phương pháp cổ điển đa bước
Phương pháp cổ điển
k
-bước cho bài toán (1.1) có dạng (xem, [3], [4]-[7])
1 1 1
0 0
( , )
k k
j i j j i j i j
j j
x h f x t
 
     
 

 
. (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
  
     
  

  

thì ta có phương pháp ẩn. 14
1.3. Mô hình thử và ổn định của phương pháp số
1.3.1. Mô hình thử
Để phân tích hiệu quả của các phương pháp, ta thường thử chúng trên mô hình G.
Dahlquist (gọi là phương trình thử hay mô hình thử)




0
, 0 , 0,1
x x x x R t


   
, (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

xx
n
n


. Ta thấy rằng nghiệm này bị chặn khi và chỉ khi
1

.
Giả sử bước 0

h cố định. Khi ấy giá trị của nghiệm chính xác tại các điểm nht
n

sẽ là
000
xxexex
nnh
t
n
n




, trong đó
h
e




là thuần ảo và
1

thì ổn định tuyệt đối được gọi là ổn định tuần hoàn (P-
ổn định).
Khi miền ổn định của phương trình sai phân đồng nhất với miền ổn định của
phương trình vi phân, lược đồ sai phân hữu hạn được gọi là ổn định - A.

15
Phương trình thử thường được sử dụng như một mô hình để dự đoán tính ổn định
của phương pháp số giải hệ dạng tổng quát (1.1)-(1.2).
Để thuận tiện, ta cũng có thể đưa ra các khái niệm ổn định tương tự như sau.
Kí hiệu
z 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

0
1 xxhx
n
n
n
n


,
trong đó
h




1
.
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




22
 hhhihh
IRIR

,

16
nghĩa là
h

nằm trong hình tròn đơn vị tâm là


0;1
(Hình1.5). Hình tròn này tiếp
xúc với trục ảo. 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).

âm và có trị tuyệt đối lớn hơn 1. Vì


0
1 xhx
n
n

 nên nghiệm số sẽ đổi dấu qua
mỗi bước. Sự thay đổi của nghiệm số mô tả khá rõ tính không ổn định.
Tương tự, ta có thể xét tính ổn định của các phương pháp Euler cải tiến. Ta đi đến
kết luận sau.
Phương pháp Euler ẩn là ổn định-L và hội tụ cấp một; Phương pháp hình thang là
ổn định-A và hội tụ cấp hai, còn phương pháp Euler hiển không phải là ổn định-A
và hội tụ cấp 1.
Hàm ổn định của các phương pháp này tương ứng là
1
1
z

,
1
2
1
2
z
z


và







nnnnnn
xhhhxxhkxhhtkxhfk

 1),(
112

và
    
nnnnn
x
h
hxhhhxkkxx










2
11

h
h


hay
02



h

.
Trường hơp 2.


i

thuần ảo,
0


.
Khi ấy 1
4
1
1
44
 h

. Phương pháp không ổn định.

22
2
1
2
1
là số phức. Đặt



i
e
h
h 
2
1
22
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

2
1
2


;

18
nnnn
x
hh
hx
h
hxh
k
xhk
















;
 
nnnn
x
hh
hhx
hh
hxhkxhk
























!4!32
122
6
1
443322
43211


.
Để phương pháp ổn định thì 1

, trong đó
!
4
!
3
2
1
443322
hhh
h


 .
Trường hợp 1.

là số thực. Khi ấy 0785.2


!
3
2
1
443322
và tìm
nghiệm phức
h

của phương trình bậc bốn 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. 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

z


với mọi nghiệm của (1.15)
và
( ) 1
z


đối với các nghiệm bội được gọi là miền ổn định của phương pháp
(1.14). 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.
Nhận xét
Với mọi phương pháp cổ điển miền ổn định chứa gốc tọa độ của mặt phẳng phức.
Tương ứng với bậc của phương pháp ta đã biết những điều sau (xem [11]):
1) Bậc của các phương pháp Runge-Kutta
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á

Phương trình có nghiệm là
kx
BeAtx

)(
, trong đó
A
và
B
là các hằng số bất kì
được xác định bởi điều kiện ban đầu hoặc điều kiện biên tương ứng.
Nếu
0

k
và


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





2
2
.
Nếu 0

k đại lượng 1
2
2



kh
kh
nên
i
kh
kh








2

sáng tỏ phương pháp, chúng tôi đã thực hiện các tính toán chi tiết (phân tích các
hàm nhiều biến dưới dạng chuỗi Taylor, ) mà trong [9]-[11] trình bày không tường
minh.
2.1. Phương pháp không cổ điển giải số hệ phương trình vi phân
phi tuyến cấp một
2.1.1. Phương pháp tổng quát
Xét bài toán Cauchy tìm nghiệm của hệ phương trình


( ) ( ( ), ), 0,1
x t f x t t t

 
(2.1)
thỏa mãn điều kiện ban đầu
0
(0)
x x

, (2.2)
trong đó




, ,
f x t x t
là các hàm vectơ
n
- chiều, hàm

,
c c
là các hằng số tùy ý,
1 2
0 , 1
c c
 
.
Mỗi công thức truy hồi (2.3) (hoặc (2.4)) cho một dãy các giá trị xấp xỉ nghiệm của
phương trình vi phân (2.1)-(2.2). Dưới đây ta cố gắng kết hợp hai phương pháp
(2.3) và (2.4) để được một phương pháp số mới giải hệ phương trình vi phân (2.1)-
(2.2).
Khai triển Taylor tại điểm
i
x x

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à



trong đó
 


,
,
f x t
J x t
x



.
Các phương pháp (2.3)-(2.4) được tuyến tính hóa như sau


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.5)
Nhận thấy rằng hệ (2.5) nói chung không có nghiệm theo nghĩa cổ điển vì số
phương trình nhiều hơn số ẩn, tức là hệ (2.3) và (2.4) nói chung không có nghiệm
trùng nhau. Để giải hệ phương trình đại số (2.5) ta nhân hai vế của hệ này với ma

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

 

 
 
   
 
 
 
    
 
 
   
   
(2.6)
Theo [11]lược đồ sai phân (2.6) là ổn định với mọi bộ hệ số
1 2 3
, ,
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


,

, , ,
, , ,
, , ,
, , ,
2
i i i i t i i
i i i i t i i
i i i i t i i
i i i i t i i
i i i i i
J x t J x t J x t h O h
J x t c h J x t J x t c h O h
f x t f x t f x t h O h
f x t c h f x t f x t c h O h
h
x x t x hx x O h


 

  

   

  

   
 
    


       
 
 
 
   
     
 
 
 
 
 
      
 

   
 
   
 

     
 
 
     
 
2 2
2 2 2
1 3 2 2
1 1 1 1 1
2 2 3
3 2 2 2


hay

24
       
 
 
 
 
 
   
 
 
 
   
 
 
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

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
hf x t c h f x t ch J x t f x t c h J x t f x t hc f x t
c c h f x t c c h J x t f x t c c h J x t f
 
 
  
       
 
 
 
 
     
 
  
 
 
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(

 
   

Từ phương trình (1.1) ta có:
( ) ( ( ), )
x t f x t t



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

3 1 2 3
1 2 2
c c c c
   . (2.7)

25
Như vậy nếu chọn
1 2 3
, ,
c c c
thỏa mãn (2.7) thì ta được công thức (2.6) có cấp hai,
bởi vì lúc này công thức (2.6) có sai số địa phương bậc
3
h
. Ta đi đến định lý sau.
Định lý 1.1
Nếu










1,0,,,
2
CttxJttxf  thì ta có đánh giá

. (2.8)
Ta có


, ,
f x t x

 do đó


, , ,
J x t x t

  
và
2 1 1 1
( , ) , ( , )
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

z h


và nhận xét rằng, với phương trình thử, ma trận
1
E

là một số, ta được
công thức
1
( ).
i i
x R z x


,
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
        

    