Xây dựng trình thực thi cho phương pháp sai phân lùi dạng khối liên tục

 ISSN: 1859-2171 
e-ISSN: 2615-9562 
TNU Journal of Science and Technology 225(06): 424 - 431 
424  Email: jst@tnu.edu.vn 
XÂY DỰNG TRÌNH THỰC THI CHO PHƯƠNG PHÁP 
SAI PHÂN LÙI DẠNG KHỐI LIÊN TỤC 
Đinh Văn Tiệp*, Phạm Thị Thu Hằng
Trường Đại học Kỹ thuật Công nghiệp - ĐH Thái Nguyên 
TÓM TẮT 
Phương pháp đa bước sai phân lùi dạng khối là một cải tiến nổi trội để tìm xấp xỉ nghiệm của 
phương trình vi với điều kiện ban đầu. Phương pháp này với miền ổn định tuyệt đối lớn nên đặc 
biệt thích hợp cho lớp các bài toán stiff. Bài báo này trình bày cách xây dựng trình thực thi cho 
phương pháp này bằng cách sử dụng phương pháp lặp Newton cho khối. Một chương trình bằng 
Ngôn ngữ Matlab được đưa ra để mô tả thuật toán và trình thực thi với số bước k cụ thể. Các ví dụ 
trình bày đưa ra sự so sánh về tính hiệu quả và chính xác của phương pháp. 
Từ khóa: phương pháp tuyến tính đa bước; phương pháp sai phân lùi; phương pháp sai phân lùi 
dạng khối liên tục; bài toán stiff; miền ổn định tuyệt đối. 
Ngày nhận bài: 19/5/2020; Ngày hoàn thiện: 27/5/2020; Ngày đăng: 29/5/2020 
CONSTRUCTING THE IMPLIMENTATION 
TO THE CONTINUOUS BLOCK BDF METHODS 
Dinh Van Tiep*, Pham Thi Thu Hang
TNU - University of Technology 
ABSTRACT 
The k-step methods of continuous block backward difference formula (or BDF) have great benefit 
to approximate the differential equation with the initial condition. These methods possessing very 
large absolute stability regions are especially useful for the stiff equations. This article presents a 
method of implementation to these k-step scheme by using the Newton’s iteration for the root 
finding problem of the non-linear multivariable case. Also, a program written in Matlab with a 
particular k is presented. The numerical experiments introduced to illustrate the efficiency and 
exactness of these scheme comparing with some conventional BDF methods. 
Keywords: linear multistep method; backward difference formula; continuous block backward 
difference formula; stiffness; absolute stability region. 
Received: 19/5/2020; Revised: 27/5/2020; Published: 29/5/2020 
* Corresponding author. Email: tiepdinhvan@gmail.com 
Đinh Văn Tiệp và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ ĐHTN 225(06): 424 - 431 
 Email: jst@tnu.edu.vn 425 
