Bài giảng Lý thuyết về số gần đúng - Chương 6: Giải gần đúng phương trình vi phân

Chương 6
GIẢI GẦN ĐÚNG
PHƯƠNG TRÌNH VI PHÂN
I. GIẢI GẦN ĐÚNG PTVP CẤP 1 : 
Xét bài toán Cauchy : tìm nghiệm y=y(x) của 
phương trình vi phân với giá trị ban đầu y0
y’ = f(x, y), ∀x ∈ [a,b]
y(a) = y0
Các phương pháp giải gần đúng :
➢ Công thức Euler
➢ Công thức Euler cải tiến
➢ Công thức Runge-Kutta
1. Công thức Euler : 
Để tìm nghiệm gần đúng của bài toán Cauchy ta 
chia đoạn [a,b] thành n đoạn nhỏ bằng nhau với 
bước h = (b-a)/n
xo= a, x1 = x0 +h, ... , xk = x0 + kh, ... , xn = b
Nghiệm gần đúng của bài toán là dãy {yk} gồm các 
giá trị gần đúng của hàm tại xk 
Ta có yk ≈ y(xk) , k =0, n
Công thức Euler : 
yk+1 = yk + h f(xk, yk) , k = 0, n-1
Ví dụ : Dùng công thức Euler tìm nghiệm gần 
đúng của bài toán Cauchy 
 y’ = y – x2 +1, 0≤x≤1
y(0) = 0.5
với n = 5 
Tính sai số biết nghiệm chính xác là : 
y(x) = (x+1)2 – 0.5ex
giải
ta có h = 0.2
x0 = 0, x1 = 0.2, x2 = 0.4, x3 = 0.6, x4 = 0.8, x5 = 1
Công thức Euler
 y0 = 0.5
 yk+1 = yk + 0.2 (yk - xk
2 +1) 
k xk yk y(xk) |y(xk) - yk |
0 0 0.5 0.5 0
1 0.2 0.8 0.8292986 0.0292986
2 0.4 1.152 1.2140877 0.0620877
3 0.6 1.5504 1.6489406 0.0985406
4 0.8 1.98848 2.1272295 0.1387495
5 1 2.458176 2.6408591 0.1826831
* Nhận xét : công thức Euler đơn gian, nhưng sai số 
còn lớn nên ít được sử dụng 
2. Công thức Euler cải tiến : 
yk+1 = yk + (k1+k2)/2 k = 0,1, ..., n-1
k1 = hf(xk, yk), 
k2 = hf(xk+h, yk + k1)
Ví dụ : 
Làm lại ví dụ trước nhưng dùng công thức Euler cải 
tiến
giải
ta có h = 0.2
x0 = 0, x1 = 0.2, x2 = 0.4, x3 = 0.6, x4 = 0.8, x5 = 1
Công thức Euler cải tiến
yo = 0.5
yk+1 = yk + (k1 +k2) /2
k1= 0.2(yk - xk
2 +1) 
k2 = 0.2(yk + k1 – (xk+0.2)
2 +1)
k xk yk y(xk) |y(xk) - yk |
0 0 0.5 0.5 0
1 0.2 0.826 0.8292986 0.0033
2 0.4 1.20692 1.2140877 0.0072
3 0.6 1.6372424 1.6489406 0.0117
4 0.8 2.1102357 2.1272295 0.0170
5 1 2.6176876 2.6408591 0.0232
3. Công thức Runge Kutta bậc 4 : 
* Chú ý : Lập công thức Runge-Kutta bằng 
máy tính casio không được vì công thức quá 
dài, không đủ bộ nhớ, ta phải tính trực tiếp 
Ví dụ : Xét bài toán Cauchy 
 y’ = 2.7xy + cos (x+2.7y), 1.2≤x
