ĐẠI HỌC THÁI NGUYÊN
TRƯỜNG ĐẠI HỌC KHOA HỌC

Ngô Thị Thu Hương

MỘT SỐ PHƯƠNG PHÁP GIẢI SỐ
PHƯƠNG TRÌNH VÀ HỆ PHƯƠNG TRÌNH
VI PHÂN CẤP CAO

LUẬN VĂN THẠC SĨ TOÁN HỌC

Thái Nguyên - 2017


ĐẠI HỌC THÁI NGUYÊN
TRƯỜNG ĐẠI HỌC KHOA HỌC

Ngô Thị Thu Hương

MỘT SỐ PHƯƠNG PHÁP GIẢI SỐ
PHƯƠNG TRÌNH VÀ HỆ PHƯƠNG TRÌNH
VI PHÂN CẤP CAO

Chuyên ngành: Toán ứng dụng
Mã số: 60.46.01.12

LUẬN VĂN THẠC SĨ TOÁN HỌC

NGƯỜI HƯỚNG DẪN KHOA HỌC

TS. Vũ Vinh Quang

1.2.2 Hàm lưới . . . . . . . . . . . . . . . . .
1.2.3 Công thức Taylor . . . . . . . . . . . .
1.2.4 Một số công thức xấp xỉ đạo hàm . . .
1.3 Thuật toán truy đuổi 3 đường chéo . . . . . . .
1.4 Phương pháp lưới giải bài toán biên cho phương

. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
trình

. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
cấp 2


7
7
7
8
10
12

2 Phương pháp số giải phương trình vi phân phi tuyến cấp cao
và hệ phương trình vi phân với hệ điều kiện đầu
16
2.1 Cơ sở lý thuyết về phương pháp Runge-Kutta . . . . . . . . . . 16
2.1.1 Phương pháp Euler 1 . . . . . . . . . . . . . . . . . . . . 17
2.1.2 Phương pháp Euler 2 . . . . . . . . . . . . . . . . . . . . 18
2.1.3 Thuật toán RK4 . . . . . . . . . . . . . . . . . . . . . . 19
2.2 Phương pháp Runge-Kutta đối với hệ phương trình vi phân phi
tuyến . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3 Phương pháp Runge-Kutta đối với phương trình vi phân cấp cao 21
2.4 Giới thiệu thư viện QH_2015 . . . . . . . . . . . . . . . . . . . 23


ii

3

Phương pháp lặp giải mô hình các bài
cấp 4
3.1 Giới thiệu . . . . . . . . . . . . . . . . .
3.2 Nghiên cứu các tính chất của nghiệm . .
3.3 Phương pháp xây dựng sơ đồ lặp . . . .
3.3.1 Cơ sở lý thuyết . . . . . . . . . .

Tôi xin chân thành cảm ơn Ban giám hiệu, phòng Đào tạo trường Đại học
Khoa học - Đại học Thái Nguyên đã quan tâm và giúp đỡ tôi trong suốt thời
gian học tập tại trường. Tôi xin cảm ơn quý thầy cô Khoa Toán - Tin và đặc
biệt là PGS.TS. Nguyễn Thị Thu Thủy, trưởng Khoa Toán - Tin, đã luôn
quan tâm, động viên, trao đổi và đóng góp những ý kiến quý báu trong suốt
quá trình học tập, nghiên cứu và hoàn thành luận văn.
Cuối cùng, tôi muốn bày tỏ lòng biết ơn sâu sắc tới những người thân trong
gia đình, đặc biệt là bố mẹ. Những người luôn động viên, chia sẻ mọi khó
khăn cùng tôi trong suốt thời gian qua và đặc biệt là trong thời gian tôi theo
học khóa thạc sỹ tại trường Đại học Khoa học - Đại học Thái Nguyên.
Thái Nguyên, tháng 10 năm 2017
Tác giả luận văn

Ngô Thị Thu Hương


1

Bảng ký hiệu
R
R+
R ∪ {±∞}
Rn
Ωh
C 1 [0; L]
A∪B
A∩B
x, y
[x, y]
l2

