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,
Dưới đây trình bày cách xây dựng các công thức Euler, Runge-Kutta, xuất phát từ
qui tắc cầu phương cơ bản (xem, thí dụ, [2]).
1.2.1. Quy tắc cầu phương cơ bản và giải số phương trình vi phân
Quy tắc cầu phương cơ bản (basic quadrature rules) có thể được coi là phương pháp
quan trọng để tính tích phân. Vì giải phương trình vi phân thường (1.1) với điều kiện
ban đầu (1.2) tương đương với việc giải phương trình tích phân
0
0
( ) ( ( ), )
t
t
x t x f x s s ds= +
∫
(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
[ ]
,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
ϕ
=
=
∑
,
trong đó
1,
( )
( )
( )
s
i
j
i i j
j i
t c
L t
c c
= ≠
−
=
−
∏
. Khi ấy
1
( ) ( )
b
s
j j
j
a
f t dt f c
ω
=
≈
∑
∫
.
5
Các trọng số
j
ω
được tính theo công thức
( ) .
b
j j
a
L t dt
ω
=
∫
Nếu
1s =
thì đa thức nội suy
1
( ) ( )t f c
ϕ
≡
và ta có:
1
( ) ( ) ( ).
b
a
f t dt b a f c≈ −
∫
Ta nói độ chính xác (precision) của quy tắc cầu phương là
p
nếu quy tắc này chính
xác cho mọi đa thức bậc nhỏ hơn
p
, tức là với mọi đa thức
( )
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( ).
p
h
+
Ta xét một số trường hợp đặc biệt.
• Nếu chọn
1s
=
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):
( ) ( ) ( ).
b
a
f t dt b a f a≈ −
∫
(1.5)
Nếu
( )x t
là nghiệm của phương trình vi phân (1.1) - (1.2) (nghiệm của phương trình
tích phân (1.4)) thì:
( ) ( ) ( ( ), )
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 thì
nghiệm được xây dựng về bên trái của
0
t
). Dưới đây ta coi
0h
>
, trường hợp
0h
<
có thể được xét tương tự.
Từ công thức (1.7) ta có
6
O
f
x
C
D
b
B
a
A
1 0 0 0
. ( , )x x h f x t= +
;
2 1 1 1
. ( , )x x h f x t= +
;
;
1
. ( , )
n n n n
x x h f x t
+
= +
.
Đây chính là công thức Euler tiến quen thuộc. Hình 1.1
• Nếu chọn
1s =
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
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
1s
=
và
1
2
a b
c
+
=
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 ABMN (Hình 1.3):
( ( ), ) . ( ( ), )
2 2
t h
t
h h
f x s s ds h f x t t
+
= + +
∫
Từ đây ta có: Hình 1.3
( ) ( ) . ( ( ), )
2 2
h h
x t h x t h f x t t+ − = + +
.
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).
7
E
b
B
a
A
F
O
f
x
2
ba
+
MN
b
B
a
A
O
f
x
• Nếu chọn
2s =
và
1 2
,c a c b= =
thì
1
( )
( )
t b
L t
a b
−
=
−
và
2
( ) .
( )
t a
L t
b a
−
=
−
Suy ra
2
1 1
1 ( )
( )
( ) ( ) 2 2
b
b b
a a
a
t b t b b a
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
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
thang ABED, Hình 1.4)
thì ta được:
( ( ), ) [ ( ( ), ) ( ( ), )]
2
t h
t
h
f x s s ds f x t h t h f x t t
+
≈ + + +
∫
.
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
2p =
.
• Nếu chọn
3s
=
và
1 2 3
, ,
2
a b
c a c c b
+
= = =
thì, đặt
h b a
= −
, ta có:
1
2
( )( )
2
2
( ) ( )( ),
2
( )( )
2
a b
t t b
a b
L t t t b
a b
h
a a b
+
− −
+
= = − −
+
− −
8
E
D
b
B
a
A
O
f
x
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 b
L t t a t
a b
h
b a b
+
− −
+
= = − −
+
− −
Suy ra
1 1
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
ω
+ −
= = − − = − − −
− − − −
= − − − = −
−
= =
∫ ∫ ∫
∫
và
2 2
2 2
3 2
2
4 4
( ) ( )( ) ( )( ( ))
4 ( ) ( ) 4
[ ( )] .
3 2 6
b b b
a a a
b
a
L t dt t a t b dt t a t a a b dt
h h
t a t a h
a b
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
( ) ( ) [ ( ( ), ) 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
+ + +
− = + + + +
.
9
Đâ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
1.2.2.1. Dẫn tới phương pháp Runge - Kutta
Vì phương pháp ẩn đòi hỏi tại mỗi bước phải giải một phương trình phi tuyến, điều
này không đơn giản, nên ta cố gắng xây dựng các công thức Runge-Kutta hiển từ
công thức hình thang ẩn, công thức điểm giữa ẩn và công thức Runge-Kutta kinh điển
cấp bốn ẩn tương ứng như sau.
• Trong công thức hình thang ẩn:
1 1 1
[ ( , ) ( , )]
2
n n n n n n
h
x x f x t f x t
+ + +
= + +
,
ta thay giá trị
1n
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
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
1
( ( , )
2 2
n n n n
h h
x x hf x t t
+
= + + +
ta nhận được phương pháp trung điểm hiển (explicit midpoint method):
1 1
2
ˆ
( , )
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
10
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
+
= + + + + =
trong đó:
1
1
2
2
3
4 3
( , );
( , );
2 2
( , );
2 2
( , ).
n n
n n
n n
n n
k f x t
hk
h
k f x t
hk h
k f x t
k f x hk t
=
= + +
= + +
= +
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
,
( )
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
+
=
= +
∑
, (1.8)
trong đó
i
X
là vectơ
n
-chiều và là nghiệm của hệ phương trình phi tuyến
++=
∑
=
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ố
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 ≠
với
j i≥
nào
đó thì ta có phương pháp Runge-Kutta ẩn. Khi ấy tại mỗi bước ta phải giải một hệ
ns
11
phương trình phi tuyến (tuyến tính nếu
( , ) ( ) ( )f x t B t x g t= +
) để tìm
s
vectơ
i
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
A
T
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)
trong đó
1 2 2 21 1
( , ); ( , ).
n n n n
k f x y k f x c h y ha k= = + +
Khái niệm
s
-nấc (
s
-stage) thể hiện rằng số lần tính các giá trị của hàm
f
(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+ = = =
.
Đâ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
0b ≠
tự do. Khi ấy các hệ số còn lại biểu diễn qua
2
0b ≠
bởi các công thức:
1 2
1b b= −
,
2
2
1
2
c
b
=
,
21
2
1
2
a
b
=
.
12
Chọn
2
1
2
b =
, thì
1
1
2
b =
và
21 2
1a c= =
. Khi ấy ta có một phương pháp Runge-Kutta
cấp hai cho phép tính
1n
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
1b =
thì
21
1
2
a =
,
1
0b =
và
2
1
2
c =
. Khi ấy ta có công thức
1 1 1 2 2
( ) . . ( ( , ), )
2 2
n n n n n n n
h h
x x h b k b k x h f x f x t t
+
= + + = + + +
.
Phương pháp tính theo công thức trên được gọi là phương pháp Euler-Cauchy.
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
ρ σ σ
+ − + − + −
= = =
=
∑ ∑ ∑
. (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
0
σ
=
thì phương pháp là phương pháp hiển. Nếu
0
0
ρ
≠
và
0
0
σ
≠
thì ta có phương pháp ẩn.
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ử)
13
( )
[ ]
0
, 0 , 0,1x 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
λλλ
+=
, 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
σ
=
+
1
,
, 2,1,0
=
n
,
trong đó
0
x
cho trước và
σ
nói chung là một số phức. Nghiệm của phương trình
này là
0
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
λ
σ
=
.
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
λ
, miền ổn định của nghiệm chính xác phải là nửa mặt phẳng mở bên
trái.
Phương pháp một bước được gọi là ổn định tuyệt đối nếu
1
≤
σ
và ổn định tương
đối nếu
h
e
λ
σ
≤
.
Nếu
λ
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.
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.
14
Không có nhận xét nào:
Đăng nhận xét