วันเสาร์ที่ 15 ธันวาคม พ.ศ. 2555

Polynomial Approximation

---m file---


%check equation
clear;
clc;

f = @(x) (x-1).^3;
%initial guess
x1 = -5;
x2 = -2;
x3 = 5;

ep = 0.001;
iter = 1;
while abs(x1-x3) >ep
y = [f(x1) f(x2) f(x3)]';
x = [1 x1 x1^2
     1 x2 x2^2
     1 x3 x3^2];
cc = inv(x)*y;
b = cc(2);
c = cc(3);
xr = -b/2/c;f
iter = iter + 1;
disp([iter x1 x2 x3 xr])
if (xr > x1) & (xr < x2)
    x3 = x2;
    x1 = x1;
    x2 = xr;
else
    x3 = x3;
    x1 = x2;
    x2 = xr;
end
end
%condition to make stop is not good so it take much calculate loops to get
%the result

--- จบ ---
หน้า CM จะแสดงคำตอบ


f =

    @(x)(x-1).^3

  1.0e+004 *

    1.5905    0.0001    0.0001    0.0004    0.0001

ใช้ ฟังก์ชั่น Newton ในการหา optimize solution

--- m file ---


clear;
clc;

fx = @(x) (x-1).^3;
df = @(x) 3*(x-1).^2;
ddf = @(x) 6*(x-1);

ep = 0.001;
x0 = 5;
xn = x0 - df(x0)/ddf(x0);
while abs(xn-x0) >ep
    x0 = xn;
    xn = x0 - df(x0)/ddf(x0);
    disp([x0 xn]);
end
%check min or max
if ddf(xn) > 0
    disp('min point');
else
    disp('max point');
end
%optimize solution
[xn fx(xn)]

--- end ---
ในหน้า CM จะแสดงคำตอบ


 3     2

    2.0000    1.5000

    1.5000    1.2500

    1.2500    1.1250

    1.1250    1.0625

    1.0625    1.0313

    1.0313    1.0156

    1.0156    1.0078

    1.0078    1.0039

    1.0039    1.0020

    1.0020    1.0010

min point

ans =

    1.0010    0.0000

การใช้ Grid search เพื่อหาคำตอบ

--- m file ---


clear;
clc;

f = @(x) (x-3).^2;
al = 0.01;       %result is up to alpha change all result will be change too
%alpha is range if it's lare it's not good, if it's small it's good
xcold = -100;     %initial guess
xc = 1;
iter = 1;
while al > 10^-4                 %cal until alpha > 0.0001 i will OK!
    xcold = xc;
    xl = xc + al*(-1);
    xr = xc + al*(1);

    if f(xl) < f(xc) & f(xl) < f(xr)
        xc = xl;
    end
    if f(xr) < f(xc) & f(xr) < f(xl)
        xc = xr;
    end
    if abs(xc-xcold) > 0           %condition for change alpha to prove result
        al = al*2;
    else
        al = al/2;
    end
    iter = iter +1;
    disp([iter xc])
end

--- end ---
ในหน้า CM จะแสดง

...    32     3

      33     3

      34     3



โจทย์อีกข้อนึง (Grid search)

--- m file ---

clear;clc ;

f = @(x)(x-3).^2;
 al = 4;
 xc_old = -100;
 xc = 1;
 iter = 1;
 while   al > 10^-4        %  abs(xc - xc_old) > 0
    xc_old = xc;
    xl = xc + al*(-1);
    xr = xc + al*(1);

    if f(xl) < f(xc) & f(xl) < f(xr)
        xc = xl;
    end
    if f(xr) < f(xc) & f(xr) < f(xl)
        xc = xr;
    end
    if abs(xc-xc_old) > 0
        al = al*2;
    else
        al = al/2;
    end
       
disp([iter xc al])
iter = iter +1;
 end
 xc

--- end ---
หน้า CM จะแสดง

. . .
   17.0000    3.0000    0.0001

   18.0000    3.0000    0.0001


xc =

     3

การใช้ วิธี??? ในการหาค่า

clear; clc; format compact;format short g ;


f = @(x)(x(1)-3).^2+(x(2)-4).^2;

al = @(x1,x2)(x1.^2+x2.^2-6*x1-8*x2+25)/(2*(x1-3).^2+2*(x2-4).^2);
df1 = @(x1,x2)2*(x1-3);
df2 = @(x1,x2)2*(x2-4);