3.4
3.5

Trường
Trường
Trường
Trường
Trường

hợp
hợp
hợp
hợp
hợp

biết trước nghiệm đúng (Tofuma_moi.m) . . . . . .
biết trước nghiệm đúng (Tofuma_moi.m) . . . . . .
không biết trước nghiệm đúng (Tofuma_moi_xx.m)
biết trước nghiệm đúng (Tofuma_tq.m) . . . . . . .
không biết trước nghiệm đúng (Tofuma_tp_xx.m)

37
38
39
40
42


3



giữa
giữa
giữa
giữa
giữa

nghiệm
nghiệm
nghiệm
nghiệm
nghiệm

đúng
đúng
đúng
đúng
đúng

và
và
và
và
và

nghiệm
nghiệm
nghiệm
nghiệm
nghiệm


.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

37
38
39
41
43


5

Chương 2: Phương pháp số giải phương trình vi phân phi tuyến cấp cao và
hệ phương trình vi phân với hệ điều kiện đầu.
Chương 3: Phương pháp lặp giải mô hình các bài toán biên phi tuyến
cấp 4.


6

Chương 1

Một số kiến thức cơ bản
Trong chương này chúng tôi trình bày một số kết quả lý thuyết về các sơ
đồ lặp, phương pháp sai phân đối với phương trình vi phân cấp 2 và thuật
toán truy đuổi 3 đường chéo. Những kết quả này là những kiến thức bổ trợ
cho việc trình bày các kết quả chính trong chương 2 và chương 3. Các kết quả
lý thuyết được tham khảo trong các tài liệu [1, 2, 3].
1.1
1.1.1

Một số khái niệm cơ bản của Giải tích hàm
Không gian metric

Định nghĩa 1.1.1 Tập X của các phần tử x, y, z . . . được gọi là không gian
metric nếu như với mọi phần x, y bất kì đều tương ứng với với 1 số không âm
d(x,y) thỏa mãn các điều kiện sau:
+ d(x,y)>0, d(x,y)=0 khi và chỉ khi x=y
+ d(x,y)=d(y,x)

1.2
1.2.1

(n ≥ 1),
(n ≥ 1).

Phương pháp sai phân
Lưới sai phân

Ta chia đoạn [a, b] thành N đoạn con bằng nhau, mỗi đoạn con dài
b−a
h=
bởi các điểm xi = a + ih, i = 0, 1, ..., N . Khi đó tập các điểm xi gọi
N
là một lưới sai phân trên [a, b] ký hiệu là Ωh , mỗi điểm xi gọi là một nút của
lưới, h gọi là bước đi của lưới.
1.2.2

Hàm lưới

Xét hàm số u(x) xác định trên đoạn [a,b], khi đó tập giá trị của hàm trên
các điểm lưới được gọi là hàm lưới ui = u (xi ).
1.2.3

Công thức Taylor

Giả sử u(x) là một hàm số xác định và có đạo hàm đến cấp m + 1 trong
một khoảng (α, β) chứa x và x + ∆x, ∆x có thể âm hay dương. Khi đó ta có
công thức khai triển:
(∆x)2

(∆x)2
(∆x)m m
u(x+∆x) = u(x)+∆xu (x)+
u (x)+...+
u (x)+O (∆x)(m+1) .
2!
m!
1.2.4

Một số công thức xấp xỉ đạo hàm

Xét không gian lưới Ωh với bước lưới h =

b−a
, chúng ta xét công thức
N

khai triển Taylo tổng quát đối với hàm u(x):
h2
h3 (3) h(4) (4)
u(x ± h) = u(x) ± hu + u ± u +
u + .. + O(hn ).
2
6
24
Xuất phát từ công thức Taylo dừng với 3 số hạng khai triển, chúng ta nhận
được các công thức xấp xỉ đạo hàm với độ chính xác cấp 1 và cấp 2 như sau:
u (xi ) =
u (xi ) =


