Bài giảng Các phương trình phi tuyến

241 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) < 0 thì có ít  nhất một nghiệm của phương trình trong đoạn [x1, x2]. Nếu đoạn [x1, x2] đủ  nhỏ  thì  trong đạon đó sẽ có một nghiệm duy nhất. Như vậy  ta có  thể phát  hiện ra nghiệm bằng cách tính trị của hàm trên các đoạn ∆x và xem chúng có  đổi dấu không.   Ta xây dựng hàm rootsearch() để tìm khoảng chứa nghiệm.   function [x1,x2] = rootsearch(func,a,b,dx)  % Tim doan chua nghiem cua ham f(x).  % Cu phap: [x1,x2] = rootsearch(func,a,d,dx)  % func = ham f(x).  % a,b = daon tim.  % dx = khoang tang  % x1,x2 = doan chu nghiem (a,b);  % dat la NaN neu khong thay nghiem  242 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 :  243      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.  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 < 4     maxiter = 100;   end  if nargin < 3      tolx = 1e‐6;   end  xx(1) = x0;  for k = 2:maxiter      xx(k) = feval(g, xx(k ‐ 1));       err = abs(xx(k) ‐ xx(k ‐ 1));       if err < tolx           break;       end  x1 xo x1 xo 244 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) < 0.  Chia  đoạn  [a,  b]    thành  2  phần  bởi  chính  điểm chia (a + b)/2.    1.   Nếu f((a+b)/2) = 0 thì ξ = (a+b)/2    2.  Nếu  f((a  +  b)/2)  ≠  0  thì  chọn   [a,(a+b)/2] hay  [(a + b)/2, b] mà giá  trị hàm  hai đầu trái dấu và kí hiệu là [a1,b1]. Đối với  [a1, b1] ta lại tiến hành như [a, b].  Ta xây dựng hàm bisection() thực hiện thuật toán trên  function [x,err,xx] = bisection(f, a, b, tolx, maxiter)  %bisection.m de giai pt f(x) = 0 bang phuong phap chia doi cung  %vao: f = ham can tim nghiem  % a/b = bien cua doan can tim nghiem  % tolx = sai so mong muon  % maxiter lan lap max  %ra: x = nghiem  % err = sai so  % xx = cac gia tri trung gian  y x a bξ  b1  245 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) < tol | abs(err) < tolx          break;      elseif fx*fa > 0          a = xx(k);           fa = fx;      else b = xx(k);      end  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ìm nghiệm của hàm f(x) = tg(π ‐ ) ‐ x ta dùng chương trình ctbisection.m  clear all, clc  f = inline(ʹtan(pi ‐ x) ‐ xʹ);  [x, ss, xx] = bisection(f, 1.6, 3, 1e‐4, 50)  §4. PHƯƠNG PHÁP DÂY CUNG Giả sử f(x) liên tục trên trên đoạn [a, b] và f(a).f(b) < 0. Cần tìm nghiệm  của  f(x) = 0. Để xác định  ta xem  f(a)  0. Khi đó  thay vì chia đôi  đoạn [a, b] ta chia [a, b] theo tỉ lệ ‐f(a)/f(b). Điều đó cho ta nghiệm gần đúng :          x1 = a + h1  246 Trong đó   1 f(a) (b a)h f(a) f(b)= − −− +   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  dấu ta được nghiệm gần đúng x2 v.v.  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ó           1 af(b) bf(a)x f(b) f(a) −= −                   (1)  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)));  a b x y x1  ξ  247     if abs(fx) < tolfun | err<tolx          break;      elseif fx*fa > 0          a = xx(k);           fa = fx;      else           b = xx(k);           fb = fx;      end  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ìm nghiệm của hàm f(x) = tg(π ‐x) ‐ x ta dùng chương trình ctchord.m  clear all, clc  f = inline(ʹtan(pi ‐ x) ‐ xʹ);  [x, ss, xx] = falsp(f, 1.7, 3, 1e‐4, 50)  §5. PHƯƠNG PHÁP NEWTON ‐ RAPHSON  Phương  pháp  lặp Newton(còn  gọi  là  phương  pháp  tiếp  tuyến)được  dùng nhiều vì nó hội tụ nhanh. Tuy nhiên phương pháp này đòi hỏi tính fʹ(x).  Công thức Newton ‐ Raphson được suy từ khai triển Taylor của f(x) lân cận x:  2 i 1 i i i 1 i i 1 if(x ) f(x ) f (x )(x x ) O(x x )+ + +′= + − + −         (1)  Nếu xi+1 là nghiệm của phương trình f(x) = 0 thì (1) trở thành:    2i i i 1 i i 1 i0 f(x ) f (x )(x x ) O(x x )+ +′= + − + −           (2)  Giả sử rằng xi gần với xi+1, ta có thể bỏ qua số hạng cuối trong (2) và có công  thức Newton ‐ Raphson:    ii 1 i i f(x )x x f (x )+ = − ′                   (3)  Nếu xi+1 là nghiệm đúng của phương trình thì sai số là ei = x ‐ xi. Khi nghiệm  được tính theo (3) thì sai số là:  248   2ii 1 i i f (x )e e 2f (x )+ ′′= − ′   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:    ‐ cho xo    ‐ tính  f(x)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  a b = xo x1  y x 249     end      dx = ‐fx/dfdx;      xx(k+1) = xx(k) + dx; %pt.(3)      fx = feval(f, xx(k + 1));      if abs(fx)<tolf | abs(dx) < tolx,           break;       end  end  x = xx(k + 1);  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  nghiệm  của  hàm  cho  trong  ví  dụ  trên  ta  dùng  chương  trình  ctnewraph.m  clear all, clc  f = inline(ʹx.^3 ‐ 10*x.^2 + 5ʹ);  [x, ss, xx] = newtonraphson(f, 0.7, 1e‐4, 50)  §6. PHƯƠNG PHÁP CÁT TUYẾN    Phương pháp cát tuyến có thể coi là biến thể của phương pháp Newton  ‐ Raphson theo nghĩa đạo hàm được thay bằng xấp xỉ:   k k 1 k k 1 f(x ) f(x )f (x) x x − − −′ ≈ −                 (1)  và tốn ít thời gian tính hơn khi dùng đạo hàm giải tích hay đạo hàm số.    Bằng cách xấp xỉ, biểu thức:    kk 1 k k f(x )x x f (x )+ = − ′   trở thành:    −+ − ⎡ ⎤−= − = −⎢ ⎥−⎣ ⎦ k k 1 k k 1 k k k k k 1 k x x f(x )x x f(x ) x f(x ) f(x ) dfdx          (2)  với  k k 1k k k 1 f(x ) f(x )dfdx x x − − −= −   250 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  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 <= 1          dfdx = (feval(f,xx(k) + h) ‐ feval(f,xx(k) ‐ h))/h2;      else           dfdx = (fx ‐ fx0)/dx;  x0 x1 x2 f(x)  x  y  a  x0 x1  f(x)  x  y  b  251     end      dx = ‐fx/dfdx;      xx(k + 1) = xx(k) + dx; %pt.(2)      fx0 = fx;      fx = feval(f, xx(k+1));      if abs(fx) < tolfun | abs(dx) < tolx,           break;       end  end  x = xx(k + 1);  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  nghiệm  của  hàm  f(x)  =  x3  ‐  10x2  +  5  ta  dùng  chương  trình  ctsecant.m  clear all, clc  f = inline(ʹx.^3 ‐ 10*x.^2 + 5ʹ);  [x, ss, xx] = secant(f, 0.5, 1, 1e‐4, 50)      §7. PHƯƠNG PHÁP BRENT    Phương  pháp  Brent  kêt  hợp  phương  pháp  chia  đôi  cung  và  phương  pháp nội  suy bậc hai  để  tạo  ra một phương pháp  tìm nghiệm  của phương  trình f(x) = 0 rất hiệu quả. Phương pháp này dùng khi đạo hàm của f(x) khó  tính hay không  thể  tính được. Giả sử  ta cần  tìm nghiêm  trong đoạn  [x1, x2].  Quá trình tìm nghiệm bắt đầu bằng việc chia đôi đoạn [x1, x2] bằng điểm x3.  x2  x x3  x x1  x x3  x x2  x1  252 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à:    2 3 1 3 1 21 2 3 1 2 1 3 2 1 2 3 3 1 3 2 (f f )(f f ) (f f )(f f ) (f f )(f f )x(y) x x x (f f )(f f ) (f f )(f f ) (f f )(f f ) − − − − − −= + +− − − − − −   Cho y = 0 ta có:    2 3 1 2 3 1 3 2 3 1 1 2 3 1 2 1 2 2 3 3 1 f f x (f f ) f f x (f f ) f f x (f f )x (f f )(f f )(f f ) − + − + −= − − − −         (1)  Độ thay đổi của nghiệm là:    3 1 2 2 3 1 2 1 2 3 1 2 3 13 3 1 2 2 3 3 1 x (f f )(f f f ) f x (f f ) f x (f f )x x x f (f f )(f f )(f f ) − − + + − + −∆ = − = − − −   (2)  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 < 4;       tol = 1.0e6*eps;   end  % lan chia doi dau tien  x1 = a;   f1 = feval(f,x 1);  if f1 == 0;       root = x1;       return;   end  x2 = b;   f2 = feval(f, x2);  if f2 == 0;       root = x2;       return;   253 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) < tol          root = x3;           return      end  % xac dinh doan chua nghiem.      if f1*f3 < 0.0;           b = x3;      else           a = x3;      end      if (b ‐ a) < tol*max(abs(b),1.0)          root = 0.5*(a + b);        return      end  % noi suy bac 2.      denom = (f2 ‐ f1)*(f3 ‐ f1)*(f2 ‐ f3);      numer = x3*(f1 ‐ f2)*(f2 ‐ f3 + f1)...              + f2*x1*(f2 ‐ f3) + f1*x2*(f3 ‐ f1);      if denom == 0;           dx = b ‐ a;      else          dx = f3*numer/denom;      end      x = x3 + dx;  %neu lap ra ngoai doan (a,b), dung cach chia doi cung      if (b ‐ x)*(x ‐ a) < 0.0          dx = 0.5*(b ‐ a);           x = a + dx;  254     end  % cho x = x3 & chon x1, x2 moi sao cho  x1 < x3 < x2.      if x < x3          x2 = x3;           f2 = f3;      else          x1 = x3;           f1 = f3;      end      x3 = x;  end  root = NaN;  Để  tìm  nghiệm  của  phương  trình    x|cos(x)|  ‐  1  =  0  ta  dùng  chương  trình  ctbrent.m  clear all, clc  f = inline(ʹx.*abs(cos(x)) ‐ 1ʹ);  x = brent(f, 0.0, 4,1e‐4)  §8. PHƯƠNG PHÁP NGOẠI SUY AITKEN  Xét phương pháp lặp:    x = f(x)                    (1)  với f(x) thoả mãn điều kiện hội tụ của phép lặp, nghĩa là với mọi x∈ [a, b] ta  có:  | f’(x) | ≤ q < 1                  (2)  Như vậy :    xn+1 = f(xn)                    (3)    xn = f(xn‐1)                    (4)  Trừ (3) cho (4) và áp dụng định lí Lagrange cho vế phả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 |  255       . . . . . . . . . . . . . . . . . . .        | 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 :                   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;  256   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ó.   257 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:  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  h1 h2 f(x)  x0, f0 x1, f1 x2, f2 258 %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.;