Bài giảng Phương pháp số - Bài 6: Nghiệm của các phương trình vi phân - Nguyễn Thị Vinh

BÀI 6
NGHIỆM CỦA CÁC 
PHƢƠNG TRÌNH VI PHÂN
MỘT SỐ KHÁI NIỆM CƠ SỞ (1)
1. PHƢƠNG TRÌNH VI PHÂN (PTVP) CẤP n
y(n)(x) = f(x, y(x), y´(x) , , y(n–1)(x))
2. NGHIỆM CỦA PTVP CẤP n: 
Nghiệm của PTVP cấp n là một hàm Φ(x) khả vi liên tục đến
cấp n trên khoảng đang xét và
Φ(n)(x) = f(x, Φ(x), Φ’(x), Φ’’(x), , Φ(n-1)(x))
Nghiệm tổng quát của PTVP cấp n thƣờng chứa n hằng số
tùy ý, và do vậy tồn tại một họ n tham số của các nghiệm
Nghiệm của bài toán giá trị đầu: Tìm nghiệm riêng suy từ
nghiệm tổng quát của PTVP sao cho nó thỏa mãn n giá trị
đầu đã biết
y(x0), y’(x0), , y
(n–1)(x0)
PHƢƠNG PHÁP SỐ-Bài 6 2
MỘT SỐ KHÁI NIỆM CƠ SỞ (2)
3. NGHIỆM CỦA PTVP TUYẾN TÍNH CẤP n: 
• PTVP là tuyến tính cấp n
y(n)(x) + p1(x)y
(n–1)(x) + + pn-1(x)y´(x) + pn(x)y(x) = q(x)
• PTVP tuyến tính thuần nhất cấp n (khi vế phải q(x) = 0)
i> Nếu y1(x), y2(x), . . . , ym(x) là các nghiệm bất kì của PTVP
tuyến tính thuần nhất thì tổ hợp tuyến tính của chúng
C1y1(x) + C2y2(x) + · · · + Cmym(x) cũng là một nghiệm
Ii> Các nghiệm y1(x), y2(x), . . . , ym(x) này đƣợc gọi là độc
lập tuyến tính nếu định thức hàm
0
yyy
yyy
yyy
1m
mmm
1m
222
1m
111
)('
)('
)('
...
............
...
...
PHƢƠNG PHÁP SỐ-Bài 6 3
MỘT SỐ KHÁI NIỆM CƠ SỞ (3)
3. NGHIỆM CỦA PTVP TUYẾN TÍNH CẤP n: 
• PTVP tuyến tính thuần nhất cấp n (tiếp)
iii> Nếu y1(x), y2(x), . . . , yn(x) là n nghiệm độc lập tuyến tính
của PTVP thuần nhất cấp n, thì tổ hợp tuyến tính của chúng
Y(x) = C1y1(x) + C2y2(x) + · · · + Cnyn(x)
đƣợc gọi là nghiệm tổng quát.
• PTVP tuyến tính thuần nhất cấp n với các hệ số hằng
y(n) + an–1y
(n–1) + · · · + a0y = 0 (*)
có nghiệm dạng eβx với β thỏa mãn PT đặc tính
βn + an–1β
n–1 + · · · + a0 = 0 (**)
PHƢƠNG PHÁP SỐ-Bài 6 4
MỘT SỐ KHÁI NIỆM CƠ SỞ (4)
3. NGHIỆM CỦA PTVP TUYẾN TÍNH CẤP n: 
• PTVP tuyến tính thuần nhất cấp n với các hệ số hằng:
- Nếu phƣơng trình (**) có n nghiệm thực phân biệt βi thì
là nghiệm tổng quát của của PTVP thuần nhất cấp n (*)
- Nếu β1 = α + iβ là một nghiệm phức của (**), 
thì β2 = α – iβ cũng là một nghiệm của PT đặc tính (**)
- Nếu β là một nghiệm kép của PT đặc tính (**), thì
cũng là các nghiệm độc lập tuyến tính của PT thuần nhất (*)
xnβ
n
xβ
2
xβ
1 eC...eCeCY(x)
21 
βx
2
βx
1 xeyvàey 
PHƢƠNG PHÁP SỐ-Bài 6 5
MỘT SỐ KHÁI NIỆM CƠ SỞ (5)
3. NGHIỆM CỦA PTVP TUYẾN TÍNH CẤP n: 
• PTVP tuyến tính không thuần nhất cấp n hệ số hằng: 
y(n) + an–1y
(n–1) + · · · + a0y = q(x) (***)
Nếu PTVP (***) có một nghiệm riêng ζ(x), thì nghiệm tổng
quát của PTVP tuyến tính không thuần nhất là
với giả thiết rằng βi là các nghiệm thực phân biệt của PTVP 
đặc tính (**) tƣơng ứng
xnβ
n
xβ
2
xβ
1 eC...eCeCζ(x)
xYζ(x)xy
21 
)()(
PHƢƠNG PHÁP SỐ-Bài 6 6
MỘT SỐ KHÁI NIỆM CƠ SỞ (6)
3. NGHIỆM CỦA PTVP TUYẾN TÍNH CẤP n: 
• PTVP tuyến tính không thuần nhất cấp n hệ số hằng:
Ví dụ: Tìm nghiệm tổng quát của phƣơng trình sau
y´´ – 4y´ + 3y = x 
Giải: Từ PT đặc tính β2 - 4β + 3 = 0 β1 = 1, β2 = 3
 nghiệm tổng quát của PT thuần nhất là
Ta tìm một nghiệm riêng y1(x) = ax + b của PT ban đầu
 a = 1/3, b = 4/9
Vậy nghiệm tổng quát của PT không thuần nhất là
x3
2
x
1 eCeCxY )(
PHƢƠNG PHÁP SỐ-Bài 6 7
x3
2
x
11 eCeC
9
4
x
3
1
xYxyxy )()()(
MỘT SỐ KHÁI NIỆM CƠ SỞ (7)
4. BÀI TOÁN CÔ SI CỦA PTVP CẤP n
Tìm nghiệm riêng của PTVP (phi tuyến) cấp n
y(n) = f(x, y, y’, , y(n-1)) 
thỏa mãn các giá trị đầu
PHƢƠNG PHÁP SỐ-Bài 6 8
)()(
)(,,')(',)(
1n
00
1n
0000 yxyyxyyxy
 
Ví dụ 1: Tìm nghiệm riêng của PTVP y’ = xy thỏa mãn
đ/k giá trị đầu y(0) = 1
Giải: . 
dxx
y
dy
xy
dx
dy
 tích phân 2 vế ta có nghiệm TQ
2
2x
AeyC
2
x
|y|lndxx
y
dy 2
Từ đ/k y(0) = 1 ta tìm đƣợc nghiệm riêng 2
2x
ey 
MỘT SỐ KHÁI NIỆM CƠ SỞ (8)
5. THUẬT TOÁN TAYLOR
Giải BT Côsi: Tìm nghiệm của PTVP cấp một
y´ = f(x,y) biết giá trị đầu y(x0) = y0 (****)
• Giả sử f là đủ khả vi theo cả x và y, (****) có nghiệm duy
nhất nếu ∂f/∂y liên tục trên miền đang xét. Khai triển
nghiệm y(x) thành chuỗi Taylor xung quanh điểm x0:
• Tính các đạo hàm bằng cách lấy đạo hàm toàn phần của
(****) theo x: y´ = f(x,y)
y´´ = f´ = fx + fyy´ = fx + fyf
y´´´ = f´´ = fxx + fxyf + fyxf + fyyf
2 + fyfx + fy
2f
= fxx + 2fxyf + fyyf
2 + fxfy + fy
2f
 
 )(''
!
)(
)(')()( 0
2
0
000 xy
2
xx
xyxxyxy
PHƢƠNG PHÁP SỐ-Bài 6 9
MỘT SỐ KHÁI NIỆM CƠ SỞ (9)
5. THUẬT TOÁN TAYLOR (tiếp)
Thuật toán Taylor bậc k tìm một nghiệm xấp xỉ của BT Côsi
y´ = f(x,y) thỏa mãn đ/k đầu y(a) = y0 trên [a, b]
1. Chọn bƣớc h = (b – a) / N, đặt xn = a + nh, n = 0, 1, . . . , N
2. Tạo ra các xấp xỉ yn cho y(xn) từ phép đệ quy
yn+1 = yn + hTk(xn,yn ) n = 0, l, . . . , N – 1
trong đó Tk(x, y) đƣợc định nghĩa bởi
Sai số địa phƣơng của thuật toán Taylor bậc k là
...,,),,(
!
),('
!
),(),(
)( 21kyxf
k
h
yxf
2
h
yxfyxT 1k
1k
k 

hxxy
1k
h
1k
yfh
E nn
1k
1kk1k


 , )(
)!()!(
))(,( )(
)(
PHƢƠNG PHÁP SỐ-Bài 6 10
MỘT SỐ KHÁI NIỆM CƠ SỞ (10)
5. THUẬT TOÁN TAYLOR (tiếp)
Ví dụ: tìm một nghiệm xấp xỉ của BT Côsi
trên [1, 2] bằng thuật toán Taylor bậc 2 với h = 0.1
Đặt xi = a + ih, i = 0, 1, . . . , 10
do đó
và
 thuật toán Taylor bậc cao hơn khó đƣợc chấp nhận
11ykđmãnathy
x
y
x
1
y 2
2
 )(/' á
'
'
),(',),( yy2
x
y
x
y
x
2
yxfy
x
y
x
1
yxf
23
2
2
'),( f
2
h
fyxT2 
 ),('),( iiiii1i yxf
2
h
yxfhyy
PHƢƠNG PHÁP SỐ-Bài 6 11
MỘT SỐ KHÁI NIỆM CƠ SỞ (11)
5. THUẬT TOÁN TAYLOR (tiếp) 
Ví dụ: (tiếp)
i xi yi y'=f(xi,yi) hf(xi,yi) f'(xi,yi) h^2*f'(x,y)/2
0 1 -1 1 0.1 -2 -0.010000
1 1.1 -0.910000 0.825619 0.082562 -1.502632 -0.007513
2 1.2 -0.834951 0.693094 0.069309 -1.157414 -0.005787
3 1.3 -0.771429 0.590020 0.059002 -0.910343 -0.004552
4 1.4 -0.716979 0.508273 0.050827 -0.728879 -0.003644
5 1.5 -0.669796 0.442349 0.044235 -0.592612 -0.002963
6 1.6 -0.628524 0.388410 0.038841 -0.488305 -0.002442
7 1.7 -0.592124 0.343718 0.034372 -0.407110 -0.002036
8 1.8 -0.559788 0.306273 0.030627 -0.342966 -0.001715
9 1.9 -0.530876 0.274588 0.027459 -0.291621 -0.001458
10 2 -0.504875 0.247539 0.024754 -0.250036 -0.001250
PHƢƠNG PHÁP SỐ-Bài 6 12
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (1)
1. PHƢƠNG PHÁP EULER GIẢI BT CÔSI CẤP 1
• Áp dụng thuật toán Taylor với k = 1 trên [a, b]
có sai số địa phƣơng
 sai số toàn cục trên [a, b] cùng bậc với h cần h rất nhỏ
• Ví dụ, giải BT Côsi y´ = y biết giá trị đầu y(0) = 1, h = 0.01
),( iii1i yxhfyy 
)('' y
2
h
E
2
i xi yi f(x,y) hf(x,y)
0 0 1 1.000000 0.010000
1 0.01 1.010000 1.010000 0.010100
2 0.02 1.020100 1.020100 0.010201
3 0.03 1.030301 1.030301 0.010303
4 0.04 1.040604 1.040604 0.010406
5 0.05 1.051010 1.051010 0.010510
PHƢƠNG PHÁP SỐ-Bài 6 13
void Euler (double a, double b, double h, double y0) {
int n = ceil((b-a)/h); vector x(n+1,0), y(n+1,0);
x[0] = a; y[0] = y0;
for (int i = 1; i <= n; i++) {
x[i] = x[i-1] + h;
y[i] = y[i-1] + h * f(x[i-1], y[i-1]);
}
cout << "Nghiem gan dung tren doan [" << x[0]
<< "," << x[n] << "]" << endl;
cout << "\t i \t x \t y \t f(x,y) \t hf(x,y)" << endl;
for (int i = 0; i <= n; i++) {
cout << "\t "<< i << "\t " << x[i] << "\t " << y[i] << "\t "
<< f(x[i], y[i]) << "\t " << h*f(x[i], y[i]) << endl;
}
}
PHƢƠNG PHÁP SỐ-Bài 6 14
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (2)
2. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 2
• Đặt vấn đề: Đối với nghiệm của BT Côsi của PTVP cấp 1 
y’ = f(x, y), y(x0) = y0
Ƣớc lƣợng hàm f(x, y) tại các điểm
đƣợc lựa chọn trên mỗi khoảng con. 
yi+l = yi + ak1(h) + bk2(h)
trong đó a = b = 1/2
k1 = hf(xi, yi)
k2 = hf(xi + h, yi + k1) 
Sử dụng CT 
Euler, tính độ
nghiêng ở xi+1
yi+l = yi + 1/2(k1 + k2)
• Sai số địa phƣơng
)()()(
42
yyx
2
yyxyxx
3
1i1i hOffff2ffff2f
12
h
yxy 
PHƢƠNG PHÁP SỐ-Bài 6 15
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (3)
2. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 2
• Ví dụ: giải PTVP y'=1-y/x biết y(2) = 2, h=0.01trên [2; 2.1]
i xi yi f(xi,yi) hf(xi,yi) zi f(xi+1,zi) hf(xi+1,zi)
0 2.00 2.000000 0.000000 0.000000 2.000000 0.004975 4.9751E-05
1 2.01 2.000025 0.004963 0.000050 2.000075 0.009864 9.8641E-05
2 2.02 2.000099 0.009852 0.000099 2.000198 0.014681 0.00014681
3 2.03 2.000222 0.014669 0.000147 2.000368 0.019427 0.00019427
4 2.04 2.000392 0.019416 0.000194 2.000586 0.024104 0.00024104
5 2.05 2.000610 0.024093 0.000241 2.000851 0.028713 0.00028713
6 2.06 2.000874 0.028702 0.000287 2.001161 0.033256 0.00033256
7 2.07 2.001184 0.033245 0.000332 2.001516 0.037733 0.00037733
8 2.08 2.001538 0.037722 0.000377 2.001916 0.042146 0.00042146
9 2.09 2.001938 0.042135 0.000421 2.002359 0.046496 0.00046496
10 2.10 2.002381 0.046485 0.000465 2.002846
PHƢƠNG PHÁP SỐ-Bài 6 16
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (4)
void Runge_Kutta2 (double h, double x, double y, unsigned n) {
double xs, z;
cout << "Nghiem cua phuong trinh " << endl;
cout << "i x y f(x,y) f(x,z)";
for (unsigned i = 0; i <= n; i++) {
z = y + h * f(x, y);
xs = x + h;
cout << "\n“ << i << "\t “ << xs << "\t “ << y 
<< "\t “ << f(x, y) << "\t “ << f(xs, z) << endl;
y = y + h * (f(x, y) + f(xs, z)) / 2;
x = xs;
}
}
PHƢƠNG PHÁP SỐ-Bài 6 17
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (5)
3. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 4
Đối với PTVP y’ = f(x, y), y(x0) = y0, tạo xấp xỉ cho y(x0+h) 
y(x0+h) = y(x0) + r1k1(h) + r2k2(h) + r3k3(h) + r4k4(h) + φ4(h)
hay y(x0+h) = y0 + ∆y0 + φ4(h)
trong đó
∆y0 = r1k1(h) + r2k2(h) + r3k3(h) + r4k4(h) 
ki(h) = hf(ξi, ૪i); 
ξi = x0 + ai h, a1 = 0
૪i = y0 + bi1k1(h) +  + bi, i-1 ki-1(h), i = 1, 2, 3, 4 
φ4(h) = y(x0+h) - y0 - ∆y0 là sai số (địa phƣơng ) 
Các hệ số ai, bij, ri đƣợc chọn sao cho các đạo hàm cấp i của 
φ4(h) thoả mãn điều kiện 
4i0)0( và )4,0 (i0)0(
)1(i
4
(i)
4  
PHƢƠNG PHÁP SỐ-Bài 6 18
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (6)
19
trong đó
 4321i1i kk2k2k
6
1
yy 
)kyh,hf(xk
k
2
1
yh,
2
1
xhfk
k
2
1
yh,
2
1
xhfk
)y,hf(xk
3ii4
2ii3
1ii2
ii1
3. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 4
 hệ 11 phƣơng trình, 13 ẩn
a2, a3, a4, b21, b31, b32, b41, b42, b43, r1, r2, r3, r4.
• Công thức Runge – Kutta bậc 4 :
PHƢƠNG PHÁP SỐ-Bài 6
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (7)
20
yi+1
xi
yi
xi+1
y
xxi+1/2
k1
k4ktổ hợp
k3
k2
3. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 4
• Sai số địa phƣơng của công thức là O(h4)
• Ƣu điểm: Tự bắt đầu (chỉ đòi hỏi giá trị y(x0) = y0) 
• Nhƣợc điểm: Mỗi bƣớc đòi hỏi 4 ƣớc lƣợng hàm 
PHƢƠNG PHÁP SỐ-Bài 6
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (8)
i xi yi k1=hf(xi,yi) k2=hf(xi+h/2,yi+k1/2) k3=hf(xi+h/2,yi+k2/2) k4=hf(xi+h,yi+k3) Delta
0 1 2 0.1
1.1 2.050000 0.122439
1.1 2.061220 0.122970
1.2 2.122970 0.145792 0.122768
3. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 4
Ví dụ: Tìm nghiệm của BT Côsi
y´ = x – 1/y với điều kiện đầu y(1) = 2 trên [1; 1.4], h = 0.2 
 y2 = y(1.4) ≈ y(1.2) + Δy1 = 2.291957
PHƢƠNG PHÁP SỐ-Bài 6 21
 y1 = y(1.2) ≈ y(1) + Δy0 = 2.122768
1 1.2 2.122768 0.145783
1.3 2.195660 0.168911
1.3 2.207224 0.169388
1.4 2.292157 0.192746 0.169188
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (9)
void Runge_Kutta4 (double h, double x, double y, unsigned n) {
double k1, k2, k3, k4, delta;
cout << " Nghiem cua phuong trinh y'= 1+y/x" << endl;
cout << "\t i x y delta“ << endl;
for (unsigned i = 1; i <= n; i++) {
k1 = h * f(x, y);
k2 = h * f(x + h / 2, y + k1 / 2);
k3 = h * f(x + h / 2, y + k2 / 2);
k4 = h * f(x+h, y+k3);
delta= (k1 + 2 * k2 + 2 * k3 + k4) / 6;
x += h;
y += delta;
cout << "\t " << i << "\t " << x << "\t " << y << "\t " << delta 
<< endl;
}
} 
PHƢƠNG PHÁP SỐ-Bài 6 22
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (10)
4. BÀI TOÁN CÔSI CỦA HỆ PTVP CẤP 1
Các phƣơng pháp trên có thể áp dụng đƣợc cho hệ PTVP 
cấp một. Chẳng hạn giải bài toán Côsi sau đây: 
Cho hệ hai PTVP
)z,y,x(g'z
)z,y,x(f'y
với điều kiện đầu y(x0) = y0 và z(x0) = z0, bằng phƣơng
pháp Runge – Kutta bậc 4. 
Khi đó công thức Runge – Kutta bậc 4 có dạng
yi+1≈yi+(k1+2k2+2k3+k4) / 6 
zi+1≈zi+(l1+2l2+2l3+l4) / 6
PHƢƠNG PHÁP SỐ-Bài 6 23
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (11)
k1 = hf(xi, yi, zi)
k2= hf(xi+h/2, yi+k1/2, zi+l1/2)
k3 = hf(xi+h/2, yi+k2/2, zi+l2/2)
k4 = hf(xi+1, yi+k3, zi+l3)
l1 = hg(xi, yi, zi)
l2= hg(xi+h/2, yi+k1/2, zi+l1/2)
l3 = hg(xi+h/2, yi+k2/2, zi+l2/2)
l4 = hg(xi+1, yi+k3, zi+l3)
4. BÀI TOÁN CÔSI CỦA HỆ PTVP CẤP 1
Các hệ số ki và li đƣợc tính nhƣ sau
Ví dụ: Giải BT Côsi của hệ dƣới đây trên [0; 0,4] với h = 0,1
với các điều kiện đầu y(0) = 0 và z(0) = 1
yzxx'z
z'y
PHƢƠNG PHÁP SỐ-Bài 6 24
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (12)
4. BÀI TOÁN CÔSI CỦA HỆ PTVP CẤP 1
PHƢƠNG PHÁP SỐ-Bài 6 25
i x y z k=hz l = h(x+xz−y)
0 0 0 1 0,1 0
0,05 0,05 1 0,1 0.005
0,05 0,05 0,0025 0,10025 0.005013
0,1 0.10025 1.005013 0,100501 0.010025
y1 = y0 + ∆y0 = 0 + 0,100167 = 0.100167
z1 = z0 + ∆z0 = 1 + 0.005008 = 1.005008
1 0,1 0,100167 1.005008 0,100501 0.010033
0,15 0,150167 1.010025 0,101003 0.015134
0,15 0,150668 1.012575 0,101258 0.015122
0,20 0,201424 1.020130 0,102013 0.020260
y2 = y1 + ∆y1 = 0.100167 + 0.101172 = 0.201339
z2 = z1 + ∆z1 = 1.005008 + 0.015134 = 1.020142
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (13)
5. BÀI TOÁN CÔSI CỦA PTVP CẤP CAO
Đƣa về bài toán Cô si của hệ PTVP cấp một
Ví dụ: Giải BT Côsi của phƣơng trình cấp 2
y´´ – y´x + y = x 
thỏa mãn các điều kiện đầu
y(0) = 0 và y´(0) = 1
Cách giải: Đƣa về BT Côsi của hệ 2 PTVP cấp một
tƣơng đƣơng
với các đ/k đầu y(0) = 0 và z(0) = 1. Giải hệ này tìm 
được y
yzxx'z
z'y
PHƢƠNG PHÁP SỐ-Bài 6 26
CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (14)
BÀI TẬP
1. Bài tập tính toán: 
8.1-1 (a, b, d), 8.1-2, 8.3-3, 8.4-1 (b và d), 8.5-1,
8.12-1.
Dùng phƣơng pháp Runge –Kutta bậc 4 giải lại
bài tập 8.1-2 (a, c) trên [0; 0.3] với h = 0.1
2. Bài tập lập trình:
Phƣơng pháp Runge –Kutta bậc 2 và bậc 4 giải
BT Côsi của hệ 2 PTVP cấp 1
PHƢƠNG PHÁP SỐ-Bài 6 27