h5 (5)
h6 (4)
ui − ui + ui −
ui +
ui + O(h6 ).
2
6
24
120
720

Từ đây suy ra
ui+1 + ui−1

h4 (4)
h6 (6)
= 2ui + h ui + ui +
ui + O(h6 ).
12
360
2


9

hay
ui =

ui−1 − 2ui + ui+1 h2 (4)
h4 (6)


(xi ) =

f (xk )Lk (xi ) + R(xi ); i = 0, 1, ..., n.
k=0

Trong đó độ chính xác của khai triển càng cao khi số mốc nội suy càng lớn.
Như vậy các đạo hàm của hàm tại điểm xi được tính qua các giá trị hàm tại
các điểm xk với độ chính xác bậc cao.
Bằng tính toán cụ thể, ta có các công thức sau đây
Đạo hàm cấp 1
1
(−25f0 + 48f1 − 36f2 + 16f3 − 3f4 ) + O(h4 ),
12h
1
f (xn ) =
(25fn − 48fn−1 + 36fn−2 − 16fn−3 + 3fn−4 ) + O(h4 ).
12h

f (x0 ) =

Đạo hàm cấp 2
f (x0 ) =

1
(35f0 − 104f1 + 114f2 − 56f3 + 11f4 ) + O(h2 ),
2
12h

f (x1 ) =

1
(f0 − 4f1 + 6f2 − 4f3 + f4 ) + O(h),
h4

f (4) (x1 ) =

1
(f0 − 4f1 + 6f2 − 4f3 + f4 ) + O(h),
h4

f (4) (xk ) =

1
(fk−2 − 4fk−1 + 6fk − 4fk+1 + fk+2 ) + O(h),
h4

f (4) (xn−1 ) =
f (4) (xn ) =

1
(fn − 4fn−1 + 6fn−2 − 4fn−3 + fn−4 ) + O(h),
h4

1
(fn − 4fn−1 + 6fn−2 − 4fn−3 + fn−4 ) + O(h).
h4

Các công thức tính xấp xỉ đạo hàm trên sẽ được sử dụng để xây dựng các
lược đồ sai phân giải các phương trình vi phân.
1.3

 ............................................................
0
0
0
...
Cn−1
Bn−1

0
0
0
...
An
Cn
dạng 3 đường chéo.
X = (x0 , x1 , ..., xn ) : ẩn phải tìm;
D = (f0 , f1 , ..., fn ) : cột vế phải.










(1.1)



; i = 1, 2, ..., n − 1.
(Ci − αi Ai )

Để có các hệ số đó, ta cần xác định α1 , β1 ;
Do phương trình đầu tiên: −C0 x0 + B0 x1 = −d0 .
B0
d0
hay: x0 =
x1 +
= α1 x1 + β1 có dạng (1.2).
C0
C0
B0
d0
Nghĩa là: α1 =
, β1 =
.
C0
C0
Nhưng để tính theo công thức (1.2), ta cần có xn . Ẩn xn sẽ được tìm nhờ


12

phương trình cuối cùng của hệ (1.1)
An xn+1 − Cn xn = −dn .
Lại theo (1.2) ta có: xn−1 = αn xn + βn .
βn An + dn
.
C n − An αn

1
Co
Co


xi = αi+1 xi+1 + βi+1 ; i = 1, 2, .., n − 1,




β A + dn

 xn = n n
.
Cn − An αn

Loại trừ xn−1 từ hệ này, ta suy ra được: xn =

(1.3)

Xuất phát từ các công thức trên, thuật toán truy đuổi 3 đường chéo được
thực hiện bằng thuật toán sau đây:
Thuật toán
B1
f1
Bước 1: Xuất phát từ α1 = − , β1 =
, xác định tất cả các giá trị
C1
C1
αi , βi , ∀i = 2, 3, ..., n theo công thức:


13

Thuật toán thông thường
Sử dụng các công thức sai phân với độ chính xác cấp 2, ta có
ui+1 − 2ui + ui−1
+ O(h2 ),
h2
u1 − u0
un−1 − un
u0 =
+ O(h), un =
+ O(h).
2h
−h
u =

