Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Untitled

By: a guest on May 2nd, 2012  |  syntax: MatLab  |  size: 2.65 KB  |  hits: 16  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. % Значения коэффициентов
  2. g = 0.2;    % постоянная шага
  3. d = 0.001;   % точность
  4. l = 0.5;    % коэффициент дробления шага
  5.  
  6. [x, y] = meshgrid(-20:0.1:20, -20:0.1:20);
  7. %Z = c1*X.^2 + c2*Y.^2 + c12*X.*Y;
  8. % func = sym('2*x1^2 + y1^2 + x1*y1')
  9. % func =  sym('x1^2 + y1^2')
  10. syms x1 y1
  11.  func = sym('(1-x1)^2 + 100*(y1-x1^2)^2');
  12. %func = sym('-2*x1^2 - y1^2 + x1*y1')
  13. %func = sym('sin(((1/2)*x1^2)-((1/4)*y1^2) + 3)*cos(2*x1 + 1 - exp(y1))');
  14.  
  15. func_hl = matlabFunction(func);
  16.  
  17. derivateX = diff(func, x1);
  18. derivateX = matlabFunction(derivateX);
  19.  
  20. derivateY = diff(func, y1);
  21. derivateY = matlabFunction(derivateY);
  22.  
  23. z = func_hl(x,y);
  24.  
  25.  % Начальная точка
  26. sx1 = 4; sx2 = 3;
  27. k = 1;  % Счетчик шагов
  28. kmax = 90; % Предельное число шагов,  
  29.             % задается для предотвращения зацикливания
  30.  % Массивы для хранения промежуточных координат
  31. x1trace = [sx1]; x2trace = [sx2];
  32. i = 2;
  33.  % Градиенты
  34.  
  35. hold on;
  36. [C, h] = contour(x, y, z);
  37. clabel(C, h);    % Отображение меток РЅР° линиях СѓСЂРѕРІРЅСЏ
  38. X = sx1; Y = sx2;
  39.  
  40. while k < kmax
  41.    
  42.     derX = derivateX(X, Y);
  43.     derY = derivateY(X, Y);
  44.      
  45.     X = X - g*derX;
  46.     Y = Y - g*derY;
  47.    
  48.     tX = X;
  49.     tY = Y;
  50.     z = func_hl(x,y);
  51.     sz = func_hl(tX, tY);
  52.    
  53.     if(sz < z)
  54.     x1trace(i) = X;
  55.     x2trace(i) = Y;
  56.     i = i + 1
  57.     plot(X, Y, '+');
  58.     if sqrt(derX^2 + derY^2) <= d
  59.         break;
  60.          end
  61.         k = k + 1
  62.     else
  63.         X = X + g*derX;
  64.         Y = Y + g*derY;
  65.         g = l*g;
  66.     end  
  67. end
  68.  
  69. plot(x1trace, x2trace, '-+');
  70. % Вывод начальной точки на график
  71. text(x1trace(1) + 0.2, x2trace(1) + 0.5, 'M0');
  72. % Вывод решения на график
  73. text(X + 2.5, Y, ...
  74.     strvcat(['x1 = ' num2str(X)], ...
  75.             ['x2 = ' num2str(Y)], ...
  76.             ['k  = ' num2str(k)]));
create a new version of this paste RAW Paste Data