y(1.2) = 5.4
Dùng công thức Runge-Kutta tính gần đúng y(1.5) 
với bước h = 0.3
xo = 1.2, yo = 5.4, y1=y(1.5)
y1 = y0 + (K1+ 2K2+ 2K3+ K4) /6
Công thức Runge-Kutta bậc 4
giải
Bấm máy (lập hàm dùng phím calc) ta được 
K1 = 4.949578057 K2 = 8.367054617
K3 = 10.33000627 K4 = 19.41193853
y(1.5) = 15.69260639 ≈ 15.6926
II. GIẢI GẦN ĐÚNG HỆ PTVP : 
Xét hệ phương trình vi phân cấp 1
y’1 = f1(x, y1, y2, ..., ym)
y’2 = f2(x, y1, y2, ..., ym)
. . .
y’m = fm(x, y1, y2, ..., ym)
với a≤ x ≤ b và thỏa điều kiện ban đầu
y1(a) = α1, y2(a) = α2, .... , ym(a) = αm
Nghiệm y = (y1, y2, , ym)
Để tìm nghiệm gần đúng, ta chia đoạn [a,b] thành n 
đoạn nhỏ bằng nhau với bước h = (b-a)/n và các 
điểm chia
xo= a, x1 = x0 +h, ... , xk = x0 + kh, ... , xn = b
Công thức Euler : 
yi k+1 = yi k + h fi(xk, y1 k,  , ym k)
∀i=1..m; k = 0.. n-1
Nghiệm gần đúng là dãy { yk=(y1 k, y2 k, , ym k)}
với yi k ≈ yi(xk), i=1,m
Công thức Euler cải tiến : 
yi k+1 = yi k + (K1 i + K2 i) / 2
K1 i = h fi(xk, y1 k,  , ym k)
K2 i = h fi(xk+h, y1 k+K1 1,  , ym k+K1 m)
∀i=1,m; k = 0, n-1
Công thức Runge-Kutta bậc 4 : 
yi k+1 = yi k + (K1 i+2K2 i+2K3 i+K4 i) / 6
K1 i = h fi(xk, y1 k,  , ym k)
K2 i = h fi(xk+h/2, y1 k+K11/2,  , ym k+K1 m/2)
K3 i = h fi(xk+h/2, y1 k+K21/2,  , ym k+K2 m/2)
K4 i = h fi(xk+h, y1 k+K31,  , ym k+K3 m)
∀i=1,m; k = 0, n-1
Ví dụ : Sử dụng công thức Euler giải gần đúng hệ 
pt vi phân
y’1 = 3y1 + 2y2 – (2x
2 +1)e2x
y’2 = 4y1 + y2 + (x
2 +2x –4) e2x
với 0 ≤x≤0.5
điều kiện ban đầu y1(0)=y2(0)=1
bước h = 0.1
Công thức Euler
y1 0 = 1
y1 k+1 = y1 k + h (3y1k + 2y2 k – (2xk
2 +1)e2xk)
y2 0 = 1
y2 k+1 = y2 k + h (4y1k + y2 k + (xk
2 +2xk –4) e
2xk)
xk y1k y2k
0 1 1
0.1 1.4 1.1
0.2 1.9154 1.3071
0.3 2.5903 1.6729
0.4 3.4870 2.2732
0.5 4.6940 3.2187
III. GIẢI GẦN ĐÚNG PTVP CẤP CAO: 
Xét phương trình vi phân cấp m 
y(m) = f(x, y, y’, ... , y(m-1)), a≤x≤b
với điều kiện ban đầu 
y(a) = α1, y’(a) = α2, .... , y
(m-1)(a) = αm
Đặt y1 = y, y2 = y’, y3 = y”, ... , ym = y
(m-1)
Ta chuyển phương trình vi phân bậc m về hệ m 
phương trình vi phân cấp 1 
với điều kiện ban đầu 
y1(a) = α1, y2(a) = α2, .... , ym(a) = αm,
y’1 = y2
y’2 = y3
. . .
y’m-1 = ym
y’m = f(x, y1, y2, ... , ym)
Ví dụ : Sử dụng công thức Euler giải gần đúng pt 
vi phân cấp 2 (tính xấp xỉ y và y’)
y “ – 2 y’ + 2y = sinx e2x , 0≤x≤0.5 
điều kiện ban đầu 
y(0) = -0.4, y’(0) = -0.6
với bước h = 0.1 
đặt y1 = y, y2 = y’ chuyển pt về hệ
y’1 = y2
y’2 = sinx e
2x– 2 y1 + 2y2
điều kiện y1(0) = -0.4, y2(0) = -0.6
Công thức Euler
y1 0 = -0.4
y1 k+1 = y1 k + 0.1 y2k 
y2 0 = -0.6
y2 k+1 = y2 k + 0.1 (sinxke
2xk - 2y1k +2y2 k)
xk y1 k=y y2 k=y’
0 -0.4 -0.6
0.1 -0.46 -0.64
0.2 -0.524 -0.6638
0.3 -0.5904 -0.6621
0.4 -0.6566 -0.6226
0.5 -0.7189 -0.5292
Ví dụ : Xét bài toán Cauchy
 x“(t) = (Mt+5) x2(t) – 2Mx’(t)+1.2t + M, 1≤t