Khi đó (1.4) được đưa về hệ phương trình sai phân sau đây


= 2hA,
 (2hα0 + α1 )u0 − α1 u1
ui+1 − 2ui + ui−1
= h2 fi ,

 −β u
1 n−1 + (2hβ0 + β1 )un = 2hB.

(1.5)


2
6
24
120
720

Hay
ui =

h4 (4)
ui−1 − 2ui + ui+1 h2
−
f
−
fi + O(h4 ).
i
2
h
12
360

Vậy ta có lược đồ sai phân với độ chính xác cấp 6
ui−1 − 2ui + ui+1
h2
h4 (4)
= fi + fi +
f ; i = 1, 2, ..., n − 1
h2
12
360 i





A=









a00 a01
a10 a11
0 a21
0 0
0 0
...

a02 a03 a04
a12 0 0
a22 a23 0
a31 a32 a33
0 0 0
...

0
0













U = (u0 , u1 , ..., un )T ; F = (F0 , F1 , ..., Fn )T
Hệ phương trình sai phân (1.3) chính là hệ phương trình sai phân tương
ứng với bài toán biên cho phương trình vi phân (1.1) với độ chính xác cấp 4.
Trong trường hợp khi điều kiện biên có dạng Dirichlet u(a) = A; u(b) = B thì
ta thu được lược đồ sai phân với đội chính xác cấp 6 (do không phải sai phân
điều kiện biên).
Bởi vì ma trận A của hệ không phải dạng 3 đường chéo, do đó hệ không
giải được bằng thuật toán truy đuổi.
Nhận xét: Có thể thấy rằng qua một số hữu hạn phép biến đổi thì ta có
thể chuyên hệ ban đầu về dạng 3 đường chéo với các hệ số được xác định bởi
các công thức:
α1
α1
a00 = α0 + ; a01 = − ;
h
h
ai−1i = 1; aii = −2; aii+1 = 1; i = 1, 2, ..., n − 1,
β1

2.1

Cơ sở lý thuyết về phương pháp Runge-Kutta

Xét bài toán Cauchy, hay còn gọi là bài toán giá trị ban đầu: Tìm y(x)
thỏa mãn điều kiện:
y = f (x, y), x0 ≤ x ≤ x
y(x0 ) = y0 , (y, f ∈ Rm ).

(2.1)

Ta quan tâm đến phương pháp tìm nghiệm bằng số của (2.1) tại điểm
x1 = x0 + h với bước h > 0. Hai nhà toán học người Đức là Runge và Kutta
đề xuất một phương pháp tìm nghiệm y1 , chỉ phải tính một hàm f (x, y) tại
một số điểm khác nhau:
Đặt y1 = y0 + ∆y0 , trong đó ∆y0 = pr1 k1 (h) + ... + prr kr (h).
ki (h) = hf (ξi , vi ) ; ξi = x0 + αi h; α1 = 0;
i = 1, r
vi = y0 + βi1 k1 (h) + ... + βi,i−1 ki−1 (h) ,
Gọi ϕr (h) := y (x0 + h) − y1 = y (x0 + h) − y (x0 ) − ∆y0 là sai số (địa


17
(s+1)

phương) của phương pháp Runge-Kutta. Nếu ϕr
r

(0) = 0 thì:


(2.2)

Xuất phát từ phương pháp Runge-Kutta tổng quát, chúng ta thu được một
số phương pháp cụ thể như sau:
2.1.1

Phương pháp Euler 1

Xét trường hợp riêng của phương pháp Runge-Kutta khi r = 1.
∆y0 = p11 k1 (h) ; k1 (h) = hf (x0 , y0 ) .
Ta có:
ϕ1 (h) := y (x0 + h) − y0 − p11 hf (x0 , y0 ) ; ϕ1 (0) = 0,
ϕ1 (0) = y0 − p11 f (x0 , y0 ) = (1 − p11 ) f (x0 , y0 ) .
Để ϕ1 (0) = 0 với mọi hàm f , ta phải có p11 = 1. Nói chung ϕ1 (0) = y0 = 0.
Vậy ∆y0 = p11 k1 (h) = hf (x0 , y0 ). Ta nhận được công thức Euler:
y1 = y0 + hf (x0 , y0 ) .