1. Mở đầu 
Xem xét bài toán giá trị ban đầu 
𝑦′ = 𝑓(𝑡, 𝑦), 𝑎 ≤ 𝑡 ≤ 𝑏, 𝑦(𝑎) = 𝛼, (1) 
với giả thiết điều kiện Lipschitz theo 𝑦 của 
hàm 𝑓. Các phương pháp xấp xỉ nghiệm của 
phương trình bao gồm chủ đạo các phương 
pháp đơn bước và đa bước. Họ các phương 
pháp sai phân lùi (ký hiệu BDF) với số bước k 
nhỏ, 𝑘 < 6, thể hiện ưu điểm lớn với miền ổn 
định tuyệt đối gần chiếm trọn mặt phẳng 
phức. Thực tế các miền đó của các phương 
pháp này chứa toàn bộ phần âm của trục thực 
trên mặt phẳng phức. Điều đó cho thấy khả 
năng áp dụng của các phương pháp này cho 
các bài toán stiff, một trong những vấn đề gây 
khó khăn lớn đối với đa số các phương pháp 
xấp xỉ khác. 
Tuy nhiên, một trong những trở ngại cho tính 
áp dụng của các phương pháp họ sai phân lùi, 
cũng giống như các phương pháp đa bước bậc 
cao khác, chẳng hạn các phương họ Adams, 
đó là chúng đều ở dạng ẩn. Cách tiếp cận để 
xử lý vấn đề này gồm các hướng sử dụng 
dạng dự báo - hiệu chỉnh và hướng sử dụng 
các phép lặp cho bản thân mỗi phương pháp. 
Với hướng tiếp cận dự báo - hiệu chỉnh, miền 
ổn định tuyệt đối khá bé. Không may, điều 
này đặc biệt đúng với các cặp dự báo - hiệu 
chỉnh với kiểu hiệu chỉnh họ sai phân lùi. Khi 
đó, các ưu điểm nổi trội của họ sai phân lùi bị 
triệt tiêu. Đối với hướng sử dụng kiểu lặp có 
hạn chế ở việc nghiệm hội tụ đến nghiệm của 
phương trình sai phân, thay vì đến nghiệm 
của phương trình vi phân. 
Hướng tiếp cận sai phân lùi dạng khối liên tục 
được đưa ra bởi N.M. Yao và cộng sự với bậc 
thấp 𝑘 ≤ 6 (xem [1], [2]) thể hiện sự nổi trội 
về tính áp dụng đối với lớp các bài toán stiff. 
Đặc biệt hơn, khả năng tự khởi động của 
phương pháp tạo nên sự thống nhất giữa khớp 
nối mà nhiều phương pháp đa bước gặp phải. 
Trình thực thi của phương pháp khối liên tục 
có hai hướng tiếp cận chủ đạo, một là sau 
bước khởi động, không dùng quy trình lặp để 
giảm khối lượng tính toán, hai là các bước lặp 
khối sẽ tiến hành cho đến cuối cùng. Với 
hướng tiếp cận thứ hai, độ chính xác tốt hơn, 
song yêu cầu số bước chia phải là một bội 
nguyên lần số bước của phương pháp. 
2. Công thức sai phân lùi dạng khối liên tục 
Trong [1], [2], N. M. Yao và cộng sự đã trình 
bày quy trình xây dựng phương pháp tổng 
quát, và đưa ra các công thức khối cụ thể cho 
các trường hợp 𝑛 = 2,3,4,6. Đồng thời, trong 
các bài báo này, các tác giả cũng đề cập đến 
vấn đề về tính ổn định, bậc và các hệ số bậc 
của từng công thức khối. Đối với các công 
thức 𝑛 = 2,3, các công thức khối dạng này là 
A-ổn định (tức A-stable). Với 𝑛 = 4,6, đặc 
tính đó có yếu đi chút ít song chúng vẫn là 
𝐿0-ổn định (tức 𝐿0-stable) theo định nghĩa của 
Cash [3]. Các tính chất miền tuyệt đối này 
được bảo tồn từ các phương pháp BDF thông 
thường. Như vậy, các ưu điểm vượt trội của 
phương pháp BDF không bị mất đi với 𝑘 ≤ 6. 
Ở đây, khác so với trình thực thi giới thiệu 
trong [1], ta sẽ xây dựng quy trình lặp tổng 
quát cho khối được tạo ra theo phương pháp 
này với chứng minh sự hội tụ của nghiệm xấp 
xỉ khi ℎ → 0+. 
Xét công thức dạng khối liên tục xây dựng 
kiểu sai phân lùi [1]: 
𝐴(1)𝑌𝑛+1 = 𝐴
(0)𝑌𝑛 + ℎ𝐵
(1)𝐹𝑛+1, (2) 
với 𝑌𝑛+1 = (𝑦𝑛+1, 𝑦𝑛+2,  , 𝑦𝑛+𝑘)
𝑇 , 
𝑌𝑛 = (𝑦𝑛−𝑘+1, 𝑦𝑛+𝑘 ,  , 𝑦𝑛)
𝑇 , 
𝐹𝑛+1 = (𝑓𝑛+1, 𝑓𝑛+2,  , 𝑓𝑛+𝑘)
𝑇 , 
𝑦𝑗 ≈ 𝑦(𝑡𝑗), 𝑓𝑗 ≈ 𝑓 (𝑡𝑗, 𝑦(𝑡𝑗)) , ∀𝑗 = 1,2, , 
𝐴(1), 𝐴(0), 𝐵(1) là các ma trận cấp 𝑘 × 𝑘 thu 
được từ xây dựng phương pháp. Ở đây, giả sử 
phương pháp BDF cho bởi công thức: 
ℎ𝑓𝑛+𝑖 = ∑ 𝑎𝑖𝑗𝑦𝑛+𝑗
𝑘−1
𝑗=0
+ ℎ𝑏𝑖𝑓𝑛+𝑘, 
∀𝑖 = 1, . . . , 𝑘 − 1. (3) 
𝑦𝑛+𝑘 = ∑ 𝑎𝑗𝑦𝑛+𝑗
𝑘−1
𝑗=0
+ ℎ𝑏𝑘𝑓𝑛+𝑘, 
Đinh Văn Tiệp và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ ĐHTN 225(06): 424 - 431 
426  Email: jst@tnu.edu.vn 
𝐴(1) có cột 𝑘 là (0, . . ,1)𝑇 , và hàng 𝑘 là 
−(𝑎1, . . , 𝑎𝑛−𝑘+1, −1), cột thứ 𝑗 của ma trận là 
−(𝑎1𝑗, 𝑎2𝑗,  , 𝑎𝑗)
𝑇
; ma trận 𝐴(0) có tất cả các 
cột bằng véctơ không, ngoại trừ cột cuối cùng 
là véctơ (𝑎10,  , 𝑎𝑘−1,0, 𝑎0)
𝑇
; ma trận 𝐵(1) 
có dạng tam giác trên với các phần tử khác 
không nằm trên đường chéo chính và cột thứ 
𝑘, trong đó đường chéo chính là véctơ 
(−1, , 𝑏𝑘), và cột 𝑘 là véctơ (𝑏1,  , 𝑏𝑘)
𝑇 . 
Các ma trận 𝐴(1) cho 𝑘 = 2,3,4,5,6 đều là các 
ma trận khả nghịch. Do đó, phương trình (2) 
có thể đưa được về dạng: 
𝑌𝑛+1 = 𝐶𝑌𝑛 + ℎ𝐵𝐹𝑛+1, (4) 
với 𝐶 là ma trận cấp 𝑘 × 𝑘 có các phần tử 
khác 0 ngoại trừ có thể các phần tử cột k, ký 
hiệu là: 
𝐴 = (𝑐1, 𝑐2,  , 𝑐𝑘)
𝑇 . (4′) 
Trong (3), 𝐹𝑛+1 = 𝐹(𝑡𝑛+1, 𝑌𝑛+1) 
= (𝑓(𝑡𝑛+1, 𝑦𝑛+1), , 𝑓(𝑡𝑛+𝑘 , 𝑦𝑛+𝑘)), 
phụ thuộc vào 𝑌𝑛+1. Chính vì thế, phương 
trình (4) có thể xấp xỉ nghiệm bằng phương 
pháp lặp Newton hoặc Newton suy rộng với 
giả thiết hàm 𝑓 đủ khả vi. Quá trình lặp 
Newton được áp dụng cho phương trình 
𝐺(𝑋) = 0 với mỗi 𝑞 ≥ 0 là: 
𝐺(𝑋) = 𝑦𝑞𝑘𝐴 + ℎ𝐵𝐹(𝑡𝑞𝑘+1, 𝑋) − 𝑋, (5) 
𝑋 = (𝑥1,  . , 𝑥𝑘)
𝑇 ∈ ℝ𝑘 , để tìm nghiệm xấp 
xỉ của phương trình này. Thực tế, khi ℎ → 0, 
nghiệm của (5) hội tụ về nghiệm đúng: 
𝑌(𝑡𝑞𝑘+1) = (𝑦(𝑡𝑞𝑘+1), , 𝑦(𝑡𝑞𝑘+𝑘))
𝑇
, 
khi giả thiết 𝑦𝑞𝑘 = 𝑦(𝑡𝑞𝑘). Điều này được 
chứng minh trong Định lý 1 dưới đây. Bậc và 
hệ số bậc của phương pháp khối là bậc nhỏ 
nhất và hệ số bậc nhỏ nhất để xấp xỉ các 
𝑦(𝑡𝑞𝑘+𝑖), ∀𝑖 = 1, , 𝑘, bằng các phương trình 
sai phân thứ 𝑖 tương ứng trong (3). Công thức 
lặp nghiệm cho (5) (tại bước thứ q) là: 
𝑋(𝑚+1) = 𝑋(𝑚) − [𝐽𝐺(𝑋
(𝑚))]
−1
𝐺(𝑋(𝑚)), (6) 
với 
𝑋(0) = {
(𝑦(𝑞−1)𝑘+1,  , 𝑦𝑞𝑘) nếu 𝑞 ≥ 1
(0, , 𝑦0), 𝑦0 = 𝛼, nếu 𝑞 = 0.
Jacobian của G là 𝐽𝐺(𝑋) = 
ℎ𝐵𝑑𝑖𝑎𝑔 (𝑓𝑦(𝑡𝑞𝑘+1, 𝑥1), , 𝑓𝑦(𝑡𝑞𝑘+𝑘 , 𝑥𝑘))
− 𝐼𝑘 , 
 trong đó 
𝑑𝑖𝑎𝑔 (𝑓𝑦(𝑡𝑞𝑘+1, 𝑥1),  , 𝑓𝑦(𝑡𝑞𝑘+𝑘 , 𝑥𝑘)) 
là ma trận đường chéo với các phần tử này 
xếp dọc đường chéo chính và 𝐼𝑘 là ma trận 
đơn vị cấp 𝑘. 
Giả sử số đoạn chia của [𝑎, 𝑏] là 𝑁 = 𝑘𝑝, tức 
𝑁 là bội nguyên lần của 𝑘. Từ điều kiện ban 
đầu 𝑦0 = 𝛼, với 𝑞 = 0, áp dụng (5) ta được 
𝑌1 = (𝑦1, 𝑦2,  , 𝑦𝑘)
𝑇 . 
Tiếp tục với 𝑞 = 1, áp dụng (5) thu được: 
𝑌𝑘+1 = (𝑦𝑘+1, 𝑦𝑘+2,  , 𝑦2𝑘)
𝑇 . 
Cứ như thế, với 𝑞 = (𝑝 − 1), thu được: 
𝑌𝑘(𝑝−1)+1 = (𝑦𝑘(𝑝−1)+1, 𝑦𝑘(𝑝−1)+2,  , 𝑦𝑁)
𝑇
. 
Như vậy, sau quá trình thu được dãy véctơ 
𝑌1, 𝑌𝑘+1,  , 𝑌𝑘(𝑝−1)+1 
xấp xỉ nghiệm của (1). Mệnh đề dưới đây chỉ 
ra một đặc trưng của trình thực thi ta đang 
xây dựng. 
Mệnh đề 1. Ma trận 𝐴 trong (4′) có tất cả các 
thành phần là 1, tức 𝐴 = (1,1, ,1)𝑇 . 
 Chứng minh. Từ (3), ta được: 
𝐴(1)
=
[
−𝑎11 −𝑎12 ⋯
−𝑎21 −𝑎22 ⋯
⋮ ⋮ ⋱
−𝑎1,𝑘−1 0
−𝑎2,𝑘−1 0
⋮ ⋮
−𝑎𝑘−1,1 −𝑎𝑘−1,2 ⋯
−𝑎1 −𝑎2 ⋯
−𝑎𝑘−1,𝑘−1 0
−𝑎𝑘−1 1]
, 
𝐴(0) = [
0 ⋯
⋮ ⋱
0
⋮
𝑎10
⋮
⋮ ⋱ ⋮ 𝑎𝑘−1,0
0 ⋯ 0 𝑎0
] , 𝑎0 = 1, 
Đinh Văn Tiệp và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ ĐHTN 225(06): 424 - 431 
 Email: jst@tnu.edu.vn 427 
𝐵(1) = 
[
−1 0 ⋯
0 −1 ⋱
⋮ 0 ⋱
0 𝑏1
⋮ 𝑏2
0 ⋮
0 ⋮ ⋱
 0 0 ⋯
−1 𝑏𝑘−1
0 𝑏𝑘 ]
. 
Viết ma trận 𝐴(1), [𝐴(1)]
−1
 dưới dạng khối: 
𝐴(1) = [
𝐷 𝟎
𝑈 1
] , [𝐴(1)]
−1
= [𝐷
−1 𝟎
𝑉 1
], 
trong đó, 𝐷−1 là nghịch đảo của 𝐷, 
𝐷 = (−𝑎𝑖𝑗)(𝑘−1)×(𝑘−1), 𝐷
−1 ∈ 𝕄(𝑘 − 1), 
𝟎 = (0, ,0)𝑇 ∈ 𝕄(𝑘 − 1,1), 
𝑈 = (−𝑎1,  , −𝑎𝑘−1), 𝑉 ∈ 𝕄(1, 𝑘 − 1). 
Giả sử 𝑉 = (𝑥1,  , 𝑥𝑘−1)
𝑇 và 𝐷−1 có các cột 
lần lượt là 𝐷1, 𝐷2,  , 𝐷𝑘−1, hay viết là: 
𝐷−1 = [𝐷1, 𝐷2,  , 𝐷𝑘−1]. 
Ký hiệu 𝐸𝑙 là véctơ cột đơn vị thứ 𝑙 của 
𝕄(𝑘 − 1,1) = ℝ𝑘−1, ∀𝑙 = 1, , 𝑘 − 1, 
𝐸𝑙 = (0, ,0,1,0, . . . ,0) 
với số 1 ở vị trí thứ 𝑙. Gọi 𝐼𝑘−1 là ma trận đơn 
vị cấp 𝑘 − 1. Ta có: 
𝐷−1𝐷 = 𝐼𝑘−1 = [𝐸1,  , 𝐸𝑘−1], 
[𝑉 1] [
𝐷
𝑈
] = (0, ,0) ∈ 𝕄(1, 𝑘 − 1) = ℝ𝑘−1, 
nên 
∑(−𝑎𝑗𝑙)𝐷
𝑗
𝑘−1
𝑗=1
= 𝐸𝑙 , ∀𝑙 = 1, , 𝑘 − 1, (7) 
∑ 𝑎𝑗𝑙𝑥𝑗
𝑘−1
𝑗=0
+ 𝑎𝑙 = 0, ∀𝑙 = 1, , 𝑘 − 1. (8) 
Vì 𝑘 phương trình sai phân ở (3) đều có bậc 
𝑟 > 0 nên hệ số bậc 0 của các phương trình 
này là 0, tức ta có: 
∑ 𝑎𝑗𝑙
𝑘−1
𝑙=0
= 0, ∀𝑗 = 1, , 𝑘 − 1, (7′) 
∑ 𝑎𝑗
𝑘−1
𝑗=0
= 1. (8′) 
Kết hợp (7) và (7′) được: 
∑ 𝑎𝑗0𝐷
𝑗
𝑘−1
𝑗=1
= ∑ (− ∑ 𝑎𝑗𝑙
𝑘−1
𝑙=1
)𝐷𝑗
𝑘−1
𝑗=1
= ∑ (∑(−𝑎𝑗𝑙)
𝑘−1
𝑗=1
𝐷𝑗) = ∑ (∑ 𝑎𝑗𝑙
𝑘−1
𝑗=1
𝐷𝑗)
𝑘−1
𝑙=1
𝑘−1
𝑙=1
= ∑ 𝐸𝑙
𝑘−1
𝑙=1
= (1, ,1)𝑇 ∈ 𝕄(𝑘 − 1,1) = ℝ𝑘−1. 
Đẳng thức này là 𝑘 − 1 thành phần đầu tiên 
của véctơ 𝐴 = (1, ,1)𝑇 ∈ ℝ𝑘. 
Kết hợp (8) và (8′) được: 
∑ 𝑎𝑗0𝑥𝑗
𝑘−1
𝑗=1
= ∑ (− ∑ 𝑎𝑗𝑙
𝑘−1
𝑙=1
)𝑥𝑗
𝑘−1
𝑗=1
= 
∑ (− ∑ 𝑎𝑗𝑙
𝑘−1
𝑗=1
𝑥𝑗) = − ∑ 𝑎𝑙
𝑘−1
𝑙=1
= 𝑎0 − 1
𝑘−1
𝑙=1
= 0. 
Đẳng thức này là thành phần thứ 𝑘 của véctơ 
𝐴 = (1, ,1)𝑇 ∈ ℝ𝑘. □ 
Định lý sau khẳng định sự hội tụ (khi ℎ →
0+) của nghiệm xấp xỉ của phương trình (5) 
về nghiệm của phương trình vi phân (1). Điều 
này cho thấy, dù phương trình phi tuyết (5), 
tại mỗi giá trị 𝑞 ≥ 0, có thể có nhiều nghiệm, 
tuy nhiên với kích thước bước ℎ đủ nhỏ, việc 
lựa chọn 𝑦𝑞𝑘 và véctơ 
𝑋(0) = 𝑌(𝑞−1)𝑘+1
= (𝑦(𝑞−1)𝑘+1, 𝑦(𝑞−1)𝑘+2,  , 𝑦𝑞𝑘) 
như là xấp xỉ ban đầu của mỗi quá trình lặp 
Newton (6) là đúng đắn. Tức là nghiệm của 
phương trình phi tuyến (5) thu được từ sự hội 
tụ của dãy (6), với mỗi 𝑞 ≥ 0, có thể gần tùy 
ý đến nghiệm đúng của (1), chỉ cần điều chỉnh 
giá trị ℎ đủ nhỏ. 
Định lý 1. Cố định mỗi 𝑞 ≥ 0. Giả sử nghiệm 
thu được ở bước thứ 𝑞 − 1 là nghiệm đúng, 
tức 𝑌𝑞𝑘+1 = 𝑌(𝑡𝑞𝑘+1) nếu 𝑞 ≥ 1 và 𝑦𝑞𝑘 =
𝑦(𝑡𝑞𝑘) nếu 𝑞 = 0. Gọi 𝑋 = 𝑌∗ là nghiệm 
đúng của phương trình (5), 𝐺(𝑋) = 0, thu 
được từ sự hội tụ của dãy (6), lim
𝑚→∞
𝑋(𝑚) = 𝑌∗, 
Đinh Văn Tiệp và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ ĐHTN 225(06): 424 - 431 
428  Email: jst@tnu.edu.vn 
tức 𝑌∗ phụ thuộc 𝑌𝑞𝑘 . Gọi 𝑌
∗ = 𝑌(𝑡𝑞𝑘+1) =
(𝑦(𝑡𝑞𝑘+1), , 𝑦(𝑡𝑞𝑘+𝑘))
𝑇
là nghiệm đúng của 
phương trình vi phân (1). Khi đó, 
lim
ℎ→0+
||𝑌∗ − 𝑌
∗|| = 0. 
 Chứng minh. Giả sử 𝝆𝒌(ℎ)ℎ
𝑠 là véc tơ có 
các thành phần là sai số làm tròn (truncation 
error) của các phương trình sai phân lần lượt 
trong (3) với s là số lớn nhất sao cho ℎ𝑠 là 
nhân tử chung của 𝑘 thành phần này. Vì các 
hệ số bậc 0 và 1 của những phương trình này 
bằng 0, nên 𝑠 ≥ 2. Ta có: 
𝑌∗ = 𝑦(𝑡𝑞𝑘)𝐴 + ℎ𝐵𝐹(𝑡𝑞𝑘+1, 𝑌
∗) + 𝝆𝒌(ℎ)ℎ
𝑠, 
𝑌∗ = 𝑦𝑞𝑘𝐴 + ℎ𝐵𝐹(𝑡𝑞𝑘+1, 𝑌∗). 
Do đó, ‖𝑌∗ − 𝑌∗‖ ≤ ℎ‖𝐵‖‖𝐹(𝑡𝑞𝑘+1, 𝑌
∗) −
𝐹(𝑡𝑞𝑘+1, 𝑌∗)‖ + 𝝆𝒌(ℎ)ℎ
𝑠 ≤ ℎ𝐿𝐹‖𝐵‖‖𝑌
∗ −
𝑌∗‖ + 𝝆𝒌(ℎ)ℎ
𝑠. (Ở đây, sử dụng điều kiện 
Lipschitz của 𝑓 suy ra sự tồn tại của hằng số 
Lipschitz 𝐿𝐹 theo 𝑌 cho hàm véc tơ 𝐹(𝑡, 𝑌)) 
Từ đó, với ℎ > 0 đủ nhỏ, sao cho 
1
2
< (1 − ℎ𝐿𝐹‖𝐵‖), 
‖𝑌∗ − 𝑌∗‖ ≤
𝝆𝒌(ℎ)ℎ
𝑠
1 − ℎ𝐿𝐹‖𝐵‖
< 2𝝆𝒌(ℎ)ℎ
𝑠. 
Do đó, ta có điều phải chứng minh. 
Trên đây là hướng tiếp cận thứ hai của trình 
thực thi như đã đề cập ở đoạn cuối phần mở 
đầu bài báo. Hướng tiếp cận thứ nhất chỉ xét 
đến quá trình lặp Newton một lần để xấp xỉ 𝑌1 
phần còn lại sẽ sử dụng đồng thời phương 
trình thứ k của (3) như là phép hiệu chỉnh, 
trong đó vế phải có chứa 𝑓𝑘+1 mà giá trị này 
sẽ được tính toán bằng cách sử dụng phương 
trình liền trước nhất thứ 𝑘 − 1, và phương 
trình này được coi như phép dự báo. Như thế, 
ở hướng tiếp cận thứ 2 này, ta có một cặp dự 
báo - hiệu chỉnh. Tuy nhiên, cặp này có thể 
làm hẹp đi đáng kể độ lớn của miền ổn định 
tuyệt đối mà các phương pháp BDF bậc thấp 
(𝑘 ≤ 6) vốn có. 
Đối với 𝑘 = 3, ta có 𝑋 = (𝑥1, 𝑥2, 𝑥3)
𝑇 , 
𝐴(1) = [
4/11 −8/11 0
14/11 −23/22 0
9/11 −18/11 1
], 
𝐴(0) = [
0 0 −4/11
0 0 5/22
0 0 2/11
], 
𝐵(1) = [
−1 0 −1/11
0 −1 2/11
0 0 6/11
] , 𝐴 = [
1
1
1
], 
𝐶 = [
0 0 1
0 0 1
0 0 1
] , 𝐵 =
[
23
12
−
4
3
5
12
7
3
−
2
3
1
3
9
4
0
3
4 ]
, 
𝐺(𝑋) = 𝑦𝑞𝑘𝐴 + ℎ𝐵𝐹(𝑡𝑞𝑘+1, 𝑋) − 𝑋 
= 𝑦𝑘 [
1
1
1
] − [
𝑥1
𝑥2
𝑥3
] + ℎ𝐵 [
𝑓(𝑡𝑛+1, 𝑥1)
𝑓(𝑡𝑛+2, 𝑥2)
𝑓(𝑡𝑛+3, 𝑥3)
], 
𝐽𝐺(𝑋)
= ℎ𝐵 [
𝑓𝑦(𝑡𝑛+1, 𝑥1) 0 0
0 𝑓𝑦(𝑡𝑛+2, 𝑥2) 0
0 0 𝑓𝑦(𝑡𝑛+3, 𝑥3)
] 
− 𝐼3. 
Với 𝑘 = 5, các ma trận khối trong (4) và (4′) 
được cho như sau, các khối này do tác giả xây 
dựng, 
𝐴 = (1,1,1,1,1)𝑇, 
𝐵 =
[
1901
720
−
1387
360
109
30
269
90
−
133
45
49
15
237
80
−
99
40
39
10
−
637
360
251
720
−
73
45
29
90
−
69
40
27
80
134
45
−
116
45
68
15
425
144
−
175
72
25
6
−
56
45
14
45
−
25
72
95
144 ]
. 
3. Trình thực thi và thuật toán 
Trình thực thi của phương pháp khối liên tục 
dạng BDF bậc 𝑘 được trình bày dạng giải mã 
sau đây. Trong trình thực thi này, hằng số 
𝑡𝑜𝑙 > 0 được đưa vào làm tiêu chuẩn dừng 
cho phép lặp Newton, đảm bảo sai số phép 
lặp mỗi bước không vượt quá số này. Ngoài 
ra, hằng số M đưa ra giới hạn cho số lần lặp 
Đinh Văn Tiệp và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ ĐHTN 225(06): 424 - 431 
 Email: jst@tnu.edu.vn 429 
được phép mà nếu vượt quá số này, quá trình 
lặp thất bại. Số đoạn chia 𝑁 được yêu cầu là 
bội nguyên lần của 𝑘. 
INPUT hàm 𝑓(𝑡, 𝑦), [𝑎, 𝑏], giá trị ban đầu 
𝑦(𝑎) = 𝛼, số đoạn chia 𝑁 = 𝑘𝑝, độ tin cậy 
𝑡𝑜𝑙 > 0, số lần lặp tối đa M mỗi quá trình lặp. 
OUTPUT dãy xấp xỉ 𝑤𝑖 của nghiệm 
𝑦(𝑡𝑖), ∀𝑖 = 0,1, . . , 𝑁, 
hoặc thông báo quá trình thất bại. 
Trình thực thi 
Bước 1. Khởi tạo: 
𝑡0 ≔ 𝑎; 
𝑤0 ≔ 𝛼; 
ℎ ≔
𝑏 − 𝑎
𝑁
; 
𝑌 ≔ (0, ,𝑤0)
𝑇; % Viết 𝑌 = (𝑦1,  , 𝑦𝑘). 
ma trận 𝐴, 𝐵; 
𝑖 ≔ 0; 
hàm 𝐺(𝑋) = 𝑤0𝐴 + ℎ𝐵𝐹(𝑡𝑖+1, 𝑋) − 𝑋; 
𝐹𝐿𝐴𝐺 ≔ 0; % Điều kiện thoát vòng While 
Bước 2. Quy trình lặp: 
While (𝑖 < 𝑁) and (𝐹𝐿𝐴𝐺 = 0) do 
 𝑘 ≔ 1; 
 While 𝑘 ≤ 𝑀 do % Vòng While trong 
 tính 𝐽𝐺(𝑌), 𝐺(𝑌); 
 giải hệ phương trình 𝐽𝐺(𝑌)𝑋 = −𝐺(𝑌); 
 cập nhật 𝑌 ≔ 𝑌 + 𝑋; 
 If ‖𝑋‖ < 𝑡𝑜𝑙 then 
 cập nhật 𝑤0 ≔ 𝑦𝑘; 
 chấp nhận xấp xỉ 𝑌𝑖+1 ≔ 𝑌; 
 exit;%Thoát vòng While 
trong 
 Else If (𝑘 = 𝑀) then 
 OUT=”fail”; 
 𝐹𝐿𝐴𝐺 ≔ 1; 
 % Thông báo số vòng lặp vượt quá 
M. 
 % Thoát vòng lặp While ngoài. 
 End If. 
 𝑘 ≔ 𝑘 + 1; 
 End do. % Kết thúc vòng While trong. 
 𝑖 ≔ 𝑖 + 3; 
End do. % Kết thúc vòng While ngoài. 
Bước 3. Xuất kết quả 𝑌1, 𝑌𝑘+1,  , 𝑌𝑘(𝑝−1)+1, 
hoặc thông báo quá trình thất bại do vượt quá 
số vòng lặp M quy định. 
Chương trình sau đây lập trình cho phương 
pháp với 𝑘 = 3. 
Chương trình Matlab 
function 
outp=BDFblock(f,a,b,alpha,N,to
l,M) 
% NOTE: N must be a multiple 
of k. 
%----Exact solution to 
compare---- 
syms p(r); 
format long; 
eqn=diff(p,r)==f(r,p); 
icd=p(a)==alpha; 
as=dsolve(eqn,icd); 
xi=double(subs(as,r,a:(b-
a)/N:b)); 
ti=a:(b-a)/N:b; 
%--------------STEP 1---------
---- 
t0=a; 
y0=alpha; 
h=(b-a)/N; 
A=[1;1;1]; 
B=[23/12,-4/3,5/12;7/3,-
2/3,1/3;9/4,0,3/4]; 
syms s0 s u z s1 s2 s3; 
K=diag([subs(diff(f(s0,z),z),[
s0,z],[s+h,s1]),subs(diff(f(s0
,z),z),[s0,z],[s+2*h,s2]),subs
(diff(f(s0,z),z),[s0,z],[s+3*h
,s3])]); 
Jg=h*B*K-eye(3);% Jacobian of 
G 
Đinh Văn Tiệp và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ ĐHTN 225(06): 424 - 431 
430  Email: jst@tnu.edu.vn 
G=u*A+h*B*[f(s+h,s1);f(s+2*h,s
2);f(s+3*h,s3)]-
[s1;s2;s3];%function G 
t=[t0]; 
w=[y0]; 
Y=[0;0;y0]; 
%-------------STEP 2----------
--- 
i=0; 
FLAG=0; 
%Use to STOP outer while-loop 
while (i<N)&(FLAG==0)% Outer 
while 
 k=1; 
while k<=M % Inner while-loop 
J=double(subs(Jg,[s,s1,s2,s3],
[t0,Y(1),Y(2),Y(3)])); 
G1=double(subs(G,[u,s,s1,s2,s3
],[y0,t0,Y(1),Y(2),Y(3)])); 
 X=linsolve(J,-G1); 
 Y=Y+X; 
 if (norm(X)<tol) 
 w=[w,Y(1),Y(2),Y(3)]; 
 y0=Y(3); 
 =[t,t0+h,t0+2*h,t0+3*h]; 
 t0=t0+3*h; 
 break 
 elseif k==M 
 fprintf('Fail!’); 
 FLAG=1; 
 end 
 k=k+1; 
end 
i=i+3; 
end 
%---------------Step 3--------
---- 
outp=[t',w',xi',abs(w-xi)']; 
end 
%---------------END-----------
---- 
4. So sánh thực nghiệm 
Phương pháp thực nghiệm được đưa ra ở đây 
để so sánh tính chính xác của phương pháp 
BDF dạng khối liên tục (với 𝑘 = 3) và 
phương pháp BDF thông thường sử dụng trực 
tiếp vòng lặp Newton bậc 4 và bậc 5. Các kết 
quả được thể hiện trong bảng 1 dưới đây. 
Ví dụ. Xét các phương trình sau (xem [4]) 
𝑦′ = 𝑦 − 𝑡2 + 1, 𝑦(0) =
1
2
, 0 ≤ 𝑡 ≤ 2, (9) 
{
𝑦′ = 5𝑒5𝑡(𝑦 − 𝑡)2 + 1
𝑦(0) = −1, 𝑡 ∈ [0,1]
 (10) 
{
𝑦′ = −20𝑦 + 20 cos 𝑡 − sin 𝑡
𝑦(0) = 0, 𝑡 ∈ [0,2]
 (11) 
{
𝑦′ = −20(𝑦 − 𝑡2) + 2𝑡
𝑦(0) = 1/3, 𝑡 ∈ [0,1]
 (12) 
Các kết quả thực nghiệm được cho trong bảng 
1 để đánh giá sai số cho thấy ưu điểm của các 
phương pháp BDF dạng khối so với các 
phương pháp BDF thông thường ở việc tạo ra 
các giá trị xấp xỉ chính xác cao hơn với kích 
thước bước ℎ vừa phải. Độ chính xác của 
BDF dạng khối với 𝑘 = 3, ký hiệu 
BDFblock3, giảm dần khi ℎ tăng lên so với 
các phương pháp BDF4 và BDF5 (có bậc 4, 
và 5 tương ứng). Điều này cũng dễ hiểu từ 
việc các phương pháp này có bậc cao hơn so 
với phương pháp được mang ra đánh giá. Tuy 
nhiên, nếu so sánh với phương pháp cùng bậc, 
thì các phương pháp BDF dạng khối liên tục 
rõ ràng nổi trội hơn. Điều này càng được thể 
hiện rõ hơn trong các bài toán stiff với độ 
“nhất thời” ( tức “transient”) càng lớn. 
5. Kết luận 
Việc xây dựng trình thực thi cho các phương 
pháp BDF dạng khối theo phương pháp lặp 
Newton giúp khai thác tốt độ chính xác của 
các phương pháp này, đặc biệt trong các bài 
toán stiff. 
Cách tiếp cận thứ hai của trình thực thi dựa 
vào việc giải phương trình phi tuyến tính (5) 
liên tiếp qua mỗi đoạn k - bước được giải 
quyết tốt nhờ thuật toán Newton. 
Đinh Văn Tiệp và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ ĐHTN 225(06): 424 - 431 
 Email: jst@tnu.edu.vn 431 
TÀI LIỆU THAM KHẢO/ REFERENCES 
[1]. O. A. Akinfenwa, S. N. Jator, and N. M. Yao, 
“Continuous block backward differentiation 
formula for solving stiff ordinary differential 
equations,” Computers & Mathematics with 
Applications, vol. 65, no. 7, pp. 996-1005, 
2013. 
[2]. O. A. Akinfenwa, S. N. Jator, and N. M. Yao, 
“On The Stability of Continuous Block 
Backward Differentiation Formula For 
Solving Stiff Ordinary Differential 
Equations,” J. of Mod. Meth. in Numer. 
Math., vol. 3, no. 2, pp. 50-58, 2012. 
[3]. J. R. Cash, “Review paper. Efficient 
numerical methods for the solution of stiff 
initial-value problems and differential 
algebraic equations,” Proc. R. Soc. Lond. Ser. 
A Math. Phys. Eng. Sci., vol. 459, no. 2032, 
pp. 797-815, 2003. 
[4]. R. L. Burden, and J. D.Faires, Numerical 
Analysis, 9th Edition, Brooks/Cole, 2010. 
Bảng 1. Các giá trị sai số tuyệt đối ở bước cuối cùng tạo ra bởi các phương pháp BDF4, BDF5 thông 
thường và phương pháp BDF3 dạng khối liên tục (ký hiệu BDFblock3). 
Equations BDF4 BDF5 BDFblock3 
(9), M=10, 
N=6,12, 30, 
tol=0.001, 
5.3× 10−3, 6.5× 10−4, 
2.4× 10−5 
1.4× 10−3, 
1.1× 10−4, 
1.8× 10−6 
6.13× 10−2, 
5.64× 10−3, 
3.05× 10−4 
(10), M=10, 
N=6, 12, 30, 
tol=0.001 
3.95× 10−3, 
8.3× 10−5, 
3.2× 10−6 
1.8× 10−3, 
4.4× 10−5, 
3.5× 10−6 
3.1× 10−4, 
2.5× 10−5, 
6.5 × 10−6 
(11), M=10, 
N=6,12, 30, 300 
tol=0.001 
7.8× 103, 
7.3× 10−2, 
9.95 × 10−8, 
1.84× 10−11 
1.84× 106, 
1.28, 
1.1 × 10−5, 
4.6× 10−14 
5.5× 10−4, 
5.7× 10−6, 
2.4 × 10−7, 
5.6× 10−10 
(12), M=10, 
N= 6, 12, 30, tol=0.001 
0.42, 
8.74 × 10−5, 
1.02× 10−9 
1.13, 
2.2 × 10−4, 
2.53× 10−7 
1.48× 10−4, 
3.79× 10−8, 
2.62× 10−10