Chương 7 : Các phương trình vi phân thường
Một phương trình vi phân cấp 1 có thể viết dưới dạng giải được y′ = f(x,y) mà ta có thể tìm được hàm y từ đạo hàm của nó. Tồn tại vô số nghiệm thoả mãn phương trình trên. Mỗi nghiệm phụ thuộc vào một hằng số tuỳ ý.
CHƯƠNG 7: CÁC PHƯƠNG TRÌNH VI PHÂN THƯỜNG
§1. BÀI TOÁN CAUCHY
Một phương trình vi phân cấp 1 có thể viết dưới dạng giải được
y′ = f(x,y) mà ta có thể tìm được hàm y từ đạo hàm của nó. Tồn tại vô số
nghiệm thoả mãn phương trình trên. Mỗi nghiệm phụ thuộc vào một hằng số
tuỳ ý. Khi cho trước giá trị ban đầu của y là yo tại giá trị đầu xo ta nhận được
một nghiệm riêng của phương trình. Bài toán Cauchy (hay bài toán có điều
kiện đầu) tóm lại như sau: cho x sao cho b ≥ x ≥ a, tìm y(x) thoả mãn điều
kiện:
⎧y ′( x) = f( x , y)
⎨ (1)
⎩ y(a ) = α
Người ta chứng minh rằng bài toán này có một nghiệm duy nhất nếu f
thoả mãn điều kiện Lipschitz:
f( x , y 1 ) − f( x , y 2 ) ≤ L y 1 − y 2
với L là một hằng số dương.
Người ta cũng chứng minh rằng nếu f′y ( đạo hàm của f theo y ) là liên
tục và bị chặn thì f thoả mãn điều kiện Lipschitz.
Một cách tổng quát hơn, người ta định nghĩa hệ phương trình bậc 1:
y′ = f1 ( x , y 1 , y 2 ,..., y n )
1
y′2 = f2 ( x , y 1 , y 2 ,..., y n )
⋅⋅⋅⋅
y′n = fn ( x , y 1 , y 2 ,..., y n )
Ta phải tìm nghiệm y1, y2,..., yn sao cho:
⎧Y′( x) = f( x , X )
⎨
⎩Y(a ) = α
với:
⎛ y′ ⎞ ⎛ f1 ⎞ ⎛ y1 ⎞
⎜ 1⎟ ⎜ ⎟ ⎜ ⎟
⎜ y′2 ⎟ ⎜ f2 ⎟ ⎜ y2 ⎟
Y′ = ⎜ .. ⎟ F = ⎜ .. ⎟ Y = ⎜ .. ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ .. ⎟ ⎜ .. ⎟ ⎜ .. ⎟
⎜ ′⎟ ⎜ ⎟ ⎜y ⎟
⎝ yn ⎠ ⎝ fn ⎠ ⎝ n⎠
360
Nếu phương trình vi phân có bậc cao hơn (n), nghiệm sẽ phụ thuộc vào
n hằng số tuỳ ý. Để nhận được một nghiệm riêng, ta phải cho n điều kiện đầu.
Bài toán sẽ có giá trị đầu nếu với giá trị xo đã cho ta cho y(xo), y′(xo), y″(xo),....
Một phương trình vi phân bậc n có thể đưa về thành một hệ phương
trình vi phân cấp 1. Ví dụ nếu ta có phương trình vi phân cấp 2:
⎧y′′ = f( x , y , y′)
⎨
⎩y(a ) = α , y′(a ) = β
Khi đặt u = y và v = y′ ta nhận được hệ phương trình vi phân cấp 1:
⎧u′ = v
⎨
⎩ v′ = g( x , u , v )
với điều kiện đầu: u(a) = α và v(a) = β
Các phương pháp giải phương trình vi phân được trình bày trong
chương này là các phương pháp rời rạc: đoạn [a, b] được chia thành n đoạn
nhỏ bằng nhau được gọi là các ʺbướcʺ tích phân h = ( b ‐ a) / n.
§2. PHƯƠNG PHÁP EULER
Giả sử ta có phương trình vi phân:
⎧y ′( x) = f( x , y)
⎨ (1)
⎩y(a ) = α
và cần tìm nghiệm của nó. Ta chia đoạn [xo,x ] thành n phần bởi các điểm
chia:
xo function [X, Y] = euler(fxy, xo, xf, yo, n)
% %Giai phuong trinh yʹ(x) = f(x,y(x)) hay y’ = f(x)
if n 1
k1 = h*feval(fxy, x, y);
else
k1 = h*feval(fxy, x);
end
y = y + k1;
x = x + h;
X(i+1) = x;
Y(i+1, :) = yʹ;
end
function dy = f1(t, y)
dy = zeros(3, 1);
dy(1) = y(2) * y(3);
dy(2) = ‐y(1) * y(3);
dy(3) = ‐0.51 * y(1) * y(2);
Để giải phương trình cho bởi hàm f1(x, y) ta dùng chương trình cteuler.m:
clear all, clc
a = 0;
362
b = 1;
y = @f1;
ya = [0 1 1]ʹ;
m = 200;
[x, y] = euler(y, a, b, ya, m)
plot(x, y);
§3. PHƯƠNG PHÁP HEUN
Phương pháp Heun còn được gọi là phương pháp hình thang hay
phương pháp. Cho phương trình:
y’ = f(t, y)
Ta có:
t k +1
∫ f(t,y)dt
t k +1
yt = y(t k+1 ) − y(t k ) =
k
tk
hay:
t k +1
y(t k+1 ) = y(t k ) + ∫ f(t,y)dt
tk
với y(t0) = y0
Nếu ta sử dụng quy tắc tích phân hình thang thì ta có:
y k+1 = y k + ⎡f(t k ,y k ) − f(t k +1 , y k+1 ) ⎤
h
2⎢⎣ ⎥
⎦
Vế phải (RHS) của phương trình này có yk+1 là gái trị chưa biết tại thời điểm tk.
Để giải quyết vấn đề này ta thay yk+1 ở RHS bằng công thức xấp xỉ:
y k+1 ≅ y k + hf(t k ,y k )
Như vậy:
h
{
y k +1 = y k + f(t k ,y k ) + f [(t k +1 , y k + hf(t k ,y k )]
2 }
Đây chính là công thức Heun.
Ta xây dựng hàm heun() để thực hiện thuật toán trên:
function [X, Y] = heun(fxy, xo, xf, yo, n)
%Giai phuong trinh yʹ(x) = f(x,y(x)) hay y’ = f(x)
%dung thuat toan Heun voi n buoc tinh
if n h = (xf ‐ xo)/n;
X = zeros(n+1, 1);
M = max(size(yo));% so phuong trinh (so cot cua ma tran Y)
Y = zeros(n+1, M);
%dat dieu kien dau
x = xo;
X(1) = x;
y = yo;
Y(1,:) = yʹ;
for i = 1:n
if nargin(fxy) > 1
f1 = feval(fxy, x, y);
f2 = feval(fxy, x+h, y+h*f1);
else
f1 = feval(fxy, x);
f2 = feval(fxy, x+h);
end
y = y + h*(f1 + f2)/2;
x = x + h;
X(i+1) = x;
Y(i+1, :) = y.ʹ;
end
Để giải phương trình ta dùng chương trình ctheun.m:
clear all, clc
a = 0;
b = 1;
y = inline(ʹ2*x + yʹ);
ya = 0.5;
n = 10;%so lan tinh chi n = 10
[x, y] = heun(y, a, b, ya, n)
plot(x, y);
§4. PHƯƠNG PHÁP RUNGE ‐ KUTTA
364
Mặc dù phương pháp Heun tốt hơn phương pháp Euler nhưng nó vẫn
chưa đủ độ chính xác đối với các bài toán thực tế.
Xét bài toán Cauchy (1). Giả sử ta đã tìm được giá trị gần đúng yi của
y(xi) và muốn tính yi+1 của y(xi+1). Trước hết ta viết công thức Taylor:
h2 h m (m) h m + 1 ( m + 1)
y( x i +1 ) = y( x i ) + hy′( x i ) + y′′( x i ) + ⋅ ⋅ ⋅ + y (x i ) + y (c) (11)
2 m! m!
với c ∈(xi, xi+1) và:
y′(x i ) = f[x i , y( x i )]
d k −1
y ( x i ) = k −1 f [x i , y( x i )]
( k)
dx
Ta viết lại (11) dưới dạng:
h2 h m (m) h m + 1 ( m + 1)
y i +1 − y i = hy ′( x i ) + y ′′( x i ) + ⋅ ⋅ ⋅ + y (x i ) + y (c) (12)
2 m! m!
Ta đã kéo dài khai triển Taylor để kết quả chính xác hơn. Để tính y′i, y″i v.v. ta
có thể dùng phương pháp Runge‐Kutta bằng cách đặt:
y i + 1 − y i = r1 k (1i ) + r2 k (2i ) + r3 k (3i ) + r4 k (4i ) (13)
trong đó:
⎧k (1i ) = hf( x i , y i )
⎪ (i)
⎪k 2 = hf( x i + ah , y i + αk 1 )
(i)
⎨ (i) (14)
⎪ k 3 = hf( x i + bh , y i + β k (1i ) + γk (2i ) )
⎪.......
⎩
và ta cần xác định các hệ số a, b,..; α, β, γ,...; r1, r2,.. sao cho vế phải của (13)
khác với vế phải của (12) một vô cùng bé cấp cao nhất có thể có đối với h.
Khi dùng công thức Runge‐Kutta bậc hai ta có:
⎧k (1i ) = hf( x i , y i )
⎪
⎨ (i) (15)
⎪k 2 = hf( x i + ah , y i + αk (1i ) )
⎩
và y i +1 − y i = r1 k (1i ) + r2 k (2i ) (16)
Ta có:
y′(x) = f[x,y(x)]
y′′( x) = fx [x , y( x)] + fy [x , y( x)]
′ ′
................
Do đó vế phải của (12) là:
hf( x i , y i ) +
h2
2
[ ′ ]
fx ( x i , y i ) + fy ( x i , y i ) y′( x) + ⋅ ⋅ ⋅
′ (17)
Mặt khác theo (15) và theo công thức Taylor ta có:
365
k (1i ) = hf( x i , y i ) = hy′
i
′ ′
k 2 = h[f( x i , y i ) + ahfx ( x i , y i ) + αk (1i ) fy ( x i , y i ) + ⋅ ⋅ ⋅]
(i)
Do đó vế phải của (16) là:
h(r1 + r2 )f( x i , y i ) + h 2 [ar2 fx ( x i , y i ) + αr2 y′ fy ( x i , y i )] + ⋅ ⋅ ⋅
′ i
′ (18)
Bây giờ cho (17) và (18) khác nhau một vô cùng bé cấp O(h3) ta tìm được các
hệ số chưa biết khi cân bằng các số hạng chứa h và chứa h2:
r1 + r2 = 1
a.r1 = 1/ 2
α.r2 = 1
Như vậy: α = a, r1 = (2a ‐ 1)/ 2a, r2 = 1/ 2a với a được chọn bất kì.
Nếu a = 1 / 2 thì r1 = 0 và r2 = 1. Lúc này ta nhận được công thức Euler. Nếu
a=1 thì r1 = 1 / 2 và r2 = 1/2. Lúc này ta nhận được công thức Euler cải tiến.
Một cách tương tự chúng ta nhận được công thức Runge ‐ Kutta bậc 4.
Công thức này hay được dùng trong tính toán thực tế :
k1 = h.f(xi, yi)
k2 = h.f(xi+h/ 2, yi + k1/ 2)
k3 = h.f(xi+h/ 2, yi + k2/ 2)
k4 = h.f(xi+h, yi + k3)
yi+1 = yi + (k1 + 2k2 + 2k3 + k4) / 6
Ta xây dựng hàm rungekutta() để thực hiện công thức Runge ‐ Kutta bậc 4:
function [x, y] = rungekutta(f, a, b, y0, n)
%Phuong phap Runge‐Kutta de giai phuong trinh yʹ(x) = f(x,y(x)) hay y’ =
%f(x)
if nargin f1 = h*feval(f, x(k), y(k, :));
f1 = f1(:)ʹ;
f2 = h*feval(f, x(k) + h/2, y(k, :) + f1/2);
f2 = f2(:)ʹ;
f3 = h*feval(f, x(k) + h/2, y(k, :) + f2/2);
f3 = f3(:)ʹ;
f4 = h*feval(f, x(k) + h, y(k, :) + f3);
f4 = f4(:)ʹ;
y(k+1, :) = y(k, :) + (f1 + 2*(f2 + f3) + f4)/6;
end
else
for k = 1:n
f1 = h*feval(f, x(k));
f1 = f1(:)ʹ;
f2 = h*feval(f, x(k) + h/2);
f2 = f2(:)ʹ;
f3 = h*feval(f, x(k) + h/2);
f3 = f3(:)ʹ;
f4 = h*feval(f, x(k) + h);
f4 = f4(:)ʹ;
y(k+1, :) = y(k, :) + (f1 + 2*(f2 + f3) + f4)/6;
end
end
Để giải phương trình ta dùng chương trình ctrungekutta.m:
clear all, clc
a = 0;
b = 1;
y = inline(ʹx + yʹ);
ya = 0.5;
n = 10;%so lan tinh chi n = 10
[x, y] = rungekutta(y, a, b, ya, n)
plot(x, y);
§5. PHƯƠNG PHÁP RUNGE ‐ KUTTA THÍCH NGHI
367
Vấn đề xác định bước tính h là rất quan trọng. Nếu muốn có độ chính
xác cao thì bước tính h phải nhỏ. Tuy nhiên khi h nhỏ, ta lại tốn thời gian tính
toán. Hơn nữa bước hằng số sẽ không thích hợp trên toàn bộ miền tìm
nghiệm. Ví dụ nếu đường cong nghiệm ban đầu thay đổi nhanh rồi sau đó
gần như không đổi thì ta phải dùng h nhỏ ở đoạn đầu và h lớn ở đoạn sau.
đây là chỗ mà các phương pháp thích nghi chiếm ưu thế. Chúng đánh giá sai
số làm tròn tại mỗi lần tích phân và tự động hiệu chỉnh độ lớn của h để sai số
nằm trong giới hạn cho phép.
Phương pháp Runge ‐ Kutta thích nghi còn gọi là phương pháp tích
phân kết hợp. Các công thức này đi thành cặp: một công thức tích phân bậc m
và một công thức tích phân bậc m+1. Ý tưởng là dùng hai công thức này cải
thiện nghiệm trong đoạn [x, x+h]. Gọi kết quả là ym(x+h) và ym+1(x+h) ta có sai
số đối với công thức bậc m là:
E(h) = ym+1(x+h) ‐ ym(x+h) (1)
Chúng ta dùng công thức kết hợp bậc 4 và 5 mà đạo hàm được tính bằng
công thức Fehlenberg. Do vậy công thức Runge ‐ Kutta thích nghi còn được
gọi là công thức Runge ‐ Kutta ‐ Fehlenberg:
K 1 = hF(x,y)
⎛ i −1 ⎞
K i = hF ⎜ x + A i h, y + ∑ Bi ,jK j ⎟ i = 1, 2,..,6 (2)
⎝ j= 0 ⎠
6
y 5 (x + h) = y(x) + ∑ CiK i (công thức bậc 5) (3)
i =1
6
y 4 (x + h) = y(x) + ∑ Di K i (công thức bậc 4) (4)
i =1
Các hệ số xuất hiện trong các công thức này không duy nhất. Bảng sau cho
các hệ số tính theo Cash và Karp:
i Ai Bi,j Ci Di
37 2825
1 ∗ ∗ ∗ ∗ ∗ ∗
378 27648
1 1
2 ∗ ∗ ∗ ∗ 0 0
5 5
3 3 9 250 18575
3 ∗ ∗ ∗
10 40 40 621 48384
3 3 9 6 125 13525
4 − ∗ ∗
5 10 10 5 594 55296
368
11 5 70 35 277
5 1 − − ∗ 0
54 2 27 27 14336
7 1631 175 575 44275 253 512 1
6
8 55296 512 13824 110592 4096 1771 4
Sai số sẽ là:
6
E(h) = y5(x+h) ‐ y4(x+h) = ∑ (Ci − Di )K i (5)
i =1
Chú ý là E(h) là một vec tơ, thành phần Ei(h) biểu diễn sai số của biến yi. Sai
số e(h) ta cần kiểm soát là:
e(h) = max E(h) (6)
i
Ta cũng có thể kiểm soát sai số trung bình bình phương:
∑ E (h)
n
1
e(h) = 2
i (7)
n
i =1
với n là số phương trình bậc 1.
Việc kiểm soát sai số đạt được bằng cách thay đổi h sao cho sai số tại mỗi
bước tính ậiphỉ cỡ sai số mong muốn ε. Sai số khi thực hiên thuật táon Runge
‐ Kutta bậc bốn là O(h5) nên:
5
e(h1 ) ⎛ h1 ⎞
≈⎜ ⎟ (8)
e(h 2 ) ⎝ h 2 ⎠
Giả sử là ta đã tính nghiệm tại bước tính với h1 và có sai số là e(h1). Tại bước
tính với h2 ta muốn có e(h2) = ε thì:
1/ 5
⎡ ε ⎤
h 2 = h1 ⎢ (9)
⎣ e(h1 ) ⎥
⎦
Để dự phòng, ta lấy:
1/ 5
⎡ ε ⎤
h 2 = 0.9h1 ⎢ (10)
⎣ e(h1 ) ⎥
⎦
Ta xây dựng hàm adaptrk() để thực hiện thuật toán này:
function [xsol, ysol] = adaptrk(f, xo, x1, y, n)
% Tich phan Runge‐Kutta bac 5 dung giap phuong trinh y’ = f(x, y) hay y’ =
%f(x).
% xo, x1 ‐ doan tim nghiem.
% y gia tri dau, n dung tim h ban dau
369
h = (x1 ‐ xo)/n;
if size(y, 1) > 1 ;
y = yʹ; % y phai la vec to hang
end
eTol = 1.0e‐9;
n = length(y);
A = [0 1/5 3/10 3/5 1 7/8];
B = [ 0 0 0 0 0
1/5 0 0 0 0
3/40 9/40 0 0 0
3/10 ‐9/10 6/5 0 0
‐11/54 5/2 ‐70/27 35/27 0
1631/55296 175/512 575/13824 44275/110592 253/4096];
C = [37/378 0 250/621 125/594 0 512/1771];
D = [2825/27648 0 18575/48384 13525/55296 277/14336 1/4];
%nghiem ban dau
xsol = zeros(2, 1);
ysol = zeros(2, n);
xsol(1) = xo;
ysol(1,:) = y;
stopper = 0;
k = 1;
for p = 2:5000
% Tinh K tu (2)
K = zeros(6, n);
if nargin(f) > 1
K(1, :) = h*feval(f, xo, y);
else
K(1, :) = h*feval(f, xo);
end
for i = 2:6
BK = zeros(1, n);
for j = 1:i‐1
BK = BK + B(i, j)*K(j, :);
end
if nargin(f) > 1
370
K(i, :) = h*feval(f, xo + A(i)*h, y + BK);
else
K(i, :) = h*feval(f, xo + A(i)*h);
end
end
% tinh su thay doi cua y theo (3) & (4)
dy = zeros(1, n);
E = zeros(1, n);
for i = 1:6
dy = dy + C(i)*K(i,:);
E = E + (C(i) ‐ D(i))*K(i,:);
end
e = sqrt(sum(E.*E)/n);
% neu sai so dat den gia tri cho phep, chap nhan ket qua
% kiem tr dieu kien ket thuc
if e 0) == (xo + hnext >= x1 )
hnext = x1 ‐ xo;
stopper = 1;
end
h = hnext;
371
end
Để tìm nghiệm của phương trình vi phân ta dùng chương trình ctadaptrk.m:
clear all, clc
a = 0;
b = 1;
y = inline(ʹx + yʹ);
ya = 0.5;
n = 10;%so lan tinh chi n = 10
%y = @f4;
[u, v] = adaptrk(y, a, b, ya, n)
plot(u, v)
§6. PHƯƠNG PHÁP BURLIRSCH ‐ STÖR
1. Phương pháp điểm giữa: Công thức điểm giữa của tích phân số của
y′ = f(x,y) là:
y(x + h) = y(x − h) + 2hf [ x,y(x)] (1)
Đây là công thức bậc 2, giống như công thức y’(x)
Euler. Ta xem xét phương pháp này vì đây là cơ
sở của phương pháp Burlisch ‐ Stör dùng tìm
nghiệm có độ chính xác cao. Hình bên minh hoạ
công thức điểm giữa đối với phương trình đơn
x
dạng y′ = f(x,y) . Sự thay đổi y trên hai phía được
xác định bằng: x‐h h x+h
x+ h
y(x + h) − y(x − h) = ∫ y′(x)dx
x−h
và bằng diện tích bên dưới đường cong. Xấp xỉ điểm giữa của diện tích này là
diện tích của hình chữ nhật có gạch chéo.
Bây giờ ta xét ưu điểm của phương pháp điểm giữa khi tìm nghiệm của
phương trình y′ = f(x,y) từ x = x0 đến x = x0 + H với công thức điểm giữa.
Ta chia đoạn tích phân thành n đoạn nhỏ có độ dài mỗi đoạn là h = H / n như
hình bên và tính:
y1 = y 0 + hf0
y 2 = y 0 + 2hf1
372
H
y 3 = y1 + 2hf2 (2)
h
M
y n = y n−2 + 2hfn −1 xo x1 x2 x3 xn‐1 xn
Ta đã dùng khái niệm yi = y(xi) và fi = f(xi, yi). Phương trình đầu tiên trong (2)
dùng công thức Euler để thay cho phương pháp điểm giữa. Các phương trình
khác là các công thức điểm giữa. Kết quả cuối cùng là trung bình cộng của yn
trong (2) và ta có:
y(xo + H) = 0.5 [ y n + (y n −1 + hfn )] (3)
2. Ngoại suy Richardson: Ta có thể thấy sai số trong (3) là:
E = c 1h 2 + c 2 h 4 + c 3 h 6 + L
Để giảm bớt sai số ta dùng phương pháp ngoại suy Richardson. Cụ thể ta tính
y(xo+H) với một giá trị nào đó của h và rồi lặp lại quá trình tính với h/2. Gọi
kết quả là g(h) và g(h1) ta có ngoại suy Richardson:
4g(h1 ) − g(h)
y(x o + H) =
3
Ta xây dựng hàm midpoint() để kết hợp phương pháp điểm giữa và phương
pháp ngoại suy Richardson. Đầu tiên phương pháp điểm giữa được dùng cho
2 tích phân. Số bước tính được tăng gấp đôi trong các lần lặp sau, mỗi lần lặp
đều dùng ngoại suy Richardson. Chương trình dừng khi sai số nhỏ hơn sai số
cho phép.
function y = midpoint(f, x, x1, y, tol)
% Phuong phap diem giua dung cho phuong trinh yʹ = f(x,y) hay y’ = f(x).
if size(y, 1) > 1 ;
y = yʹ;
end % y phai la vec to hang
if nargin for k = 2:kmax
% Tang gap doi so buoc va tinh chinh ket qua
% ngoai suy Richardson
nsteps = 2*k;
r(k, 1:n) = mid(f, x, x1, y, nsteps);
r = richardson(r, k, n);
% kiem tra hoi tu.
dr = r(1, 1:n) ‐ rold;
e = sqrt(dot(dr, dr)/n);
if e 1
y1 = y0 + h*feval(f, x, y0);
else
y1 = y0 + h*feval(f, x);
end
for i = 1:nsteps‐1
x = x + h;
if nargin(f) > 1
y2 = y0 + 2.0*h*feval(f, x, y1);
else
y2 = y0 + 2.0*h*feval(f, x);
end
y0 = y1;
y1 = y2;
end
if nargin(f) > 1
374
y = 0.5*(y1 + y0 + h*feval(f, x, y2));
else
y = 0.5*(y1 + y0 + h*feval(f, x));
end
function r = richardson(r, k, n)
% Richardson extrapolation.
for j = k‐1:‐1:1
c =(k/(k‐1))^(2*(k‐j));
r(j, 1:n) =(c*r(j+1, 1:n) ‐ r(j, 1:n))/(c ‐ 1.0);
end
3. Thuật toán Burlisch ‐ Stör: Phương pháp điểm giữa có nhược điểm là
nghiệm nằm tại điểm giữa của khoảng tìm nghiệm không được tinh chỉnh
bằng phương pháp ngoại suy Richardson. Khuyết điểm này được khác phục
trong phương pháp Burlisch ‐ Stör. Ý tưởng của phương pháp này là áp dụng
phương pháp điểm giữa trên từng đoạn. Ta xây dựng hàm burlischstoer() để
thực hiện thuật toán này:
function [xout, yout] = burlischstoer(f, x, x1, y, H, tol)
% Phuong phap Bulirsch‐Stoer giai phuong trinh yʹ = F(x, y) hay y’ = f(x).
%[x, x1] la khoang tim nghiem.
% H = do tang sau moi lan tinh
if size(y, 1) > 1
y = yʹ;
end % y phai la vec to hang
if nargin k = k + 1;
H = min(H, x1 ‐ x);
y = midpoint(f, x, x + H, y, tol);
x = x + H;
xout(k) = x;
yout(k, :) = y;
end
Để giải phương trình ta dùng chương trình ctburlischstoer.m:
clear all, clc
a = 0;
b = 1;
y = @f3;
ya = 1;
H = .1;
[u, v] = burlischstoer(y, a, b, ya, H)
plot(u, v)
§7. PHƯƠNG PHÁP CHUỖI TAYLOR
Phương pháp chuỗi Taylor đơn giản về ý tưởng và có độ chính xác cao.
Cơ sở của phương pháp này là cắt chuỗi Taylor của y theo x:
1 1 1 (m)
y(x + h) ≈ y(x) + y′(x)h + y′′(x)h 2 + y′′(x)h 3 + L + y (x)h m (1)
2! 3! m!
Do phương trình (1) dự đoán trước y tại (x + h) từ các thông tin có tại x, nó
cũng là công thức tích phân. Số hạng cuối trong (1) là bậc của tích phân. Như
vậy (1) là tích phân bậc m. Sai số là:
1
E= y(m +1) (ξ)h m +1 x < ξ < x + h
(m + 1)!
Dùng xấp xỉ đạo hàm:
y(m) (x + h) − y(m) (x)
y(m +1) (ξ) ≈
h
ta có:
h m ⎡ (m)
E= y (x + h) − y(m) (x) ⎤ (2)
(m + 1)! ⎣⎢ ⎥
⎦
Ta xây dựng hàm taylor() để giải bài toán trên:
376
function [xout, yout] = taylor(deriv, x, y, x1, h)
% Tich phan chuoi Taylor bac 4.
% x, y = cac gia tri dau; i la vec to hang.
% x1 = gia tri cuoi cua x
if size(y,1) > 1;
y = yʹ;
end
xout = zeros(2, 1);
yout = zeros(2, length(y));
xout(1) = x;
yout(1, :) = y;
k = 1;
while x plot(x, y);
§8. PHƯƠNG PHÁP DỰ ĐOÁN ‐ HIỆU CHỈNH
1. Phương pháp Adam ‐ Bashfort ‐ Moulton: Năm 1855, nhà toán học người
Anh Adams đề xuất một phương pháp đa bước giải bài toán Cauchy theo yêu
cầu của ông Bashforth, một chuyên gia kỹ thuật pháo binh Anh. Kết quả của
Adams sau này bị quên lãng. Mãi đến đầu thế kỷ 20, nhà toán học Nauy khi
tính quỹ đạo của hạt điện tích rời xa mặt trời với vận tốc lớn đã phát minh lại
công thức Adams. Sau này viện sỹ Krylov đã hoàn thiện phương pháp
Adams. Phương pháp Adams ‐ Bashfort ‐ Moulton (ABM) gồm hai bước.
Bước dầu tiên là xấp xỉ f(x, y) bằng một đa thức(ví dụ đa thức Lagrange) bậc 4
qua 4 điểm:
{( t k−3 ,fk−3 ) , ( t k−2 ,fk−2 ) , ( t k−1 ,fk−1 ) , ( t k ,fk )}
và thay thế đa thức này vào phương trình vi phân để có được giá trị dự báo
yk+1:
( )
h
h
pk +1 = y k + ∫ l 3 (t)dt = y k + −9fk −3 + 37fk −2 − 59fk −1 + 55fk (1a)
0
24
Bước thứ hai là lặp lại công việc với 4 điểm được cập nhật:
{( t k −2 ,fk−2 ) , ( t k −1 ,fk−1 ) , ( t k ,fk ) , ( t k+1 ,fk+1 )}
và nhận giá trị hiệu chỉnh:
fk+1 = f(t k+1 ,pk+1 )
( )
h
h
c k = y k + ∫ l′ (t)dt = y k +
3 fk−2 − 5fk −1 + 19fk + 9fk +1 (1b)
0
24
Ta viết khai triển Taylor của yk+1 lân cận tk và của yk lân cận tk+1:
h2 h3
′ ′′
y k+1 = y k + hfk + fk + fk + L (2a)
2 3!
h2 h3
′ ′′
y k = y k+1 − hfk+1 + fk+1 − fk+1 + L
2 3!
h2 h3
′ ′′
y k+1 = y k + hfk+1 − fk+1 + fk+1 + L (2b)
2 3!
và thay thế các đạo hàm bậc 1, 2, 3 bằng các xấp xỉ
⎛ 1 3 11 ⎞
2 − fk − 3 + fk−2 − 3fk−1 + fk
h ⎜ 3 2 6 + h f + L⎟
1 3 (4)
y k+1 = y k + hfk + ⎜ k ⎟
2 ⎜ h 4 ⎟
⎝ ⎠
378
h 3 ⎛ −fk−3 + 4fk−2 − 5fk−1 + 2fk 11 2 (4)
+ ⎜ + h fk + L ⎞ ⎟
3! ⎝ h 2
12 ⎠
h ⎛ −fk−3 + 3fk−2 − 3fk−1 + fk 3 (4)
4
⎞ h (4)
5
+ ⎜ + hfk + L ⎟ + fk + L
4! ⎝ h3 2 ⎠ 120
= y k + ( h
24
−9fk−3 + 37fk−2 − 59fk−1 + 55fk + ) 251 5 (4)
720
h fk + L
251 5 (4)
≈ pk+1 + h fk (3a)
720
⎛ − 1 f + 3 f − 3f + 11 f ⎞
h 2 ⎜ 3 k − 2 2 k −1 k
6
k +1 1 3 (4) ⎟
y k+1 = y k + hfk+1 − ⎜ + h fk+1 + L ⎟
2 ⎜ h 4 ⎟
⎝ ⎠
h ⎛ −f + 4fk−1 − 5fk + 2fk+1 11 2 (4)
3
⎞
+ ⎜ k −2 + h fk+1 + L ⎟
3! ⎝ h 2
12 ⎠
h ⎛ −fk−2 + 3fk−1 − 3fk + fk+1 3 (4)
4
⎞ + h f (4) + L
5
− ⎜ + hfk+1 + L ⎟ k +1
4! ⎝ h3 2 ⎠ 120
= y k + ( h
24 )
fk−2 − 5fk−1 + 19fk + 9fk+1 −
19 5 (4)
720
h fk+1 + L
19 5 (4)
≈ c k +1 − h fk +1 (3b)
720
Từ các phương trình này và giả sử rằng fk +1 ≅ fk ≅ K ta có thể viết các sai số
(4) (4)
dự đoán/hiệu chỉnh:
251 5 (4) 251 5
EP ,k +1 = y k +1 − pk+1 ≈ h fk ≅ Kh (4a)
720 720
19 5 (4) 19
EC ,k +1 = y k +1 − c k+1 ≈ − h fk+1 ≅ − Kh 5 (4b)
720 720
Do K chưa biết nên ta phải tìm nó. Ta có;
270 270 270
EP ,k +1 − EC ,k +1 = c k +1 − pk+1 ≅ Kh 5 ≡ EP ,k+1 ≡ − EC ,k+1 (5)
720 251 19
Do vậy ta có các công thức dùng để đánh giá sai số:
251
EP ,k +1 = y k +1 − pk+1 ≅ ( c k+1 − pk+1 ) (6a)
720
19
EC ,k +1 = y k +1 − c k+1 ≅ − ( c k +1 − p k +1 ) (6b)
270
Tóm lại, thuật toán Adams ‐ Bashforth ‐ Moulton gồm:
Dự đoán: p k +1 = y k + ( h
24 )
−9fk−3 + 37fk−2 − 59fk−1 + 55fk (7a)
379