(2.3)

Nói chung chúng ta có công thức tổng quát
yk+1 = yk + hf (xk , yk ) , xk = x0 + kh.
Sai số địa phương của phương pháp:
R1 (h) =

ϕ1 (ξ) 2 y (ξ) 2
h =
h = O h2 .
2!
2


(l)

y0 = p21 k1 (0) + p22 k2 (0); l = 1, 2 .

(2.4)

Vì k1 (h) = hf (x0 , y0 ) nên k1 (0) = 0; k1 (0) = f (x0 , y0 ) và k1 (0) = 0. Tiếp
theo k2 (h) = hf (ξ2 , v2 ) trong đó ξ2 = x0 + α2 h; v2 = y0 + β21 k1 (h). Ta có:
k2 (h) = f (ξ2 , v2 ) + h

∂f
∂f
α2 +
β21 k1 (h) ≡ f (ξ2 , v2 ) + h [.] .
∂x
∂y

Dễ thấy k2 (0) = 0; k2 (0) = f (x0 , y0 ).
∂f0
∂f0
+ β21 f0
.
∂x
∂y
∂f0 ∂f0
∂f ∂f
;
là đạo hàm
;
Ở đây chúng ta dùng kí hiệu f0 := f (x0 , y0 ) ;

∂y

(2.6)

Vì công thức Runge-Kutta (ứng với r = 2) đúng cho mọi hàm f nên để
(2.6) nghiệm đúng, cần 1 − 2α2 p22 = 1 − 2p22 β21 = 0.
1
1
Như vậy: α2 = β21 = 1; p21 = p22 = và ∆y0 = h {f (x0 , y0 ) + f (x0 + h, y0 + hf (x
2
2
Ta nhận được công thức gọi là công thức Euler 2 sau:

 y0 := y0 + hf (x0 , y0 ) ,
(2.7)
1
 y1 = y0 + h [f (x0 , y0 ) + f (x1 , y0 )] .
2
Có thể thấy rằng phương pháp Euler 2 có độ chính xác O(h2 ).
2.1.3

Thuật toán RK4
4

Khi r = 4, ϕ4 (h) = y (x0 + h) − y (x0 ) −

p4i ki (h).
i=1

4

h
k1
k2 = hf (x0 + , y0 + ),
2
2
h
k2
k3 = hf (x0 + , y0 + ),
2
2
k4 = hf (x0 + h, y0 + k3 ),

(2.8)


20
(5)
5 ϕ4 (ξ)

Có thể chứng minh rằng R4 (h) = h
2.2

là sai số địa phương của RK4.
120
Phương pháp Runge-Kutta đối với hệ phương trình vi phân
phi tuyến

Xét hệ phương trình vi phân cấp 1

dx1


(2.9)

Nghiệm của hệ cần thỏa mãn hệ điều kiện đầu:
x1 (0) = x10 , ..., xm (0) = xm0 .
Đặt
U = (x1 , x2 , ..., xm )T ;
F (t, U ) = (f1 (t, U ), ..., fm (t, U ))T ,
U0 = (x10 , ..., xm0 )T .
Khi đó hệ (2.9) tương với phương trình vi phân cấp 1 dạng vecto
dU
= F (t, U ); U (t = 0) = U0 .
dt

(2.10)

Sử dụng kết quả của phương pháp Runge Kutta với độ chính xác bậc 4,
chúng ta có thể xây dựng lược đồ giải số hệ phương trình vi phân cấp bằng
sơ đồ dạng vecto như sau:
Sơ đồ QH_m
1
Uk+1 = Uk − [K1 + 2K2 + 2K3 + K4 ],
6
h
K1
K1 = hF (tk , Uk ), K2 = hF (tk + , Uk +
),
2
2
h