điều kiện ban đầu 
x(1) = 1.3M, x’(1) = 1.8M
Dùng công thức Euler cải tiến, xấp xỉ giá trị của 
hàm x(t) và x’(t) tại điểm t = 1.2 với bước h = 0.2 
và M = 2.7
đặt y1 = x, y2 = x’ chuyển pt về hệ
y’1 = y2
y’2 = (Mt+5)y1
2-2My2+1.2t+M
điều kiện y1(1) = 1.3M, y2(1) = 1.8M
giải
Công thức Euler cải tiến
 y1 0 = 1.3M
 y1 1 = y1 0 + (K11+K21)/2
 y2 0 = 1.8M
 y2 1 = y2 0 + (K12+K22)/2
K11= 0.2*y2 0
K21= 0.2*(y2 0+K12)
K12= 0.2( (Mt0+5)y10
2 -2My20+ 1.2t0+M)
K22= 0.2( (M(t0+h)+5)(y10 +K11)
2 -2M(y20+K12)
 + 1.2(t0+h)+M)
K11 = 0.972 K12 =14.504154 
K21 = 3.8728308 K22 = 13.02027163
x(1.2) = y1 1 =5.9324154
x’(1.2) = y2 1=18.622212816
IV. GIẢI PTVP TUYẾN TÍNH CẤP 2 BẰNG 
PP SAI PHÂN HỮU HẠN : 
Xét phương trình vi phân tuyến tính cấp 2 với điều 
kiện biên
p(x)y” + q(x)y’ + r(x)y = f(x), a≤x≤b
y(a) = α, y(b) = β
❖ PP sai phân hữu hạn :
▪ Chia đoạn [a,b] thành n đoạn bằng nhau với 
bước h=(b-a)/n và các điểm nút
x0 = a, x1 = x0 +h, ... , xk = x0 + kh, ... , xn = b
▪ sử dụng các công thức sai phân hướng tâm ta 
xấp xỉ
y’(xk) ≈ (yk+1 – yk-1) /2h
y”( xk) ≈ (yk+1 – 2yk + yk-1)/h
2
với yk là giá trị xấp xỉ của hàm tại điểm xk.
▪ thay xk vào phương trình ta được 
pk (yk+1–2yk + yk-1)/h
2+ qk (yk+1–yk-1)/(2h) +rkyk= fk
với pk = p(xk), qk = q(xk), rk = r(xk), fk = f(xk),
▪ biến đổi phương trình trên ta thu được hệ 
phương trình sau :
Đây chính là hệ phương trình tuyến tính 
Ay = b
Với A là ma trận 3 đường chéo
Ví dụ : Giải gần đúng pt vi phân cấp 2
y“ - y’ – (x+1)y = x-1, 0≤x≤1
y(0) = 1, y(1) = 0
với bước h = 0.25
giải
n = 4
x0 = 0, x1 = 0.25, x2 = 0.5, x3 = 0.75, x4 = 1
k 1 2 3
qk -1 -1 -1
rk -1.25 -1.5 -1.75
fk -0.75 -0.5 -0.25
Vậy nghiệm gần đúng
y(0) =1, y(0.25) = 0.8093, y(0.5) = 0.5827, y(0.75)=0.3182, y(1)=0
giải hệ phương trình tuyến tính: Ay = b
Ví dụ : Giải gần đúng pt vi phân cấp 2
Mx2y”+xy’-8My = -3M2x2, 1.4≤x≤1.8
y(1.4) = 0.5M, y(1.8) = 1.5M
với bước h = 0.1, M= 2.7
giải
n = 4
x0 = 1.4, x1 = 1.5, x2 = 1.6, x3 = 1.7, x4 = 1.8
k 1 2 3
pk 6.075 6.912 7.803
qk 1.5 1.6 1.7
rk -21.6 -21.6 -21.6
fk -49.2075 -55.9872 -63.2043
giải hệ phương trình tuyến tính : Ay = b
Vậy nghiệm gần đúng
y(1.4) =1.35, y(1.5) = 2.0499, y(1.6) = 2.7247, y(1.7)=3.3882, y(1.8)=4.05