Các lệnh trong Matlab

   KC1 = 112/121;   KC2 = 9/121;  h312 = 3*h*[‐1 2 1];  for k = 4:n      p1 = y(k ‐ 3, :) + h34*(2*(F(1, :) + F(3, :)) ‐ F(2, :)); %Pt.(8a)      m1 = p1 + KC1*(c ‐ p); %Pt.(8b)      if nargin(f) > 1          c1 = (‐y(k ‐ 2, :) + 9*y(k, :) + h312*[F(2:3, :);          feval(f, t(k + 1), m1)ʹ])/8; %Pt.(8c)      else          c1 = (‐y(k ‐ 2, :) + 9*y(k, :) + h312*[F(2:3, :);          feval(f, t(k + 1))ʹ])/8; Pt.(8c)      end      y(k+1,:) = c1 ‐ KC2*(c1 ‐ p1); %Pt.(8d)      p = p1; c = c1; %cap nhat cac gia tri du bao/hieu chinh      if nargin(f) > 1          F = [F(2:3, :); feval(f, t(k + 1),y(k + 1,:))ʹ];      else  383         F = [F(2:3, :); feval(f, t(k + 1))ʹ];      end  end  Để giải phương trình ta dùng chương trình cthamming.m:  clear all, clc  a = 0;  b = 1;  y = @f1;  ya = [0 1 1]ʹ;  n = 10;  tic  [t, y] = hamming(y, a, b, ya, n);  toc  plot(t, y)  §9. PHƯƠNG PHÁP MILNE    Quá trình tính toán nghiệm được thực hiện qua ba bước:    ‐ Tính gần đúng ym+1 theo công thức (dự đoán):    ( )+ − − −′ ′ ′= + − +(1)m 1 m 3 m 2 m 1 m4hy y 2y y 2y3           (1)    ‐ Dùng  + (1) m 1y  để tính:    + + +′ = (1)m 1 m 1 m 1y f(x ,y )                 (2)    ‐ Dùng  +′m 1y  vừa tính được để tính gầm đúng thứ 2 của ym+1(hiệu chỉnh)     ( )+ − − +′ ′ ′= + + +(2)m 1 m 1 m 1 m m 1hy y y 4y y3             (3)  Ta xây dựng hàm milne() để thực hiện thuật toán trên:  function [t, y] = milne(f, to, tf, y0, n)  h = (tf ‐ to)/n;  y(1, :) = y0ʹ;  ts = to + 3*h;  [t, y] = rungekutta(f, to, ts, y0, 3);  t = [t(1:3)ʹ  t(4):h:tf]ʹ;  for i = 2:4  384     if nargin(f) > 1          F(i ‐ 1, :) = feval(f, t(i), y(i, :));       else          F(i ‐ 1, :) = feval(f, t(i));       end  end  for i = 4:n      p = y(i ‐ 3, :) + (4*h/3)*(2*F(1, :) ‐ F(2, :) + 2*F(3, :)); %Pt.(1)      if nargin(f) > 1          F(4, :) = f(t(i+1), p);%Pt.(2)                else          F(4, :) = f(t(i+1));      end      y(i+1, :) = y(i‐1, :) + (h/3)*(F(2, :) + 4*F(3, :) + F(4, :));%Pt.(3)      F(1, :) = F(2, :);      F(2, :) = F(3, :);      if nargin(f) > 1          F(3, :) = f(t(i+1), y(i+1, :));      else          F(3, :) = f(t(i+1));      end  end  Để giải phương trình ta dùng chương trình ctmilne.m:  clear all, clc  a = 0;  b = 1;  y  = @f2;  ya = 1;  n = 10;  [t, y] = milne(y, a, b, ya, n);  plot(t, y)  §10. BÀI TOÁN GIÁ TRỊ BIÊN  1. Khái niệm chung: Ta xét bài toán tìm nghiệm của phương trình:  385   ′′ ′= = α = βy f(x,y,y ) y(a) ,y(b)   Đây  là bài  toán  tìm nghiệm của phương  trình vi  phân  khi  biết  điều  kiện  biên  và  được  gọi  là  bài  toán giá  trị biên hai  điểm. Trong bài  toán giá  trị  đầu  ta có  thể bắt đầu  tìm nghiệm  ở điểm có các  giá  trị  đầu  đã  cho và  tiếp  tục  cho  các  thời  điểm  sau. Kỹ  thuật này  không  áp dụng  được  cho  bài  toán giá  trị biên vì không có đủ điều kiện đầu ở  các biên để có nghiệm duy nhất. Một cách để khác  phục  khó  khăn  này  là  cho  các  giá  trị  còn  thiếu.  Nghiệm  tìm  được  dĩ  nhiên  sẽ  không  thoả mãn  điều kiện  ở các biên. Tuy nhiên  ta có  thể căn cứ  vào  đó  để  thay  đổi  điều kiện  đầu  trước khi  tích  phân  lại. Phương pháp này cũng giống như bắn  bia. Trước hết  ta bắn  rồi xem  có  trúng  đích hay  không,  hiệu  chỉnh  và  bắn  lại.  Do  vậy  phương  pháp này gọi là phương pháp bắn.    Một phương pháp khác để giải bài toán giá trị biên là phương pháp sai  phân hữu hạn  trong các đạo hàm được  thay bằng các xấp xỉ bằng sai phân  hữu hạn tại các nút lưới cách đều. Như vậy ta sẽ nhận được hệ phương trình  đại số đối với các sai phân.       Cả  hai  phương  pháp  này  có một  vấn  đề  chung:  chúng  làm  tăng  số  phương trình phi tuyến nếu phương trình vi phân là phi tuyến. Các phương  trình này được giải bằng phương pháp lặp nên rất tốn thời gian tính toán. Vì  vạy việc giải bài toán biên phi tuyến rất khó. Ngoài ra, đối với phương pháp  lặp,  việc  chọn  giá  trị  đầu  rất  quan  trọng. Nó  quyết  định  tính  hội  tụ  của  phương pháp lặp.    2. Phương pháp shooting: Ta xét bài toán biên là phương trình vi phân cấp 2  với điều kiện đầu tại x = a và x = b. Ta xét phương trình:  ′′ ′= = α = βy f(x,y,y ) y(a) ,y(b)             (1)  Ta tìm cách đưa bài toán về dạng bài toán giá trị đầu:    ′′ ′ ′= = α =y f(x,y,y ) y(a) ,y (a) u             (2)  Chìa khoá  thành công  là  tìm  ra giá  trị đúng u.  ta có  thể  thực hiện việc này  bằng phương pháp “thử và sai”: cho một giá trị u và giải bài toán giá trị đầu  bằng cách đi từ a đến b. Nếu nghiệm giống với điều kiện biên mô tả y(b) = β  386 thì ta đã có nghiệm của bài toán. Nếu không ta phải hiệu chỉnh u và làm lại.  Tuy nhiên  làm như vậy  chưa hay. Do nghiệm  của bài  toán giá  trị  đầu phụ  thuộc u nên giá trị biên tính được y(b) là hàm của u, nghĩa là:    y(b) = θ(u)                      Do đó u là nghiệm của phương trình:    r(u) = θ(u) ‐ β = 0                  (3)  Trong đó θ(u) gọi  là số dự biên(hiệu số giữa giá  trị  tính  được  và  giá  trị  biên  cho  trước). Phương  trình  (3)  có  thể  gải  bằng  các  phương  pháp  tìm  nghiệm  trong  chương  trước.  Tuy  nhiên  phương  pháp  chia  đôi  cung  đòi  hỏi  tính  toán  lâu  còn  phương pháp Newton ‐ Raphson đòi hỏi tìm đạo  hàm dθ/dt. Do vậy cúng ta sẽ dùng phương pháp  Brent. Tóm lại thuật oán giải bài toán giá trị biên  gồm các bước:    ‐ Mô tả giá trị u1 và u2 vây nghiệm u của (3)    ‐ Dùng phương pháp Brent tìm nghiệm u của (3). Chú ý là mỗi bước lặp  đòi hỏi  tính θ(u) bằng cách giải phương  trình vi phân như  là bài  toán điều  kiện đầu.    ‐ Khi đã có u, giải phương trình vi phân lần nữa để tìm nghiệm   Ta xây dựng hàm bvp2shoot() để giải phương trình bậc 2:  function [t, x] = bvp2shoot(f, t0, tf, x0, xf, n, tol, kmax)  %Giai phuong trinh: [x1, x2]ʹ = f(t, x1, x2) voi x1(t0) = x0, x1(tf) = xf  if nargin < 8      kmax = 10;   end  if nargin < 7      tol = 1e‐8;   end  if nargin < 6      n = 100;   end  dx0(1) = (xf ‐ x0)/(tf ‐ t0); % cho gia tri dau cua xʹ(t0)  y0 = [x0 dx0(1)]ʹ;  [t, x] = rungekutta(f, t0, tf , y0, n); % khoi gan bg RK4  387 e(1) = x(end,1) ‐ xf;   dx0(2) = dx0(1) ‐ 0.1*sign(e(1));  for k = 2: kmax‐1      y1 = [x0 dx0(k)]ʹ;      [t, x] = rungekutta(f, t0, tf, y1, n);      %sai so giua gia tri cuoi va dich      e(k) = x(end, 1) ‐ xf; % x(tf)‐ xf      ddx = dx0(k) ‐ dx0(k ‐ 1);       if abs(e(k))< tol | abs(ddx)< tol          break;       end      deddx = (e(k) ‐ e(k ‐ 1))/ddx;       dx0(k + 1) = dx0(k) ‐ e(k)/deddx;   end  Để giải phương trình:    ′′ ′= +2y 2y 4xyy  với điều kiện biên: y(0) = 0.25, y(1) = 1/3   Đặt:  ′ = 1y y ,  ′′ = 2y y  ta đưa phương trình về hệ phương trình vi phân cấp 1:  =⎧⎨ = +⎩ 1 2 2 1 2 1 y y y 2y 4xy y và biểu diễn nó bằng hàm f7():  function dx = f7(t, x) %Eq.(6.6.5)  dx(1) = x(2);   dx(2) = (2*x(1) + 4*t*x(2))*x(1);  Để giải phương trình ta dùng chương trình ctbvp2shoot.m:  clear all, clc  t0 = 0;   tf = 1;   x0 = 1/4;   xf = 1/3; % thoi gian dau/cuoi va cac vi tri  n = 100;   tol = 1e‐8;   388 kmax = 10;  y = @f7;  [t, x] = bvp2shoot(y, t0, tf, x0, xf, n, tol, kmax);  xo = 1./(4 ‐ t.*t);   err = norm(x(:,1) ‐ xo)/(n + 1)  plot(t,x(:, 1))   3. Phương pháp sai phân hữu hạn: Ta xét phương trình:  ′′ ′= = α = βy f(x,y,y ) y(a) ,y(b)   Ý tưởng của phương pháp này là chia đoạn [a, b] thành n đoạn nhỏ có bước h  và xấp xỉ các đạo hàm bằng các sai phân:    + −′ = i 1 iy yy 2h   − +− +′′ = i 1 i i 12y 2y yy h   Như vậy:    − + +− + −⎛ ⎞′′ = = ⎜ ⎟⎝ ⎠ i 1 i i 1 i 1 i i i2 y 2y y y yy f x ,y , h 2h   i = 1, 2, 3,...  Phương  trình bnày sec đưa đến hệ phương  trình đối với xi, yi. Ta xây dựng  hàm bvp2fdf():  function [t, x] = bvp2fdf(a1, a0, u, t0, tf, x0, xf, n)  % Giai pt : xʺ + a1*x’ + a0*x = u with x(t0) = x0, x(tf) = xf  % bang pp sai phan huu han  h = (tf ‐ t0)/n;   h2 = 2*h*h;  t = t0 + [0:n]ʹ*h;  if ~isnumeric(a1)      a1 = a1(t(2: n));   else      length(a1) == 1      a1 = a1*ones(n ‐ 1, 1);  end  if ~isnumeric(a0)      a0 = a0(t(2:n));   else length(a0) == 1  389     a0 = a0*ones(n ‐ 1, 1);  end  if ~isnumeric(u)      u = u(t(2:n));   elseif length(u) == 1      u = u*ones(n‐1,1);  else       u = u(:);  end  A = zeros(n ‐ 1, n ‐ 1);   b = h2*u;  ha = h*a1(1);   A(1, 1:2) = [‐4 + h2*a0(1) 2 + ha];  b(1) = b(1) + (ha ‐ 2)*x0;  for m = 2:n ‐ 2       ha = h*a1(m);       A(m, m ‐ 1:m + 1) = [2‐ha   ‐4+h2*a0(m)   2+ha];  end  ha = h*a1(n ‐ 1);   A(n ‐ 1, n ‐ 2:n ‐ 1) = [2‐ha   ‐4+h2*a0(n ‐ 1)];  b(n ‐ 1) = b(n‐1) ‐ (ha+2)*xf;  x = [x0 trid(A,b)ʹ xf]ʹ;    function x = trid(A, b)  % giai he pt trdiagonal  n = size(A, 2);  for m = 2:n       tmp = A(m, m ‐ 1)/A(m ‐ 1, m ‐ 1);      A(m, m) = A(m, m) ‐A(m ‐ 1, m)*tmp;      A(m, m ‐ 1) = 0;      b(m, :) = b(m, :) ‐b(m ‐ 1, :)*tmp;  end  x(n, :) = b(n, :)/A(n, n);  for m = n ‐ 1: ‐1: 1       x(m, :) = (b(m, :) ‐A(m, m + 1)*x(m + 1))/A(m, m);  end  390 Để giải phương trình:    ′′ ′+ − =22 2y y 0x x   với y(1) = 5 và y(2) = 3 ta dùng chương trình ctbvp2fdf.m:   clear, clc  x0 = 1;%toa do bien dau   y0 = 5; %gia tri bien dau  xf = 2; %toa do bien cuoi  yf = 3; %gia tri bien cuoi  n = 100;%so buoc tinh  a1 = inline(ʹ2./xʹ, ʹxʹ);   a0 = inline(ʹ‐2./x./xʹ,ʹxʹ);   u = 0;%ve phai cua phupng trinh   [x, y] = bvp2fdf(a1, a0, u, x0, xf, y0, yf, n);  plot(x, y)   §11. PHƯƠNG PHÁP LẶP PICARD    Việc giải bài toán Cauchy:    ′ =y f(t,y) ,  =o oy(t ) y                 (1)  hoàn toàn tương đương với việc tìm nghiệm y(t) của phương trình tích phân:    [ ]= + ∫1 o t o t y(t) y f z,y(z) dz               (2)  Nội dung  của phương  trình Picard  là  thay  cho  việc  tìm  nghiệm  đúng  của  phương trình (2) ta tìm nghiệm gần đúng theo công thức:    [ ]+ = + ∫1 o t n 1 o t y y f z,y(z) dz               (3)  Trong đó để đơn giản ta chọn nghiệm gần đúng đầu tiên là:    yo(x) = yo  Ta xây dựng hàm picard() để thực hiện thuật toán trên:  function g = picard(f, y0, maxiter)  syms x y  g = subs(f, y0);  391 g = int(g, x);  for i = 1:maxiter      g = subs(f, g);      g = int(g, x);  end  g = sort(g);   Để  giải  phương  trình  ta  dùng  chương  trình  ctpicard.m:  clear all, clc  syms x y  %f = 1 + y^2;%Phuong trinh yʹ = 1+y^2  %y0 = 0;%dieu kien dau  f = y ‐ x;  y0 = 2;  g = picard(f, y0, 4)  §12. PHƯƠNG PHÁP GILL    Phương pháp Gill cũng tương tự như phương pháp Runge ‐ Kutta, giá  trị nghiệm tại yn+1 được tính bằng:    + ⎡ ⎤= + + − + + +⎢ ⎥⎦⎣n 1 n 1 2 4 4 1y y k (2 2)k (2 2)k k 6 Với:    =1 n nk hf(x ,y )          ⎡ ⎤= + +⎢ ⎥⎣ ⎦2 n n 1 1 1k hf x h,y k 2 2   ( )⎡ ⎤⎛ ⎞= + + − + + −⎢ ⎥⎜ ⎟⎝ ⎠⎣ ⎦3 n n 1 21 1 2k hf x h,y 1 2 k 1 k2 2 2   ⎡ ⎤⎛ ⎞= + − + +⎢ ⎥⎜ ⎟⎝ ⎠⎣ ⎦4 n n 2 3 2 2k hf x h,y k 1 k 2 2 Ta xây dựng hàm gill() để thực hiên thuật toán trên:  function [x, y] = gill(f, a, b, y0, n)  %Phuong phap Gill de giai phuong trinh yʹ(x) = f(x,y(x))  if nargin < 4 | n <= 0  392     n = 100;   end  if nargin < 3      y0 = 0;   end  y(1,:) = y0(:)ʹ;   h = (b ‐ a)/n;   x = a + [0:n]ʹ*h;  if nargin(f) >1      for k = 1:n          k1 = h*feval(f, x(k), y(k, :));           k1 = k1(:)ʹ;           k2 = h*feval(f, x(k)+h/2, y(k, :)+k1/2);           k2 = k2(:)ʹ;                 k3 = h*feval(f, x(k)+h/2, y(k, :)+(sqrt(2)‐1)*k1/2+(1‐sqrt(2)/2)*k2);           k3 = k3(:)ʹ;          k4 = h*feval(f, x(k)+h, y(k, :) ‐sqrt(2)/2*k2+(1+sqrt(2)/2)*k3);           k4 = k4(:)ʹ;           y(k+1, :) = y(k, :) + (k1 + (2 ‐ sqrt(2))*k2 + (2+sqrt(2))*k3 + k4)/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 ‐ sqrt(2))*f2 + (2+sqrt(2))*f3 + f4)/6;       end  end  Để giải phương trình ta dùng chương trình ctgill.m:  393 clear all, clc  a = 0;  b = 1;  y = inline(ʹx + yʹ);  ya = 0.5;  n = 10;%so lan tinh chi n = 10  [t, u] = gill(y, a, b, ya, n);  plot(t, u);  [l, v] = rungekutta(y, a, b, ya, n);  hold on  plot(l, v, ʹ.rʹ)  §13. PHƯƠNG PHÁP RUNGE – KUTTA – FEHLBERG    Một   phương pháp để bảo đảm độ chính xác của nghiệm của phương  trình vi phân là giả bài toán hai lần với bươc tính là h và 0.5h rồi so sánh kết  quả tại các nút  ứng với h. Tuy nhiên điều này đòi hỏi tăng số lần tính và phải  tính  lại khi giá trị tại các nút khác nhau nhiều. Dùng phương pháp Runge –  Kutta – Fehlberg có thể tránh được khó khăn này. Nó có thủ tục để xác định  liệu kích thước bước tính h đã thích hợp chưa. Tại mỗi bươc tính, hai xấp xỉ  khác nhau của nghiệm được  tính và so sánh với nhau. Nếu chúng sai khác  trong phạm vi sai số cho trước thì nghiệm được chấp nhận. Nếu sai khác còn  lơn hơn sai số cho phép thì bước h được giảm và ta lại tìm nghiệm với h mới.  Mỗi bước tính theo phương pháp Runge – Kutta – Fehlberg đòi hỏi 6 giá trị:    1 i ik hf(t ,y )=      2 i i 1 1 1k hf t h,y k 4 4 ⎛ ⎞= + +⎜ ⎟⎝ ⎠     3 i i 1 2 3 3 9k hf t h,y k k 8 32 32 ⎛ ⎞= + + +⎜ ⎟⎝ ⎠     4 i i 1 2 3 12 1932 7200 7296k hf t h,y k k k 13 2197 2197 2197 ⎛ ⎞= + + − +⎜ ⎟⎝ ⎠     5 i i 1 2 3 4 439 3680 845k hf t h,y k 8k k k 216 513 4104 ⎛ ⎞= + + − + −⎜ ⎟⎝ ⎠     5 i i 1 2 3 4 5 1 8 3544 1859 11k hf t h,y k 2k k k k 2 27 2565 4104 40 ⎛ ⎞= + − + − + −⎜ ⎟⎝ ⎠   Xấp xỉ nghiệm theo phương pháp Runge – Kutta bậc 4:  394 i 1 i 1 3 4 5 25 1408 2197 1y y k k k k 216 2565 4104 5+ = + + + −   Và nghiệm tốt hơn dùng phương pháp Runge – Kutta bậc 5:    i 1 i 1 3 4 5 6 16 6656 28561 9 2z y k k k k k 135 12825 56430 50 55+ = + + + − +   Bước tính tối ưu được xác định bằng sh với s là:  1 1 4 4 i 1 i 1 i 1 i 1 h hs 0.840896 2 z y z y+ + + + ⎛ ⎞ ⎛ ⎞ε ε= =⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟− −⎝ ⎠ ⎝ ⎠ Ta xây dựng hàm rkf() để thực hiện thuật toán trên:  function [t, y] = rkf ( f, t0, tf, x0, parms )  % tim nghiem cua phuong trinh  yʹ(t) = f( t, y ), y(t0) = x0  % dung phuong phap Runge‐Kutta‐Fehlberg bac 4, bac 5  neqn = length(x0);  hmin = parms(1);  hmax = parms(2);  tol  = parms(3);  t(1) = t0;  y(1:neqn, 1) = x0ʹ;  count = 0;  h = hmax;  i = 2;  while( t0 < tf )      if nargin(f) > 1          k1 = h*feval(f, t0, x0);          k2 = h*feval(f, t0 + h/4, x0 + k1/4);          k3 = h*feval(f, t0 + 3*h/8, x0 + 3*k1/32 + 9*k2/32);          k4 = h*feval(f, t0 + 12*h/13, x0 + 1932*k1/2197 ‐...                 7200*k2/2197 + 7296*k3/2197);          k5 = h*feval(f, t0 + h, x0 + 439*k1/216 ‐ 8*k2 +...                 3680*k3/513 ‐ 845*k4/4104);          k6 = h*feval(f, t0 + h/2, x0 ‐ 8*k1/27 + 2*k2 ‐...                 3544*k3/2565 + 1859*k4/4104 ‐ 11*k5/40);      else          k1 = h * feval(f, t0);          k2 = h * feval(f, t0 + h/4);  395         k3 = h * feval(f, t0 + 3*h/8);          k4 = h * feval(f, t0 + 12*h/13);          k5 = h * feval(f, t0 + h);          k6 = h * feval(f, t0 + h/2);        end        r = max(abs(k1/360 ‐ 128*k3/4275 ‐ 2197*k4/75240 +...             k5/50 + 2*k6/55)/h);        q = 0.84*(tol/r)^(1/4);        if ( r < tol )           x0 = x0 + 16*k1/135 + 6656*k3/12825 + 28561*k4/56430 ‐...                   9*k5/50 + 2*k6/55;            t0 = t0 + h;            t(i) = t0;            y(1:neqn, i) = x0ʹ;            i = i + 1;         end;        h = min(max(q, 0.1), 4.0)*h;       if (h > hmax)           h = hmax;         end;        if (t0 + h > tf)                      h = tf ‐ t0;       elseif (h < hmin)                     disp(ʹCan giam kich thuoc buoc tinhʹ);           return;       end;  end;  Để giải phương trình ta dùng chương trình ctrkf.m:  clear all, clc  a = 0;  b = 1;  y  = @f2;  ya = 0.5;  p = [1e‐5 1e‐3 1e‐8];%[hmin  hmax tol]  396 [t, y] = rkf(y, a, b, ya, p)  plot(t, y);  370 CHƯƠNG 8: TỐI ƯU HOÁ §1. KHÁI NIỆM CHUNG VỀ TỐI ƯU HOÁ    Tối ưu hoá là thuật ngữ thường được dùng để cực tiểu hoá hay cực đại  hoá một hàm. Thông thường ta chỉ cần tìm cực tiểu một hàm là đủ. Việc tìm  cực đại của f(x) thực hiện một cách đơn giản bằng cách tìm cực tiểu của hàm  −f(x) . Hàm f là hàm giá trị hay hàm đối tượng, cần được giữ cực tiểu. Biến x  là biến có thể hiệu chỉnh tự do.    Các thuật toán cực tiểu hoá là các thủ thuật lặp đòi hỏi một giá trị ban  đầu của biến x. Nếu f(x) có nhiều cực tiểu địa phương, việc chọn giá trị đầu sẽ  xác định cực tiểu nào được tính. Ta không có cách nào bảo đảm là tìm được  cực tiểu toàn cục.    Các biến  có  thể bị  ràng buộc bằng  các  đẳng  thức hay bất  đẳng  thức.  Phần lớn các phương pháp là tìm cực tiểu không ràng buộc, nghĩa là không có  hạn chế nào đối với biến x. Các bài toán này bao gồm tìm cực tiểu của hàm,  tìm điểm tĩnh ‐ điểm có gradient triệt tiêu. Các bài toán tìm cực tiểu có ràng  buộc khó hơn và thuật toán khá phức tạp.     Trong chương này chúng ta sẽ  lần  lượt xét các thuật toán tìm cực tiểu  không ràng buộc và có ràng buộc.  §2. PHƯƠNG PHÁP TIẾT DIỆN VÀNG   Ta xét bài toán tìm cực tiểu của hàm một biến f(x). Điểm cực tiểu được  xác định  theo điều kiện df/dx = 0. Do có  thể có nhiều điểm cực  tiểu nên  ta  phải vây điểm cực tiểu(xác định lân cận chứa điểm cực tiểu) trước. Thủ thuật  vây điểm cực tiểu khá đơn giản: cho điểm đầu x0 và tính giá trị của hàm đang  đi  xuống  tại  các  điểm  tiếp  theo  x1,  x2,  x3,...  cho  đến  tại  xn hàm  tăng  lại  thì  dừng. Điểm cực tiểu bị vây trong khoảng (xn‐2, xn). Khoảng (xi+1, xi) không nên  chọn là hằng số vì như vậy cần nhiều bước tính. Hợp lí nhất là nên tăng kích  thước bước tính để đạt được cực tiểu nhanh hơn, ngay cả khi cực tiểu bị vây  trong  một  đoạn  khá  rộng.  Ta  chọn  kích  thước  tăng  theo  dạng  hằng  số:  + =i 1 ih ch với  >c 1.  ta xây dựng hàm goldbracket() để vây điểm cực  tiểu của  hàm:  function [a, b] = goldbracket(func, x1, h)  % vay diem cuc tieu cua f(x).  c = 1.618033989;  371 f1 = feval(func, x1);  x2 = x1 + h;   f2 = feval(func,x2);  if f2 > f1      h = ‐h;      x2 = x1 + h;       f2 = feval(func, x2);      if f2 > f1          a = x2;           b = x1 ‐ h;           return      end  end  for i = 1:100      h = c*h;      x3 = x2 + h;       f3 = feval(func, x3);      if f3 > f2          a = x1;           b = x3;           return      end      x1 = x2;       f1 = f2;       x2 = x3;      f2 = f3;  end  error(ʹGoldbracket khong tim thay diem cuc tieuʹ)    Tiết diện vàng là một biến thể của phương pháp chia đôi dùng khi tìm  nghiệm của phương trình f(x) = 0. Giả sử điểm cực tiểu bị vây trong khoảng  (a,  b)  có  độ  dài  h.  Để  thu  nhỏ  khoảng  (a,  b)  ta  tính  giá  trị  của  hàm  tại  = −1x b rh  và  = +2x a rh  như hình vẽ. Nếu f1 = f(x1) lớn hơn f2 = f(x2) như hình  a  thì  cực  tiểu nằm  trong  khoảng  (x1,  b) nếu ngược  lại  cực  tiểu nằm  trong  khoảng (a, x2).  372   Giả sử f1 > f2 , ta đặt a = x1 và và x1 =  x2 và có khoảng (a, b) mới như hình b. Để  thực  hiện  bước  thu  gọn  tiếp  theo  ta  lại  tính giá trị của hàm tại x2 = a + rh’ và lặp  lại quá  trình. Quá  trình  làm việc chỉ nếu  hình a và hình b tương tự, nghĩa  là hằng  số r không đổi khi xác định x1 và x2 ở cả  hai hình. Từ hình a ta thấy:    − = −2 1x x 2rh h   Cùng một khoảng cách đó từ hình b ta có:    x1 ‐ a = h’ ‐ rh’  Cân bằng các khoảng này ta được:    2rh ‐ h = h’ ‐ rh’  Thay h’ = rh và khử h:    2r ‐ 1 = r(1 ‐ r)  Giải phương trình này ta nhận được tỉ lệ vàng:    −= =5 1r 0.618033989... 2 Chú ý là mỗi lần thu gọn khoảng chứa điểm cực tiểu thì khoảng (a, b) giảm tỉ  lệ với r. Điều này làm số lần tính lớn hơn phương pháp chia đôi. Tuy nhiên  phương pháp tỉ lệ vàng chỉ đòi hỏi tính giá trị hàm một lần trong khi phương  pháp chia đôi cần tính giá trị hàm 2 lần. Số lần tính xác định bằng:    − = εnb a r   hay:  ε − ε= = − − ln b an 2.078087n ln b a     h = b ‐ a = 0.382  Ta xây dựng hàm golden() để thực hiện thuật toán này:  function [xmin, ymin] = golden(f, a, b, delta, epsilon)  % a va b la doan tim cuc tieu  % delta sai so cua x  % epsilon sai so cua y  r1 = (sqrt(5) ‐ 1)/2;  r2 = r1^2;  h = b ‐ a;  a bx1  x2  h  rh  2rh‐h  rh  a  a bx1  x2  h’  rh’  rh’  b  373 fa = f(a);  fb = f(b);  x1 = a + r2*h;  x2 = a + r1*h;  f1 = f(x1);  f2 = f(x2);  k = 1;  while (abs(fb‐fa) > epsilon)|(h > delta)     k = k + 1;     if (f1 < f2)        b = x2;        fb = f2;        x2 = x1;        f2 = f1;        h = b ‐ a;        x1 = a + r2*h;        f1 = f(x1);     else        a = x1;        fa = f1;        x1 = x2;        f1 = f2;        h = b ‐ a;        x2 = a + r1*h;        f2 = f(x2);     end  end  dp = abs(b ‐ a);  dy = abs(fb ‐ fa);  p = a;  yp = fa;  if (fb < fa)     p = b;     yp = fb;  end  xmin = p;  374 ymin = yp;  Để tìm cực tiểu của hàm ta dùng chương trình ctgolden.m:  clear all, clc  f = inline(ʹ1.6*x^3 + 3*x^2 ‐2*xʹ);  x = 0;   delta = 1e‐8;  epsilon = 1e‐10;  [a, b] = goldbracket(f, x, 0.2);  [xmin, ymin] = golden(f, a, b, delta, epsilon)  §3. PHƯƠNG PHÁP XẤP XỈ BẬC HAI    Ý tưởng của phương pháp này là:    ‐ xấp xỉ hàm đối  tượng  f(x) bằng một hàm bậc 2 p2(x) qua 3 điểm cho  trước    ‐ cập nhật 3 điểm này bằng cách  thay một  trong 3 điểm bằng cực  tiểu  của hàm p2(x)    Qua 3 điểm:   [ ] [ ] [ ]{ }0 0 1 1 2 2(x ,f(x ) , (x ,f(x ) , (x ,f(x )   x0 < x1 < x2  ta tìm được đa thức nội suy p2(x) và điểm có đạo hàm bằng zero:    − + − + −= = ⎡ ⎤− + − + −⎣ ⎦ 2 2 2 2 2 2 0 1 2 1 2 0 2 0 1 3 0 1 2 1 2 0 2 0 1 f (x x ) f (x x ) f (x x )x x 2 f (x x ) f (x x ) f (x x )         (1)  Đặc biệt, nếu các  điểm tìm được trước đây phân bố đều với khoảng cách h(  nghĩa là x2 ‐ x1 = x1 ‐ x0 = h) thì (1) trở thành:    − + − + − − += = + − + −⎡ ⎤− + − + −⎣ ⎦ 2 2 2 2 2 2 0 1 2 1 2 0 2 0 1 0 1 2 3 0 0 1 20 1 2 1 2 0 2 0 1 f (x x ) f (x x ) f (x x ) 3f 4f fx x h 2( f 2f f )2 f (x x ) f (x x ) f (x x )   (2)  Ta  cập  nhật  3  điểm  theo  cách  này  cho  đến  khi  − ≈2 0x x 0   hay  − ≈2 0f(x ) f(x ) 0  và cực tiểu là x3. Quy tắc cập nhật 3 điểm là:    ‐ Trong trường hợp  < <0 3 1x x x  ta dùng  0 3 1(x , x , x )  hay  3 1 2(x , x , x )  làm 3  điểm mới tuỳ theo f(x3) < f(x1) hay không    ‐ Trong trường hợp  < <1 3 2x x x  ta dùng  1 3 2(x , x , x )  hay  0 1 3(x , x , x )  làm 3  điểm mới tuỳ theo f(x3) ≤ f(x1) hay không.  Quá trình tìm cực tiểu được mô tả trên hình sau:  375 Ta xây dựng hàm optquad() để thực hiện thuật toán này.  function [xo, fo] = optquad(f, x0, tolx, tolfun, maxiter)  %Tim cuc tieu cua f(x) bang phuong phap xap xi bac 2  if length(x0) > 2      x012 = x0(1:3);  else      if length(x0) == 2          a = x0(1);           b = x0(2);      else           a = x0 ‐ 10;           b = x0 + 10;      end      x012 = [a (a + b)/2  b];  end  f012 = f(x012);  [xo, fo] = optquad0(f, x012, f012, tolx, tolfun, maxiter);  function [xo, fo] = optquad0(f, x012, f012, tolx, tolfun, k)  x0 = x012(1);   x1 = x012(2);   x2 = x012(3);  376 f0 = f012(1);   f1 = f012(2);   f2 = f012(3);  nd = [f0 ‐ f2 f1 ‐ f0 f2 ‐ f1]*[x1*x1 x2*x2 x0*x0; x1 x2 x0]ʹ;  x3 = nd(1)/2/nd(2);   f3 = feval(f, x3); %Pt.(1)  if k <= 0 | abs(x3 ‐ x1) < tolx | abs(f3 ‐ f1) < tolfun      xo = x3;       fo = f3;      if k == 0          fprintf(ʹDay co the xem la diem cuc tieuʹ)      end  else      if x3 < x1          if f3 < f1               x012 = [x0 x3 x1];               f012 = [f0 f3 f1];          else               x012 = [x3 x1 x2];               f012 = [f3 f1 f2];          end      else          if f3 <= f1              x012 = [x1 x3 x2];               f012 = [f1 f3 f2];          else               x012 = [x0 x1 x3];               f012 = [f0 f1 f3];          end      end      [xo, fo] = optquad0(f, x012, f012, tolx, tolfun, k ‐ 1);  end  Để tìm điểm cực tiểu ta dùng chương trình ctoptquad.m:  clear all, clc  377 %f = inline(ʹ(x.^2 ‐ 4).^2/8 ‐ 1ʹ);  a = 0;   b = 3;  delta = 1e‐8;  epsilon = 1e‐10;  maxiter = 100;  [xmin, ymin] = optquad(f, [a  b], delta, epsilon, maxiter)  §4. PHƯƠNG PHÁP NELDER ‐ MEAD    Phương pháp Nelder ‐ Mead có thể dùng để tìm cực tiểu của hàm nhiều  biến mà phương pháp  tiết diện vàng hay phương pháp xấp xỉ bậc 2 không  dùng được. Thuật toán Nelder ‐ Mead gồm các bước như sau:  Bước 1: Cho 3 điểm đầu tiên a, b, c với f(a) < f(b) < f(c)  Bước 2: Nếu 3 điểm và các giá trị tương ứng của hàm đủ gần nhau thì ta  coi a là điểm cực tiểu và kết thúc quá trình tính  Bước  3:  Nếu  không  ta  coi  điểm  cực  tiểu  nằm  đối  diện  với  điểm  c  trên  đường  ab(xem  hình  vẽ) và lấy:    e = m + 2(m ‐ c)   với   m = (a + b)/2  và nếu f(e) < f(b) thì lấy:    r = (m + e)/2 = 2m ‐ c  và nếu f(r) <