CHƯƠNG 5: CÁC PHƯƠNG TRÌNH PHI TUYẾN
Nếu phương trình đại số hay siêu việt khá phức tạp thì ít khi tìm được nghiệm đúng. Bởi vậy việc tìm nghiệm gần đúng và ước lượng sai số là rất cần thiết. Ta xét phương trình : f(x) = 0 (1) với f(x) là hàm cho trước của biến x. Chúng ta cần tìm giá trị gần đúng của nghiệm của phương trình này. Quá trình giải thường chia làm hai bước: bước sơ bộ và bước kiện toàn nghiệm. Bước giải sơ bộ có 3 nhiệm vụ: vây nghiệm, tách nghiệm và thu hẹp khoảng chứa nghiệm. Vây nghiệm là tìm xem các nghiệm của phương trình có thể nằm trên ...
CHƯƠNG 5: CÁC PHƯƠNG TRÌNH PHI TUYẾN
§1. KHÁI NIỆM CHUNG
Nếu phương trình đại số hay siêu việt khá phức tạp thì ít khi tìm được
nghiệm đúng. Bởi vậy việc tìm nghiệm gần đúng và ước lượng sai số là rất
cần thiết.
Ta xét phương trình :
f(x) = 0 (1)
với f(x) là hàm cho trước của biến x. Chúng ta cần tìm giá trị gần đúng của
nghiệm của phương trình này.
Quá trình giải thường chia làm hai bước: bước sơ bộ và bước kiện toàn
nghiệm.
Bước giải sơ bộ có 3 nhiệm vụ: vây nghiệm, tách nghiệm và thu hẹp
khoảng chứa nghiệm.
Vây nghiệm là tìm xem các nghiệm của phương trình có thể nằm trên
những đoạn nào của trục x. Tách nghiệm là tìm các khoảng chứa nghiệm sao
cho trong mỗi khoảng chỉ có đúng một nghiệm. Thu hẹp khoảng chứa
nghiệm là làm cho khoảng chứa nghiệm càng nhỏ càng tốt. Sau bước sơ bộ ta
có khoảng chứa nghiệm đủ nhỏ. Để xác định khoảng chứa nghiệm ta có thể
dùng phương pháp đồ thị. Ngoài ra ta cũng có thể tìm nghiệm bằng phương
pháp tìm tăng dần. Ý tưởng của phương pháp này là nếy f1(x).f2(x) x1 = a; f1 = feval(func,x1);
x2 = a + dx; f2 = feval(func,x2);
while f1*f2 > 0.0
if x1 >= b
x1 = NaN; x2 = NaN;
return
end
x1 = x2; f1 = f2;
x2 = x1 + dx; f2 = feval(func,x2);
end
Khi phát hiện thấy khoảng chứa nghiệm, hàm trả về giá trị biên của đoạn.
Nếu không có nghiệm, x1 = x2 = NaN. Ta gọi rootsearch() nhiều lần để phát
hiện hết các đoạn chứa nghiệm. Với ví dụ tìm khoảng chứa nghiệm của hàm
f(x) = x3 ‐ 10x2 + 5 ta dùng chương trình ctrootsearch.m
clear all, clc
f = inline(ʹx^3 ‐ 10*x^2 + 5ʹ);
[x1, x2] = rootsearch(f,2,10,.2)
Bước kiện toàn nghiệm tìm các nghiệm gần đúng theo yêu cầu đặt ra.
Có rất nhiều phương pháp xác định nghiệm của (1). Sau đây chúng ta
xét từng phương pháp.
§2. PHƯƠNG PHÁP LẶP ĐƠN
Giả sử phương trình (1) được đưa về dạng tương đương:
x = g(x) (2)
từ giá trị xo nào đó gọi là giá trị lặp đầu tiên ta lập dãy xấp xỉ bằng công thức:
xn = g(xn‐1) (3)
với n = 1,2,....
Hàm g(x) được gọi là hàm lặp. Nếu dãy xn → α khi n →∝ thì ta nói phép lặp
(3) hội tụ.
Ta có định lí: Xét phương pháp lặp (3), giả sử:
‐ [a, b] là khoảng chứa nghiệm α của phương trình (1) tức là của (2)
‐ mọi xn tính theo (3) đều thuộc [a, b]
‐ g(x) có đạo hàm thoả mãn :
242
g′(x) ≤ q < 1 a < x < b (4)
trong đó q là một hằng số thì phương pháp lặp (3) hội tụ
Ta có thể minh hoạ phép lặp trên bằng hình vẽ sau.
x1 xo xo x1
Ta xây dựng hàm simpiter() để lặp
function [x, err, xx] = simpiter(g, x0, tolx, maxiter)
% giai pt x = g(x) tu x0 bang cah lap
%vao : g, x0 = ham va gia tri dau
% tolx = sai so mong muon
% maxiter = so lan lap max
%ra: x = nghiem
% err = sai so |x(k) ‐ x(k ‐ 1)|
% xx = cac gia tri trung gian
if nargin end
x = xx(k);
if k == maxiter
fprintf(ʹKhong hoi tu sau %d lan lap\nʹ, maxiter)
else
fprintf(ʹHoi tu sau %d lan lap\nʹ,k)
end
Để tính lại ví dụ trên ta dùng chương trình ctsimpiter4_2.m
clear all, clc
f = inline(ʹ‐0.5*((x ‐ 1).^2 ‐ 3)ʹ);
[x, ss, xx] = simpiter(f, 0.5,.00001,200)
§3. PHƯƠNG PHÁP CHIA ĐÔI CUNG
Giả sử cho phương trình f(x) = 0 với
f(x) liên tục trên đoạn [a, b] và f(a).f(b) tol = eps;
fa = feval(f, a);
fb = feval(f, b);
if fa*fb > 0
error(ʹNghiem khong o trong doan nayʹ);
end
for k = 1: maxiter
xx(k) = (a + b)/2;
fx = feval(f, xx(k));
err = (b ‐ a)/2;
if abs(fx) Trong đó y
= −f(a)
1 h −f(a)+ f(b) (b−a)
Tiếp theo dùng cách đó với đoạn [ a, x1]
hay [x1, b] mà hai đầu hàm nhận giá trị trái a x1 ξ x
dấu ta được nghiệm gần đúng x2 v.v. b
Về mặt hình học, phương pháp này có
nghĩa là kẻ dây cung của đường cong f(x)
qua hai điểm A[a, f(a)] và B[b, f(b)] hay nói cách khác là tuyến tính hoá hàm
f(x) trong đoạn [a, b].
Thật vậy phương trình dây cung AB có dạng:
f(a) − f(b) af(b) − bf(a)
y = x+
a−b a−b
Cho x = x1, y = 0 ta có
af(b) − bf(a)
x1 = (1)
f(b) − f(a)
Ta xây dựng hàm chord() để thực hiện thuật toán trên
function [x, err, xx] = chord(f, a, b, tolx, maxiter)
%giai pt f(x) = 0 bg phuong phap day cung.
%vao : f ‐ ham can tim nghiem
% a/b ‐ khoang tim nghiem
% tolx ‐ sai so mong muon cua nghiem
% maxiter lan lap max
%ra: x ‐ nghiem
% err ‐ sai so
% xx ‐ cac gia tri trung gian
tolfun = eps;
fa = feval(f, a);
fb = feval(f, b);
if fa*fb > 0
error(ʹNghiem khong o trong doan nay !ʹ);
end
for k = 1: maxiter
xx(k) = (a*fb ‐ b*fa)/(fb ‐ fa); %pt.(1)
fx = feval(f, xx(k));
err = min(abs(xx(k) ‐ a), abs(b ‐ xx(k)));
246
if abs(fx) f′′(x i ) 2
e i +1 = − ei y
2f′(xi )
Minh hoạ hình học của thuật toán
Newton ‐ Raphson như hình bên.
Thuật toán được tóm lược như sau:
a x
‐ cho xo
x1 b = xo
f(x)
‐ tính ∆x = −
f ʹ(x)
‐ cho x = x + ∆x
‐ lặp lại bước 2 và 3 cho đến khi |∆x| ≤ ε
Ta xây dựng hàm newtonraphson() để thực hiện thuật toán trên.
function [x, fx, xx] = newtonraphson(f, df, x0, tolx, maxiter)
%giai pt f(x) = 0 bang pp Newton‐Raphson.
%vao: f = ftn to be given as a string ’f’ if defined in an M‐file
% df = df(x)/dx (neu khong cho se dung dao ham so.)
% x0 = gia tri ban dau
% tolx = sai so mong muon
% maxiter = so lan lap max
%ra: x = nghiem
% fx = f(x(last)), xx = cac gia tri trung gian
h = 1e‐4;
h2 = 2*h;
tolf = eps;
if nargin == 4 & isnumeric(df)
maxiter = tolx;
tolx = x0;
x0 = df;
end
xx(1) = x0;
fx = feval(f,x0);
for k = 1: maxiter
if ~isnumeric(df)
dfdx = feval(df, xx(k)); %dao ham cua ham
else
dfdx = (feval(f, xx(k) + h)‐feval(f, xx(k) ‐ h))/h2; %dao ham so
248
end
dx = ‐fx/dfdx;
xx(k+1) = xx(k) + dx; %pt.(3)
fx = feval(f, xx(k + 1));
if abs(fx)Phương trình (2) chính là biểu thức tổng quát của phép lặp. Hai giá trị đầu
tiên x1 và x2 cần để khởi động phép lặp. Quá trình lặp được minh hoạ bằng
hình a
y y
f(x)
x2 x1 x0 x1 x0
f(x) x x
a b
Phép lặp có thể không hội tụ (hình b). Tuy nhiên khi hội tụ, nó hội rất nhanh.
Ta xây dựng hàm secant() để thực hiện thuật toán trên.
function [x, fx, xx] = secant(f, x0, x1, tolx, maxiter)
% giai pt f(x) = 0 bg phuong phap day cung
%vao : f ‐ ham can tim nghiem
% x0, x1 ‐ gia tri khoi dong phep lap
% tolx ‐ sai so mong muon
% maxiter ‐ so lan lap max
%ra: x = nghiem
% fx = f(x(last)), xx = cac gia tri trung gian
h = (x1 ‐ x0)/100;
h2 = 2*h;
tolfun = eps;
xx(1) = x0;
fx = feval(f, x0);
for k = 1: maxiter
if k end
dx = ‐fx/dfdx;
xx(k + 1) = xx(k) + dx; %pt.(2)
fx0 = fx;
fx = feval(f, xx(k+1));
if abs(fx) Trong quá trình này ta tính được f(x1), f(x2) và f(x3). Qua 3 điểm này ta có một
đường cong bậc 2 và tìm được nghiệm x của đường cong bậc 2 này. Nếu x
nằm trong đoạn [x1, x2] như hình trên thì giá trị này được chấp nhận. Tiếp
theo ta tìm nghiệm trong đoạn [x1, x3] hay [x3, x2] tuỳ theo vị trí của x.
Đa thức nội suy Lagrange qua 3 điểm là:
(f − f2 )(f − f3 ) (f − f1 )(f − f3 ) (f − f1 )(f − f2 )
x(y) = x1+ x2+ x3
(f1 − f2 )(f1 − f3 ) (f2 − f1 )(f2 − f3 ) (f3 − f1 )(f3 − f2 )
Cho y = 0 ta có:
f f x (f − f ) + f1f3 x 2 (f3 − f1 ) + f1f2 x 3 (f1 − f2 )
x=− 2 3 1 2 3 (1)
(f1 − f2 )(f2 − f3 )(f3 − f1 )
Độ thay đổi của nghiệm là:
x (f − f )(f − f + f ) + f2 x1 (f2 − f3 ) + f1x 2 (f3 − f1 )
∆x = x − x 3 = f3 3 1 2 2 3 1 (2)
(f1 − f2 )(f2 − f3 )(f3 − f1 )
Ta xây dựng hàm brent() để thực hiện thuật toán trên
function root = brent(f, a, b, tol)
% giai pt f(x) = 0 bang thuat toan Brent.
% Cu phap: root = brent(f, a, b, tol)
% vao: f = ham can tim nghiem
% a, b = doan chua nghiem
% tol = sai so cho truoc
if nargin end
if f1*f2 > 0.0
error(ʹNghiem khong nam trong doan nayʹ)
end
x3 = 0.5*(a + b);
% bat dau lap.
for i = 1:30
f3 = feval(f, x3);
if abs(f3) end
% cho x = x3 & chon x1, x2 moi sao cho x1 . . . . . . . . . . . . . . . . . . .
| xn+1 ‐ xn | ≤ q | xn ‐ xn‐1 |
Điều này có nghĩa là dãy xi+1 ‐ xi , một cách gần đúng, là một cấp số nhân. Ta
cũng coi rằng dãy xn ‐ y với y là nghiệm đúng của (1), gần đúng như một cấp
số nhân có công sai q . Như vậy :
x n +1 − y
= q < 1 (7)
xn − y
hay : x n +1 − y = q(x n − y) (8)
Tương tự ta có :
x n + 2 − y = q(x n +1 − y) (9)
Từ (8) và (9) ta có :
x − x n +1
q = n+2 (10)
x n +1 − x n
Thay giá trị của q vừa tính ở (10) vào biểu thức của q ở trên ta có :
( x n − x n +1 )
2
y = xn − (11)
x n − 2x n +1 + x n + 2
Công thức (11) được gọi là công thức ngoại suy Adam. Như vậy theo (11)
trước hết ta dùng phương pháp lặp để tính giá trị gần đúng xn+2, xn+1, xn của
nghiệm và sau đó theo (11) tìm được nghiệm với sai số nhỏ hơn.
Khi phương pháp lặp được kết hợp với phương pháp Aitken ta có
phương pháp Steffensen. Bắt đầu lặp từ x0, hai bước lặp được dùng để tính:
x1 = f(x0) x2 = f(x1)
và sau đó dùng thuật toán Aitken để tinh y0 theo (11). Để tiếp tục lặp ta cho
x0=y0 và lặp lại bước tính trước.
Ta xây dựng hàm aitstef() để thực hiện hai thuật toán trên
function [x, y] = aitstef(g, x0, tolx, maxiter)
% phuong phap Aitken ‐ Steffenson
% giai pt x = g(x)
xstart = x0;
f0 = 0;
f0old = 1.;
count = 0;
while ((count .0000001))
count = count + 1;
f0old = f0;
255
x1 = feval(g, x0);
x2 = feval(g, x1);
f0 = x0 ‐ (x1 ‐ x0)*(x1 ‐ x0) / (x2 ‐ 2.*x1 + x0);
x0 = x1;
end
x = f0;
fprintf(ʹ\nʹ);
fprintf(ʹPhuong phap Aitken‐Steffensonʹ);
x0 = xstart;
count = 0;
f0 = 0;
x2 = 1.;
while ((count .0000001))
count = count+1;
x1 = feval(g, x0);
x2 = feval(g, x1);
f0 = x0 ‐ (x1 ‐ x0)*(x1 ‐ x0) / (x2 ‐ 2.*x1 + x0);
x0 = f0;
end
y = f0;
Để tìm nghiệm của phương trình x = (2 ‐ ex + x2)/3 = g(x) ta dùng chương trình
ctaitstef.m
clear all, clc
f = inline(ʹ(2.‐exp(x)+x.^2)/3ʹ);
[x, y] = aitstef(f,0., 1e‐4, 50)
§9. PHƯƠNG PHÁP MUELLER
Trong phương pháp dây cung khi tìm nghiệm trong đoạn [a, b] ta xấp
xỉ hàm bằng một đường thẳng. Tuy nhiên để giảm lượng tính toán và để
nghiệm hội tụ nhanh hơn ta có thể dùng phương pháp Muller. Nội dung của
phương pháp này là thay hàm trong đoạn [a, b] bằng một đường cong bậc 2
mà ta hoàn toàn có thể tìm nghiệm chính xác của nó.
256
Gọi các điểm đó có hoành độ lần lượt là a = x2, b = x1 và ta chọn thêm
một điểm x0 nằm trong đoạn [x2, x1]. Gọi
h1 = x1 ‐ x0 x0, f0
h2 = x0 ‐ x2
v = x ‐ x0 x1, f1
f(x0) = f0
f(x1) = f1
f(x2) = f2 f(x) x2, f2 h2 h1
h2
γ=
h1
Qua 3 điểm này ta có một đường parabol:
y = av2 + bv + c
Ta tìm các hệ số a, b, c từ các giá trị đã biết v:
v = 0 (x = x 0 ) a(0)2 + b(0) + c = f0
v = h1 (x = x1 ) ah1 + bh1 + c = f1
2
v = − h 2 (x = x 2 ) ah 2 − bh 2 + c = f2
2
Từ đó ta có :
γf − f (1 + γ ) + f2
a = 1 02
γh 1 (1 + γ )
f1 − f0 − ah 1
2
b=
h1
c = f0
Sau đó ta tìm nghiệm của phương trình av2 + bv + c = 0 và có :
2c
n = x0 −
b ± b 2 − 4ac
Dấu của mẫu số được chọn sao cho nó có giá trị tuyệt đối lớn nhất, nghĩa là
khi b > 0, ta lấy dấu +, khi b %vao ‐ f la ham can tim nghiem
% ‐ a, b la doan chua nghiem
% ‐ maxiter so lan lap max
%ra ‐ p xap xi Muller cua f
% ‐ y la gia tri y = f(p)
% ‐ err sai so thuc cua nghiem.
%khoi gan a,b,x0 va cac gia tri tuong ung cua ham
x0 = (a + b)/2.;
P = [x0 a b];
Y = feval(f, P);
delta = 1e‐4;
%tinh cac he so cua pt bac hai
for k = 1:maxiter
h0 = P(1) ‐ P(3);
h1 = P(2) ‐ P(3);
e0 = Y(1) ‐ Y(3);
e1 = Y(2) ‐ Y(3);
c = Y(3);
denom = h1*h0^2 ‐ h0*h1^2;
a = (e0*h1 ‐ e1*h0)/denom;
b = (e1*h0^2 ‐ e0*h1^2)/denom;
%neu nghiem phuc
if b^2‐4*a*c > 0
disc = sqrt(b^2 ‐ 4*a*c);
else
disc = 0;
end
%tim nghiem nho nhat
if b P = Q;
Y = feval(f, P);
end
if abs(p‐P(3)) function x = halley1(f, x, maxiter)
%ham dung de tim nghiem cua pt f(x) = 0
%vao: ‐ ham can tim nghiem f
% ‐ gia tri dau x0
% ‐ so lan lap cuc dai maxiter
%ra ‐ nghiem x
% dung dao ham so
i = 0;
h = 1e‐4;
while (i