x0 = [-10 -10]';
xn = [1 1]';
ep = 10^-2; 
iter = 1;
while min(abs(xn-x0)) > ep
    x0 = xn;
    s = -[df1(x0(1),x0(2)) 
          df2(x0(1),x0(2))];

%s = -[df1; df2];
iter = iter +1;
xn = x0 + (al(x0(1),x0(2)))*s;

disp(x0)
end

--- จบ ---
ในหน้า CM จะแสดง

     1
     1
     3
     4

การใช้ ฟังก์ชั่น Newton ในการแก้ปัญหา multi condition

--- m file ---
clear;
clc;

x0 = [10 10]';
xn = [5 5]'; % ' is transpost

iter = 1; maxiter = 50;  
ep = 0.01;

while abs(min(xn-x0)) > ep
    x0 = xn;
    x1 = x0(1);
    x2 = x0(2);
    
    df = [8*x1 - 2*x2
          2*x2 - 2*x1];
    H = [ 8  -2
         -2   2 ];
    
     xn = x0 - inv(H)*df;
    
    iter = iter + 1;
    disp([x0 xn]);
end

--- end ---
ในหน้า CM จะแสดง

     5     0
     5     0
     0     0
     0     0

หรือสามารถเขียน code ได้อีกแบบนึง ซึ่งก็คือ

--- m file---
clear;
clc;

f = @(x1,x2) 4*x1^2 + x2^2 - 2*x1*x2;
x0 = [10 10]';
xn = [5 5]'; % ' is transpost

iter = 1; maxiter = 50;  
epx = 0.01;
epf = 0.01;
cond1 = 1;
cond2 = 1;
while cond1 | cond2
    x0 = xn;
    x1 = x0(1);
    x2 = x0(2);
    
    df = [8*x1 - 2*x2
          2*x2 - 2*x1];
    H = [ 8  -2
         -2   2 ];
    
     xn = x0 - inv(H)*df;
    
    iter = iter + 1;
    disp([iter x0' xn']);
    
    cond1 = abs(min(xn-x0)) > epx;
    cond2 = abs(min(f(x0(1),x0(2))...
        - f(xn(1),xn(2)))) > epf
    disp([cond1 cond2]);
end

--- end ---
ในหน้า CM จะแสดง

cond2 =
     1
     1     1
     3     0     0     0     0
cond2 =
     0
     0     0

 การใช้ function  buildin ของ Matlab

--- m file ---
function fx = fobj(x)
x1 = x(1);
x2 = x(2);
fx = 4*x1^2 + x2^2 - 2*x1*x2;

--- end ---
ในหน้า CM พิม...

>> [x,fval,exitflag]=fminsearch(@fobj,[1 1])
x =
 -2.7148e-006 -3.6605e-005
fval =
  1.1707e-009
exitflag =
     1

หรือ

>> [x,fval,exitflag]=fminunc(@fobj,[1 1])
Warning: Gradient must be provided for trust-region method;
  using line-search method instead. 
> In fminunc at 354
Optimization terminated: relative infinity-norm of gradient less than options.TolFun.
x =
  -8.056e-008 -1.6227e-009
fval =
  2.5701e-014
exitflag =
     1

การนำไปใช้ในการทำ Regression

--- m file ---
function sr = fobj1(x)
a0 = x(1);
a1 = x(2);
x = [ 10 20 30 40 50 60 70 80];
y = [ 25 70 380 550 610 1220 830 1450];

sr = 0;
n = length(y); 
for i = 1:n
sr = sr + (y(i)-a0-a1*x(i))^2;
end


%sr = sum((y-a0-a1*x).^2  can take place in line 7 to line 11

--- end ---
ในหน้า CM พิม...

>> [x,fv,ex] = fminsearch(@fobj1,[1 1])
x =
      -234.29        19.47
fv =
  2.1612e+005
ex =
     1

ถ้าโจทย์มี 3 ตัวแปร

--- m file ---
function sr = fobj2(x)
a0 = x(1);
a1 = x(2);
a2 = x(3);
x = [ 10 20 30 40 50 60 70 80];
y = [ 25 70 380 550 610 1220 830 1450];

sr = 0;
n = length(y); 
for i = 1:n
sr = sr + (y(i)-a0-a1*x(i)-a2*x(i)^2)^2;
end

--- end ---
ในหน้า CM พิม...

>> [x,fv,ex] = fminsearch(@fobj2,[1 1 1])
x =
      -178.48       16.122     0.037202
fv =
  2.1379e+005
ex =
     1

จบ บบ บบบ เย่!!!