Wednesday, May 09, 2012 Người viết: Ôn Ngũ Minh
GIẢI GẦN ĐÚNG PHƯƠNG TRÌNH VI PHÂN
1. Phương trình vi phân cấp 1
Tìm nghiệm của phương trình y' = f(x, y) trên với a ≤ x ≤ b, với điều kiện y(a) = y0.
Cách giải xấp xỉ: Cho n là số tự nhiên dương.
Ký hiệu: h =(b – a)/n, xk = x0 + kh, fk = f(xk, yk) , k = 1..n.
Ví dụ: y' = cosx + 2siny, 0 ≤ x ≤ 0.4, y(0) = 1, n = 4 ⇒ h = 0.1.
1.1 Phương pháp Euler
a. Phương pháp Euler
y0 là giá trị khởi tạo đã cho,
Với k = 0..n-1: yk+1 = yk + hfk. Sai số là O(h).
Với ví dụ trên, ta lập bảng như sau:
k 0 1 2 3 4
x 0 0.1 0.2 0.3 0.4
y 1 1.268 1.559 1.857 2.144
Sử dụng máy tính Casio: Đặt cấu hình về radian.
0.1→M
0→X
1→Y
Y+M(cosX+2sinY)→Y (1) Ghi kết quả vào bảng mỗi khi thực hiện lệnh này,
X+M→X (2)
Lặp lần lượt hai lệnh (1) và (2) cho tới khi điền đầy bảng.
b. Phương pháp Euler cải tiến
y0 là giá trị khởi tạo đã cho,
Với k = 0..n-1, lặp quá trình sau:
yk+1(0) = yk + hfk,
yk+1(m+1) = yk + h/2*(fk + f(xk+1, yk+1(m))), m = 0, 1, 2, ...
Dừng khi | yk+1(m+1) – yk+1(m)| < ε. Sai số làm tròn là O(h2).
Với ví dụ trên, ta lập bảng như sau:
k 0 1 2 3 4
x 0 0.1 0.2 0.3 0.4
y 1.000 1.280 1.575 1.867 2.141
f(x,y) 2.683 2.911 2.980 2.868 2.605
1.268 1.571 1.873 2.154
1.279 1.575 1.867 2.140
1.280 1.575 1.867 2.141
1.280 2.141
Sử dụng máy tính Casio: Đặt cấu hình về radian.
0.1→M Giá trị của h = 0.1
0→X
1→Y
cosX+2sinY→F (1) Giá trị f(x, y) này được dùng nhiều lần
Y+MF→D (2) D là y
k+1(0)
Y+M÷2(F+cos(X+M)+2sinD)→D (3) D là yk+1(m+1), lặp liên tiếp
D→Y (4) Giá trị yk+1(m+1) đã thỏa mãn sai số
X+M→X (5) Tăng X
Lệnh (3) được lặp cho đến khi 2 giá trị liên tiếp của M bằng nhau. Sau khi thực hiện lệnh (5) thì
quay về lệnh (1).
1.2 Phương pháp Runge – Kutta
a. Phương pháp Runge – Kutta thứ nhất
y0 là giá trị khởi tạo đã cho,
Với k = 0..n-1:
a = hfk, b = hf(xk + h/2, yk + a/2) , yk+1 = yk + b.
Sai số làm tròn là O(h2).
Với ví dụ trên, ta lập bảng như sau:
k 0 1 2 3 4
x 0 0.1 0.2 0.3 0.4
y 1.000 1.281 1.578 1.872 2.146
a 0.268 0.291 0.298 0.287
b 0.281 0.297 0.294 0.274
Sử dụng máy tính Casio: Đặt cấu hình về radian.
0.1→M Giá trị của h = 0.1
0→X
1→Y
M(cosX+2sinY)→A (1)
M(cos(X+M÷2)+2sin(Y+A÷2))→B (2)
Y+B→Y (3)
X+M→X (4)
Lặp lại từ lệnh (1).
b. Phương pháp Runge – Kutta thứ hai
y0 là giá trị khởi tạo đã cho,
Với k = 0..n-1:
a = hfk, b = hf(xk+ h, yk + a), yk+1 = yk + (a + b)/2, k = 0..n-1.
Sai số làm tròn là O(h2).
Với ví dụ trên, ta lập bảng như sau:
k 0 1 2 3 4
x 0 0.1 0.2 0.3 0.4
y 1.000 1.281 1.577 1.871 2.146
a 0.268 0.291 0.298 0.287
b 0.281 0.297 0.294 0.275
Sử dụng máy tính Casio: Đặt cấu hình về radian.
0.1→M Giá trị của h = 0.1
0→X
1→Y
M(cosX+2sinY)→A (1)
M(cos(X+M)+2sin(Y+A))→B (2)
Y+(A+B) ÷2→Y (3)
X+M→X (4)
c. Phương pháp Runge – Kutta thứ ba
y0 là giá trị khởi tạo đã cho,
Với k = 0..n-1:
a = hfk, b = hf(xk + h/2, yk + a/2), c = hf(xk + h, yk + 2b - a)
yk+1 = yk + (a + 4b + c)/6, k = 0..n-1.
Sai số làm tròn là O(h3).
Với ví dụ trên, ta lập bảng như sau:
k 0 1 2 3 4
x 0 0.1 0.2 0.3 0.4
y 1.0000 1.2808 1.5769 1.8707 2.1451
a 0.2683 0.2911 0.2980 0.2866
b 0.2811 0.2968 0.2945 0.2746
c 0.2919 0.2980 0.2868 0.2613
Sử dụng máy tính Casio: Đặt cấu hình về radian.
0.1→M Giá trị của h = 0.1
0→X
1→Y
M(cosX+2sinY)→A (1)
M(cos(X+M÷2)+2sin(Y+A÷2))→B (2)
M(cos(X+M)+2sin(Y+2B-A))→C (3)
Y+(A+4B+C)÷6→Y (4)
X+M→X (5)
d. Phương pháp Runge – Kutta thứ tư
y0 là giá trị khởi tạo đã cho,
Với k = 0..n-1:
a = hfk, b = hf(xk + h/2, yk + a/2), c = hf(xk + h/2, yk + b/2)
d = hf(xk + h, yk + c), yk+1 = yk + (a + 2b + 2c + d)/6, k = 0..n-1.
Sai số làm tròn là O(h4).
Với ví dụ trên, ta lập bảng như sau:
k 0 1 2 3 4
x 0 0.1 0.2 0.3 0.4
y 1.00000 1.28084 1.57693 1.87070 2.14504
a 0.26829 0.29115 0.29800 0.28661
b 0.28111 0.29680 0.29449 0.27461
c 0.28165 0.29688 0.29454 0.27512
d 0.29120 0.29800 0.28656 0.25994
Sử dụng máy tính Casio: Đặt cấu hình về radian.
0.1→M Giá trị của h = 0.1
0→X
1→Y
M(cosX+2sinY)→A (1)
M(cos(X+M÷2)+2sin(Y+A÷2))→B (2)
M(cos(X+M÷2)+2sin(Y+B÷2))→C (3)
M(cos(X+M)+2sin(Y+C))→D (4)
Y+(A+2B+2C+D) ÷6→Y (5)
X+M→X (6)
Chú ý: Có thể không đủ bộ nhớ với một số máy, nên cần thử trước. Có thể ghép trực tiếp biểu
thức A vào lệnh (2) và lệnh (5) như sau:
0.1→M Giá trị của h = 0.1
0→X
1→Y
M(cos(X+M÷2)+2sin(Y+ M(cosX+2sinY)÷2))→B (2)
M(cos(X+M÷2)+2sin(Y+B÷2))→C (3)
M(cos(X+M)+2sin(Y+C))→D (4)
Y+(M(cosX+2sinY)+2B+2C+D) ÷6→Y (5)
X+M→X (6)
2. HỆ 2 PHƯƠNG TRÌNH VI PHÂN CẤP 1
y' = f(x, y, z) y(x0) = y0
z' = g(x, y, z) z(x0) = z0
Ký hiệu: xk = x0 + kh, fk = f(xk, yk, zk), gk = g(xk, yk, zk), k = 0, 2, ..., n
2. 1. Phương pháp Euler
y0, z0, là các giá trị khởi tạo đã cho,
Với k = 0..n-1:
yk+1 = yk + hfk,j, zk+1 = zk + hgk,
Sai số làm tròn là O(h).
Ví dụ:
y' = x + sinz
z' = y + cosz y(1) = 2, z(1) = 3, h = 0.1, n = 4.
Ta lập bảng như sau:
k 0 1 2 3 4
x 1 1.1 1.2 1.3 1.4
y 2 2.11 2.22 2.33 2.44
z 3 3.10 3.21 3.33 3.46
Sử dụng máy tính Casio: Đặt cấu hình về radian.
0.1→E
1→X
2→Y
3→M Thay cho biến Z
Y+E(X+sinM)→A (1)
M+E(Y+cosM)→M (2)
A→Y (3)
X+E→X (4)
Nếu lệnh (1) viết là X+sinM→Y thì giá trị Y mới này sẽ tham gia vào lệnh (2), không
đúng với giải thuật trên.
1.3 Phương pháp Runge – Kutta
a. Phương pháp Runge – Kutta thứ nhất
y0, z0, là các giá trị khởi tạo đã cho,
Với k = 0..n-1:
a = hf(xk, yk, zk), c = hg(xk, yk, zk),
b = hf(xk+ h/2, yk + a/2, zk + c/2),
d = hg(xk+ h/2, yk + a/2, zk + c/2)
yk+1 = yk + b, zk+1 = zk + d
Sai số làm tròn là O(h2).
Với ví dụ trên, ta lập bảng như sau:
k 0 1 2 3 4
x 1 1.1 1.2 1.3 1.4
y 2 2.114 2.227 2.338 2.446
z 3 3.106 3.223 3.352 3.495
a 0.114 0.114 0.112 0.109
b 0.114 0.113 0.111 0.108
c 0.101 0.111 0.123 0.136
d 0.106 0.117 0.129 0.143
Sử dụng máy tính Casio: Đặt cấu hình về radian.
0.1→E
1→X
2→Y
3→M Thay cho biến Z
E(X+sinM)→A (1)
E(Y+cosM)→C (2)
E(X+E÷2+sin(M+C÷2))→B (3)
E(Y+A÷2 + cos(M +C÷2))→D (4)
Y+B→Y (5)
M+D→M (6)
X+E→X (7)
Nếu số phép gán không được lưu đầy đủ, ta ghép trực tiếp (1) với (6) vào (4):
0.1→E
1→X
2→Y
3→M Thay cho biến Z
E(Y+cosM)→C (2)