Giáo trình Matlab cơ bản - Trần Văn Chính

ải với c ∈ [a, b] ta có : xn+1- xn = f(xn) - f(xn-1) = (xn - xn-1)f’(c) (5) Vì phép lặp (1) nên : | xn+1- xn | ≤ q | xn - xn-1 | (6) Do (6) đúng với mọi n nên cho n = 1 , 2 , 3 , . . . ta có : | x2 - x1 | ≤ q | x1 - xo | | x3 - x2 | ≤ q | x2 - 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ó 254 công sai q . Như vậy : n 1 n x y q 1 x y + − = <− (7) hay : n 1 nx y q(x y)+ − = − (8) Tương tự ta có : n 2 n 1x y q(x y)+ +− = − (9) Từ (8) và (9) ta có : n 2 n 1 n 1 n x xq x x + + + −= − (10) Thay giá trị của q vừa tính ở (10) vào biểu thức của q ở trên ta có : ( )2n n 1n n n 1 n 2 x x y x x 2x x + + + −= − − + (11) 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; 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'); 255 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ó. 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 h2 = x0 - x2 v = x - x0 f(x0) = f0 f(x1) = f1 f(x2) = f2 1 2 h h=γ 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: h1 h2 f(x)  x0, f0 x1, f1 x2, f2 256 2 0 0 2 1 1 1 1 1 2 2 2 2 2 2 v 0 (x x ) a(0) b(0) c f v h (x x ) ah bh c f v h (x x ) ah bh c f = = + + = = = + + = = − = − + = Từ đó ta có : 0 1 2 101 2 1 201 fc h ahffb )1(h f)1(ffa = −−= γ+γ +γ+−γ= Sau đó ta tìm nghiệm của phương trình av2 + bv + c = 0 và có : 0 2 2cn x b b 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 < 0 ta lấy dấu -. Tiếp đó ta chọn x0 làm một trong 3 điểm để tính xấp xỉ mới. Các điểm này được chọn gần nhau nhất, nghĩa là nếu nghiệm n ở bên phải x0 thì ba điểm tính mới là x0, x1 và n; nếu n nằm bên trái x0 thì 3 điểm tính mới là x0, x2 và nghiệm. Tiếp tục quá trình tính đến khi đạt độ chính xác yêu cầu thì dừng lại. Ta xây dựng hàm muller() để thực hiện thuật toán trên function p = muller(f, a, b, maxiter) % giai pt f(x) = 0 %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); 257 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 < 0 disc = -disc; end z = -2*c/(b + disc); p = P(3) + z; %sap xep lai cac tri tinh lap if abs(p - P(2)) < abs(p -P(1)) Q = [P(2) P(1) P(3)]; P = Q; Y = feval(f, P); end if abs(p-P(3)) < abs(p-P(2)) R = [P(1) P(3) P(2)]; P = R; Y = feval(f, P); end %cac tri tinh lan moi P(3) = p; Y(3) = feval(f, P(3)); y = Y(3); %dieu kien dung lap err = abs(z); relerr = err/(abs(p) + delta); if (err < delta)|(relerr < delta)|(abs(y) < eps) break end end Để giải phương trình sin(x) - 0.5*x = 0 ta dùng chương trình ctmuller.m clear all, clc format long f = inline('sin(x) - 0.5*x'); 258 x = muller(f, 1.8, 2.2, 50) §10. PHƯƠNG PHÁP HALLEY Phép lặp Newton tìm nghiệm của hàm phương trình x = g(x) là: ′ f(x)g(x) = x ‐  f (x) (1) Tốc độ hội tụ tăng đáng kể khi hàm có nghiệm đơn. Để tăng tốc độ hội tụ, Edmon Halley đưa ra công thức lặp: [ ] −⎧ ⎫′′⎪ ⎪−⎨ ⎬′ ′⎪ ⎪⎩ ⎭ 1 2 f(x) f(x)f (x)h(x) = x ‐  1 f (x) 2 f (x) (2) Ta xây dựng hàm halley1() để thực hiện thuật toán trên 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 < maxiter) f1 = feval(f, x); df1 = (feval(f, x + h)-feval(f, x - h))/(2*h); ddf1 = (feval(f, x + h)-2*feval(f, x)+feval(f, x - h))/(h*h); hx = x - (f1/df1)*1./(1 - (f1*ddf1)/(2*(df1)^2)) x = hx; i = i + 1; if (abs(f1) < eps) break; end end hay dùng hàm halley2() function x = halley2(f, x, maxiter) %ham dung de tim nghiem cua pt f(x) = 0 %vao: - ham can tim nghiem f % - gia tri dau x % - so lan lap cuc dai maxiter %ra - nghiem x % dung dao ham symbolic 259 df = diff(f, x); ddf = diff(f, x, 2); i = 0; while (i < maxiter) f1 = eval(f); 262 df1 = eval(df); ddf1 = eval(ddf); hx = x - (f1/df1)*1./(1 - (f1*ddf1)/(2*(df1)^2)); x = hx; i = i + 1; if (abs(f1) < eps) break; end end Để giải phương trình f(x) = x3 - 3x + 2 = 0 ta dùng chương trình cthalley.m: clc, clear all %f = inline('x.^3 - 3*x + 2');%khi dung halley1() %x = halley1(f, -3, 50); syms x f = x^3 - 3*x + 2;%khi dung halley2() x = halley2(f, -3, 50) §11. PHƯƠNG PHÁP CHEBYSHEV Khi tìm nghiệm của phương trình đại số tuyến tính hay phương trình siêu việt f(x) = 0 ta có thể dùng một hàm có 4 thông số để xấp xỉ hàm f(x) y(x) = p1 + p2(x - p3)e (1) Các thông số p1 và p3 tạo sự chuyển dịch theo các trục; thông số p xác định biên độ và e cung cấp độ cong của hàm. Ta khảo sát hàm f(x) trên đoạn [a, b] trong đó f(a).f(b) < 0, nghĩa là trong đoạn [a, b] tồn tại nghiệm của phương trình f(x) = 0. Ta có thêm điều kiện f'(x).f''(x) ≠ 0 ∀x ∈ [a, b]. Gọi xi ∈ [a, b] là lần xấp xỉ thứ i của nghiệm thì nghiệm lần thứ i + 1 theo công thức Popovski là: 1 e i 1 i 2 f e f.fx x (e 1) 1 1 f e 1 f+ ⎡ ⎤′ ′′⎛ ⎞⎢ ⎥− = − − −⎜ ⎟′′ ′−⎢ ⎥⎝ ⎠⎣ ⎦ (2) Khi e = -1 i 1 i 2 f.fx x 0.5f.f f+ ′− = ′′ ′− và đó là phép lặp Halley Khi e → 1: i 1 i fx x f+ − = − ′ và đó là phép lặp Newton Khi e = 0.5 263 22 i i i i 1 i 3 3 i i 1 0.5f.ff f(x ) f (x ) f (x )fx x f f (x ) 2f (x )+ ′′+⎛ ⎞− ⎜ ⎟ ′′×′⎝ ⎠− = = − −′ ′ ′ và ta có phép lặp Chebyshev. Ta xây dựmg hàm chebyiter() để thực hiện thuật toán trên function [x, fx, xx] = chebyiter(f, x0, tol, maxiter) %giai pt f(x) = 0 bang pp Chebyshev. %vao: f = ham can tim nghiem % x0 = gia tri ban dau % tol = sai so mong muon % maxiter = so lan lap max %ra: x = nghiem % fx , xx = cac gia tri trung gian if nargin < 4 maxiter = 200; end if nargin < 3 maxiter = 100; tol = 1e-4; end h = 1e-4; h2 = 2*h; xx(1) = x0; fx = feval(f, x0); for k = 1:maxiter df = (feval(f, xx(k) + h) - feval(f, xx(k) - h))/h2; %dao ham so d2f = (feval(f, xx(k) + h) - 2*feval(f, xx(k)) + feval(f, xx(k) - h))/h^2; dx = - fx/df^3 - 0.5*fx^2*d2f/df^3; xx(k+1) = xx(k) + dx; fx = feval(f, xx(k+1)); if abs(fx) < tol | abs(dx) < tol break; end end x = xx(k + 1); Để giải phương trình ta dùng chương trình ctchebyiter clear all, clc f = inline('x.^3 - 10*x.^2 + 5'); x = chebyiter(f, -3, 1e-4) 264 §12. PHƯƠNG PHÁP NEWTON DÙNG CHO HỆ PHI TUYẾN Phương pháp Newton có thể được tổng quát hoá để giải hệ phương trình phi tuyến dạng : 1 1 2 3 n 2 1 2 3 n n 1 2 3 n f (x ,x ,x ,...,x ) 0 f (x ,x ,x ,...,x ) 0 f (x ,x ,x ,...,x ) 0 =⎧⎪ =⎪⎨⋅ ⋅ ⋅ ⋅ ⋅⎪⎪ =⎩ hay viết gọn hơn dưới dạng : F(X) = 0 Trong đó : X = (x1, x2, x3,....., xn) Với một phương trình một biến, công thức Newton là : )x(f )x(fxx i i i1i ′−=+ hay : f'(xi).Δx = -f(xi) với Δx = xi+1 - xi Đối với hệ phương trình, công thức lặp là : J(Xi)ΔX = -F(Xi) Trong đó J(Xi) là toán tử Jacobi. Nó là một ma trận bậc n ( n - tương ứng với số thành phần trong vectơ X) có dạng : ⎟⎟ ⎟⎟ ⎟⎟ ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎜ ⎝ ⎛ ∂ ∂⋅⋅⋅∂ ∂ ∂ ∂ ∂ ∂ ⋅⋅⋅⋅⋅⋅ ⋅⋅⋅⋅⋅⋅ ∂ ∂⋅⋅⋅∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂⋅⋅⋅∂ ∂ ∂ ∂ ∂ ∂ = n n 3 n 2 n 1 n n 2 3 2 2 2 1 2 n 1 3 1 2 1 1 1 i x f x f x f x f x f x f x f x f x f x f x f x f )X(J và ΔX = Xi+1 - Xi Phương pháp Newton tuyến tính hoá hệ và như vậy với mỗi bước lặp cần giải một hệ phương trình tuyến tính (mà biến là ΔXi) xác định bởi công thức lặp cho tới khi vectơ X(x1, x2, x3,....., xn) gần với nghiệm. Ta xây dựng hàm new4sys() để thực hiện thuật toán này function [P, iter, err] = new4sys(f, jf, P, max1) %vao -F la he pt luu trong M-file f.m % -JF la ma tran jacobi luu trong M-file jf.m 265 % -P vec to nghiem ban dau % -max1 so lan lap cuc dai %ra -P la ve to nghiem % -iter so lan lap thuc te % -err sai so Y = f(P); for k = 1:max1 J = jf(P); Q = P - (J\Y')'; Z = f(Q); err = norm(Q - P); relerr = err/(norm(Q) + eps); P = Q; Y = Z; iter = k; if (err<eps)|(relerr<eps)|(abs(Y)<eps) break end end Để giải hệ phương trình: 2 2 x xy 2 0 2xy y 3 0 ⎧ + − =⎪⎨ + − =⎪⎩ ta dùng chương trình ctnew4sys.m clear all, clc format long p = [.5, .5]; [n ,ll, ss] = new4sys(@f, @jf, p, 50) Nội dung của f.m: function f = f(p) f = [(p(1)^2 + p(1)*p(2) - 2), (2*p(1)*p(2) + p(2)^2 - 3)]; Nội dung của jf.m: function jf = jf(p) jf = [(2*p(1) + p(2)) p(1) (2*p(1)) (2*p(1) + 2*p(2))]; Ta có thể dùng hàm new4sys2() để thực hiện thuật toán: 266 function root = new4sys2(func, x, maxiter) % Phuong phap Newton-Raphson de tim nghiem % cua he pt fi(x1,x2,...,xn) = 0, i = 1, 2,..., n. % Cu phap: root = new4sys2(func, x, tol) % vao: % func = con tro ham, tra ve [f1, f2,..., fn]. % x = vec to ban dau [x1, x2,..., xn]. % tol = sai so mong muon % ra: % root - nghiem cua he if size(x,1) == 1; x = x'; end % x phai la vec to i = 0; while (i < maxiter) [jac, f0] = jacobi(func, x); dx = jac\(-f0); x = x + dx; root = x; if max(abs(dx)) < eps break; else i = i + 1; end end if i == maxiter fprintf('Khong hoi tu sau %d lan lap\n', maxiter); else fprintf('Hoi tu sau %d lan lap\n', i); end Hàm jacobi() gồm các lệnh: function [jac, f0] = jacobi(func, x) % Tinh ma tran Jacobi va ham f(x). h = 1.0e-4; n = length(x); jac = zeros(n); f0 = feval(func, x); for i =1:n temp = x(i); x(i) = temp + h; 267 f1 = feval(func, x); x(i) = temp; jac(:,i) = (f1 - f0)/h; end Hàm t() gồm các lệnh: function x = t(p) x = [(p(1)^2 + p(2)^2 + p(3)^2 - 14) (p(1)^2 + 2*p(2)^2 - p(3) - 6 (p(1) -3*p(2 )^2 + p(3)^2 + 2)]; Để giải hệ phương trình ta dùng chương trình ctnew4sys2.m: clear all, clc format long p = [1 1 1 ]; r = new4sys2(@t, p, 50) §13. PHƯƠNG PHÁP BROYDEN DÙNG CHO HỆ PHI TUYẾN 1. Phương pháp Broyden: Để giải hệ phương trình phi tuyến tính F([X]) = [0] bằng phương pháp lặp Newton ta cho vec tơ nghiệm ban đầu [P0] và tạo ra dãy [Pk] hội tụ về nghiệm [P], nghĩa là F([P]) = [0]. Khi này ta cần tính ma trận Jacobi của hệ. Việc tính ma trận Jacobi đòi hỏi tính n2 đạo hàm riêng. Đạo hàm của hàm f(x) tại pk có thể tính gần đúng bằng công thức: k k 1 k k k 1 f(p ) f(p )f (p ) p p − − −′ = − nên: k k k 1 k k 1f (p )(p p ) f(p ) f(p )− −′ − = − Mở rộng cho hệ phương trình ta có thể viết: J([Pk])([Pk]-[Pk-1]) = F([Pk]) - F([Pk-1]) Phương pháp Broyden bắt đầu bằng việc tính ma trận Jacobi A0 = J([P0]). Sau đó trong quá trình lặp liên tiếp ta dùng ma trận Jacobi được cập nhật Ak: Ak([Pk] - [Pk-1]) = F([Pk]) - F([Pk-1]) Thuật toán Broyden như vậy bao gồm các bước: Bước 0: Tính ma trận Jacobi ban đầu A0 = J([P0]). Sử dụng nó để tính lần lặp đầu tiên theo phương pháp Newton [P1] = [P0] - (A0)-1F([P0]) Với k ≥ 1, giả sử đã biết [Pk] sử dụng các bước sau tính [Pk+1] Bước 1: Tính hàm [ ] 1 k kk 2 k k f (p ,q ) F( P ) f (p ,q ) ⎡ ⎤= ⎢ ⎥⎣ ⎦ Bước 2: Cập nhật ma trận Jacobi bằng cách dùng 268 [S] = [Pk] - [Pk-1] và [Yk] = F([Pk]) - F([Pk-1]) và có: [ ] [ ] [ ] [ ] [ ] [ ][ ]( )[ ]Tk k 1 k k 1T1A A Y A S SS S− −= + − Bước 3: Giải hệ phương trình tuyến tính [Ak]Δ[Pk] = -F([Pk]) để tìm Δ[Pk] Bước 4: Lặp tiếp với [Pk+1] = [Pk] + Δ[Pk] Để thực hiện thuật toán này ta xây dựng hàm broyden() function root = broyden(g, p0, maxiter) if nargin == 2 maxiter = 50; end if size(p0, 1) == 1; p0 = p0'; end [a0, f0] = jacobi(g, p0); p1 = p0 - inv(a0)*f0; i = 0; while i < maxiter s = p1 - p0; d = (s')*s; f1 = feval(g, p1); y = feval(g, p1) - feval(g, p0); m = y - a0*s; a1 = a0 - (1/d)*(y - a0*s)*s'; dp = -inv(a1)*f1; p2 = p1 + dp; err = max(abs(dp)); if err < eps root = p2; break; end p0 = p1; p1 = p2; f0 = f1; a0 = a1; i = i + 1; end Để giải hệ phương trình 269 2 3 4 4 2x 2y 4x 1 0 x 4y 4y 4 0 ⎧ − − + =⎪⎨ + + − =⎪⎩ ta dùng chương trình ctbroyden.m clear all, clc format long p = [.1 .7 ]; root = broyden(@g, p, 50) Ta có thể dùng hàm broyden1( ): function [t, k] = broyden1(fcn1, fcn2, x0, maxits, tol) % Tim nghiem cua he pt phi tuyen % cu phap [t, k] = broyden1(fcn1, fcn2, x0, maxits, tol) % vao % - fcn1: ham thu nhat % - fcn2: ham thu hai % - x0: nghiem ban dau % - maxiter: so lan lap max % - tol sai so mong muon % ra % - x: nghiem % - k: so lan lap tol = 1e-8; maxiter = 50; if size(x0, 1) == 1 x0 = x0'; end syms x y B = [diff(fcn1, x) diff(fcn1, y);diff(fcn2, x) diff(fcn2, y)]; x = x0(1); y = x0(2); h = inline(fcn1); g = inline(fcn2); f(1:2,1) = [h(x, y);g(x, y)]; B = eval(B); t = [x y]'; for k = 1:maxiter s = B\(-f); t = t + s; fnew = [h(t(1), t(2));g(t(1), t(2))]; 270 u = fnew-f; if abs(fnew-f) < tol break end f = fnew; B = B + ((u - B*s)*s')/(s'*s); end và dùng chương trình ctbroyden1.m clc, clear all syms x y f1 = 2*x^2 - 2*y^3 - 4*x + 1; f2 = x^4 + 4*y^4 + 4*y - 4; [n, l] = broyden1(f1, f2, [.1 .7 ], 50) Ngoài ra ta có một phiên bản khác là hàm broyden2(): function [t, k] = broyden2(f, x0, maxiter, tol) % Tim nghiem cua he pt phi tuyen % cu phap [t, k] = broyden2(fcn1, fcn2, x0, maxits, tol) % vao % - fcn1: ham thu nhat % - fcn2: ham thu hai % - x0: nghiem ban dau % - maxiter: so lan lap max % - tol sai so mong muon % ra % - x: nghiem % - k: so lan lap tol = eps; maxiter = 50; if size(x0, 1) == 1 x0 = x0'; end syms x y B = [diff(f(1), x) diff(f(1), y); diff(f(2), x) diff(f(2), y)]; x = x0(1); y = x0(2); h = inline(f(1)); g = inline(f(2)); 271 f(1:2,1) = double( [h(x, y);g(x, y)]); B = double(eval(B)); t = [x y]'; for k = 1:maxiter s = double(B\(-f)); t = t + s; fnew = [h(t(1),t(2)); g(t(1),t(2))]; u = fnew - f; if abs(double(u)) < tol break end f = fnew; B = B + ((u - B*s)*s')/(s'*s); end và dùng với chương trình ctroyden2.m clc, clear all syms x y f = [ 2*x^2 - 2*y^3 - 4*x + 1; x^4 + 4*y^4 + 4*y - 4]; [n, l] = broyden2(f, [.7 .7 ]) 2. Phương pháp Broyden cải tiến: Ma trận nghịch đảo tính toán rất phức tạp. Vì vậy ta dùng công thức Sherman - Morison để giảm bớt khối lượng tính toán khi nghịch đảo ma trận. Nếu [A]-1 không suy biến và [U], [V] là 2 vec tơ sao cho [V]T[A]-1[U] ≠ -1 thì: [ ] [ ][ ]( ) [ ] [ ] [ ][ ] [ ]( )[ ] [ ] [ ]( )[ ] 1 T 1 1T 1 T 1 1 A U V A A U V A 1 V A U A − − − − − −+ = − + hay [ ] [ ][ ]( ) [ ] [ ] [ ]( )[ ][ ] [ ] [ ]( ) [ ] 1 T 1T 1 T 1 A U V A U V E A 1 V A U − − − − ⎧ ⎫⎪ ⎪+ = −⎨ ⎬+⎪ ⎪⎩ ⎭ Để giải hệ phương trình phi tuyến F([X]) = [0] ta cho vec tơ nghiệm ban đầu [P0] và tạo ra dãy [Pk] hội tụ về [P], nghĩa là F([P]) = [0]. Trước hết ta tính ma trận Jacobi A0 = J([P0]) và dùng nó tính lần lặp thứ nhất theo phương pháp Newton. [P1] = [P0] - [A]-1F([P0]) Giả sử đã có [Pk] với k ≥ 1 ta dùng các bước sau để tính [Pk+1] 272 Bước 1: Tính hàm [ ]k kF F( P )= Bước 2: Cập nhật ma trận Jacobi bằng cách dùng [V] = [Pk] - [Pk-1] [ ] [ ] [ ] [ ] [ ] [ ][ ]( )k k 1 k 1T1U F F A VV V − −= − − và có: [ ] [ ] [ ][ ]Tk k 1A A U V−= + Bước 3: Tính [Bk] = [Ak]-1 bằng cách dùng công thức Sherman - Morison [ ] [ ] [ ] [ ] [ ]( )[ ][ ] [ ] [ ]( ) [ ] 1 T k 11 1 k k k 11T k 1 A U V B A E A 1 V A U − −− − −− − ⎧ ⎫⎪ ⎪= = −⎨ ⎬+⎪ ⎪⎩ ⎭ Bước 4: Lặp tiếp với [Pk+1] = [Pk] - [Bk]F([Pk]) Để thực hiện thuật toán trên ta xây dựng hàm improvedbroyden() function [m, ll] = improvedbroyden(g, p0, maxiter); % cu phap [n, ll] = improvedboyden(g, x0, maxiter); % vao: % - g la fi chua ham can tim nghiem % - p0 la vec to nghiem ban dau % - maxiter so lan lap max %ra: % - m la nghiem % - ll so lan lap thuc te if size(p0, 1) == 1 p0 = p0'; end n = length(p0); [a0, f0] = jacobi(g, p0); b0 = inv(a0); p1 = p0 - b0*g(p0); for k = 1: maxiter f1 = g(p1); v = p1 - p0; d = v'*v; u = (1/d)*(f1 - f0 - a0*v); a1 = a0 + u*v'; e = eye(n); ts = (b0*u)*v'; ms = 1 + v'*b0*u; b1 = (e - ts/ms)*b0; p2 = p1 - b1*g(p1); 273 if abs(max(v)) < eps break; end a0 = a1; b0 = b1; p0 = p1; p1 = p2; f0 = f1; end m = p2; ll = k; Ta giải hệ phương trình: 2 3 4 4 2x 2y 4x 1 0 x 4y 4y 4 0 ⎧ − − + =⎪⎨ + + − =⎪⎩ bằng cách dùng chương trình ctimprovedbroyden.m clear all, clc format long p = [.1 .7 ]; [r, s] = improvedbroyden(@g, p, 50) §14. TÍNH TRỊ SỐ CỦA ĐA THỨC 1. Sơ đồ Horner: Giả sử chúng ta cần tìm giá trị của một đa thức tổng quát dạng: P(x) = a0xn + a1xn - 1 + a2xn - 2 +....+ an (1) tại một trị số x nào đó. Trong (1) các hệ số ai là các số thực đã cho. Chúng ta viết lại (1) theo thuật toán Horner dưới dạng: P(xo) = (...((a0x + a1)x+ a2x)+...+ an -1 )x + an (2) Từ (2) ta có : P0 = a0 P1 = P0x + a1 P2 = P1x + a2 P3 = P2x + a3 .................. P(x) = Pn = Pn-1x + an Tổng quát ta có : Pk = Pk-1x + ak với k = 1, 2...n ; P0 = a0 Do chúng ta chỉ quan tâm đến trị số của Pn nên trong các công thức truy hồi về sau chúng ta sẽ bỏ qua chỉ số k của P và viết gọn P := Px + ak với k = 0...n. Khi ta tính tới k = n thì P chính là giá trị cần tìm của đa thức khi đã cho x. Chúng ta thử các bước tính như sau : Ban đầu P = 0 274 Bước 0 k = 0 P = ao Bước 1 k = 1 P = aox + a1 Bước 2 k = 2 P = (aox + a1)x + a2 ................................. Bước n-1 k = n - 1 P = P(xo) = (...((aox + a1)x+a2x)+...+an-1)x Bước n k = n P = P(xo) = (...((aox + a1)x+a2x)+...+an-1)x + an Ta xây dựng hàm horner() để tính trị của đa thức tại x: function p = horner(a, x) % Tinh tri so da thuc % p = a(1)*xˆn + a(2)*xˆ(n-1) + ... + a(n+1) % cu phap: p = horner(a,x) n = length(a) - 1; p = a(1); for i = 1:n p = p*x + a(i+1); end Để tính trị số của đa thức P3(x) = x3 + 3x2 + 2x - 5 tại x = 1 ta dùng chương trình cthorner.m clear all, clc a = [1 3 2 -5]; p = horner(a, 1) 2. Sơ đồ Horner tổng quát: Giả sử chúng ta có đa thức : Pn(x) = a0xn + a1xn - 1 + a2xn - 2 +....+ an (1) Khai triển Taylor của đa thức tại x = xo có dạng : n 0 0 )n( 2 0 0 0 0 0nn )xx(!2 )x(P)xx( !2 )x(P)xx( !1 )x(P)x(P)x(P −+⋅⋅⋅+−′′+−′+= (2) Mặt khác chúng ta có thể biến đổi đa thức về dạng : Pn(x) = (x - xo)Pn-1(x) + Pn(xo) (3) Trong đó Pn-1(x) là đa thức bậc n - 1 và có dạng : Pn-1 (x) = boxn-1 + bo-1xn - 2 + b2xn - 3 +....+ bn-1 (4) Thuật toán để tìm các hệ số nhận được bằng cách so sánh (1) và (3) : bo = ao bi = ai + bi-1xo bn = Pn(xo) So sánh (2) và (3) ta có : 275 20 0 0 n 1 0 n 0 n 0 0 0 (n) n0 0 P (x ) P (x )(x x )P (x ) P (x ) P (x ) (x x ) (x x ) 1! 2! P (x ) (x x ) 2! − ′ ′′− + = + − + − + ⋅ ⋅ ⋅ + − hay : n 0 0 )n( 2 0 0 0 0 1n0 )xx(!2 )x(P)xx( !2 )x(P)xx( !1 )x(P)x(P)xx( −+⋅⋅⋅+−′′+−′=− − và khi chia hai vế cho (x - x0) ta nhận được : 1n00 )n( 0 00 1n )xx(!2 )x(P)xx( !2 )x(P !1 )x(P)x(P −− −+⋅⋅⋅+−′′+′= (5) So sánh (4) và (5) ta nhận được kết quả : !1 )x(P)x(P 001n ′=− Trong đó Pn-1(x) lại có thể phân tích giống như Pn(x) dạng (3) để tìm ra Pn-1(xo). Quá trình này được tiếp tục cho đến khi ta tìm hết các hệ số của chuỗi Taylor của Pn(x). Tổng quát thuật toán thể hiện ở bảng sau: Pn(x) ao a1 a2 a3 ... an-1 an x = xo 0 boxo b1xo b2xo ... bn-2xo bn-1xo Pn-1(x) bo b1 b2 b3 ... bn-1 bn = Pn(xo) Ta xây dựng hàm genhorner() để thực hiện thuật toán trên function b = genhorner(a, x) % tra ve he so cua da thuc khai trien % c(1)(x-x0)^n + c(2)(x-x0)^(n-1) + ...+ c(n+1) m = length(a) x = 2; for k = 1:m b(1) = a(1); for i = 2:m-k+1 b(i) = b(i - 1)*x + a(i); end c(m-k+1) = b(m-k+1) a = b(1:m-k); end Để khai triển đa thức P(x) = x5 - 2x4 + x3 - 5x + 4 tại x0 = 2 ta dùng chương trình ctgenhorner.m clear all, clc 276 a = [1 -2 1 0 -5 4]; c = genhorner(a, 2) §15. PHƯƠNG PHÁP LAGUERRE Ta xét đa thức bậc n: Pn(x) = a1xn + a2xn-1 + ⋅⋅⋅ + an+1 (1) Nếu đa thức có nghiệm là r thì ta có: Pn(x) = (x - r)Pn-1(x) (2) Trong đó: Pn-1(x) = b1xn-1 + b2xn-2 + ⋅⋅⋅ + bn Cân bằng (1) và (2) ta có: b1 = a1 b2 = a2 + rb1 . . . . . . bn = an +rbn-1 Ta xây dựng hàm deflpoly() để tính các hệ số của đa thức Pn-1(x) function b = deflpoly(a, r) % ha bac da thuc n = length(a) - 1; b = zeros(n, 1); b(1) = a(1); for i = 2:n b(i) = a(i) + r*b(i-1); end Bây giờ ta xét đa thức Pn(x) có nghiệm đơn x = r và (n-1) nghiệm trùng nhau x = q. Đa thức như vậy sẽ được viết thành: Pn(x) = (x - r)(x - q)n-1 (3) Bài toán của ta là cho đa thức (3) dưới dạng: Pn(x) = a1xn + a2xn-1 + ⋅⋅⋅ + an+1 và cần tìm r(chú ý là q cũng chưa biết). Đạo hàm (3) theo x ta có: n 1 n 2nP (x) (x q) (n 1)(x q) (x r) − −′ = − + − − − n 1 n 1P (x) x r x q ⎛ ⎞−= −⎜ ⎟− −⎝ ⎠ Như vậy: n n P (x) 1 n 1 P (x) x r x q ′ −= −− − (4) Từ đó ta có: 277 2 n n 2 2 n n P (x) P (x) 1 n 1 P (x) P (x) (x r) (x q) ′′ ⎡ ′ ⎤ −− = − −⎢ ⎥ − −⎣ ⎦ (5) Ta đặt: n n P (x)G(x) P (x) ′= 2 n n P (x)H(x) G (x) P (x) ′′= − (6) Như vậy (4) và (5) trở thành: n n P (x) 1 n 1G(x) P (x) x r x q ′ −= = +− − (7) 2 n 2 2 n P (x) 1 n 1H(x) G (x) P (x) (x r) (x q) ′′ −= − = +− − (8) Nếu ta giải (7) theo (x - q) và thay kết quả vào (8) ta nhận được phương trình bậc 2 đối với (x - r). Nghiệm của phương trình này là công thức Laguerre: 2 nx r G(x) (n 1) nH(x) G (x) − = ⎡ ⎤± − −⎣ ⎦ (9) Như vậy thuật toán để tìm điểm zero của đa thức tổng quát theo Laguerre là: - cho giá trị ban đầu của nghiệm của đa thức là x - tính Pn(x), nP (x)′ và nP (x)′′ - tính G(x) và H(x) theo (6) - tính nghiem theo (9). Chọn dấu sao cho mẫu số lớn nhất. - cho x = r và tiếp tục lặp từ bước 2 cho đến 5 đến khi |Pn(x)| < ε hay |x - r| < ε. Ta xây dựng hàm laguerre() để thực hiện thuật toán trên function x = laguerre(a, tol) % tim nghiem cua da thuc x = randn; % cho x ngau nhien n = length(a) - 1; for i = 1:30 [p, dp, ddp] = evalpoly(a, x); if abs(p) < tol; return; end g = dp/p; h = g*g - ddp/p; f = sqrt((n - 1)*(n*h - g*g)); if abs(g + f) >= abs(g - f) dx = n/(g + f); else dx =