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

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

การplot กราฟ ที่ง่ายสุดในโลกกก 

----- ใน m-file -----

 fx = @(x) (x-1).^2;

 x = -5:0.2:6;

 y = fx(x);

 plot(x,y)  %คำสั่ง plot กราฟ

----- จบ -----

โปรแกรมมจะแสดงรูป

ต่อไปเป็นการ plot กราฟ 3D และการคอนทัวร์นะฮ้าฟฟ

----- ใน m-file -----

clear;

clc;

fx = @(x1,x2) (x1-3).^2 + (x2-4).^2;

x1 = 0:0.2:6;

x2 = 0:0.2:8;

[x1m,x2m] = meshgrid(x1,x2);

zf = fx(x1m,x2m);

mesh(x1,x2,zf)   % สั่ง plot กราฟ 3D

colorbar    %plot กราฟ หลายๆสี

figure

[cf,hf] = contour(x1,x2,zf,[2 3 4])

clabel(cf,hf)

----- จบ -----

โปรแกรมมจะแสดงรูป 2 รูป พร้อมกันน

การ plot หลายกราฟ ในรูปเดียวของ Matlab

------ ใน m-file ------

clear;
clc;

fx = @(x1,x2) (x1-2).^2 + (x2-2).^2;
x1 = 0:0.2:6;
x2 = 0:0.2:8;

[x1m,x2m] = meshgrid(x1,x2);
zf = fx(x1m,x2m);

figure  %ให้แสดงรูป
[cf,hf] = contour(x1,x2,zf,[2 3 4 5 6 7 8 9 10])
clabel(cf,hf)  % ให้ตั้งชื่อแกน x y

g1 = @(x1,x2)5 - x1 -x2;
zg1 = g1(x1m,x2m);
hold on  % ให้แสดงรูปของโค้ดก่อนหน้าค้างไว้
[cg1,hg1] = contour(x1,x2,zg1,[0 0],'--c')

g2 = @(x1,x2)-2.5 + x1 -x2;
zg2 = g2(x1m,x2m);
[cg2,hg2] = contour(x1,x2,zg2,[0 0],'--m')

----- จบ -----

เมื่อ RUN แล้วในหน้า command window  จะแสดง 

 Columns 28 through 36

    5.1000    5.2000    5.3000    5.4000    5.5000    5.6000    5.7000    5.8000    5.9000
    2.6000    2.7000    2.8000    2.9000    3.0000    3.1000    3.2000    3.3000    3.4000

  Column 37

    6.0000
    3.5000


hg2 =

  460.0011

ซึ่งเป็นค่าที่คำนวณจาก column ต่างๆ และจบที่ค่าสุดท้ายคือ ค่า hf

และจะมีกราฟแสดงขึ้นมา

โจทย์อีกข้อ

------ ใน m-file ------

clear;

clc;

fx = @(x1,x2) (x1-3).^2 + (x2-4).^2;

x1 = 0:0.2:6;

x2 = 0:0.2:8;

[x1m,x2m] = meshgrid(x1,x2);

zf = fx(x1m,x2m);

figure

[cf,hf] = contour(x1,x2,zf,[2 3 4 5 6 7 8 9 10])

clabel(cf,hf)

g1 = @(x1,x2)5 - x1 -x2;

zg1 = g1(x1m,x2m);

hold on

[cg1,hg1] = contour(x1,x2,zg1,[0 0],'--c')

g2 = @(x1,x2)-2.5 + x1 -x2;

zg2 = g2(x1m,x2m);

[cg2,hg2] = contour(x1,x2,zg2,[0 0],'--m')

----- จบ -----

เมื่อ RUN แล้วในหน้า command window  จะแสดง 

Columns 28 through 36

    5.1000    5.2000    5.3000    5.4000    5.5000    5.6000    5.7000    5.8000    5.9000

    2.6000    2.7000    2.8000    2.9000    3.0000    3.1000    3.2000    3.3000    3.4000

  Column 37

    6.0000

    3.5000

hg2 =

  243.0027

และแสดงกราฟ 

จบแบ้วววจ้าา าาา  ^^

งานเก่า MATLAB

MATLAB  23/09/55
MATLAB::fNewton
MAIN FUNCTION::

function [xn,fxn] = fnewton (fx,dfx,x0)

k = 1;
xn = x0 - fx(x0)/dfx(x0);
disp([k x0 xn fx(xn)]);

ep = 0.001;

maxiter = 100;
while abs(xn-x0) > ep & k <= maxiter
%while k <= maxiter
x0 = xn;
k = k+1;
xn = x0 - fx(x0)/dfx(x0);
disp([k x0 xn fx(xn)]);
end
fxn = fx(xn);

RUN IN NEXT THIS::
clear
clc
fx = @(x) x^6 - x -1;
dfx = @(x) 6*x^5 - 1;

x0 = 1.5;
[x,y] = fNewton(fx,dfx,x0);

ต่อ

การบ้าน clear
clc
K = 0.016;
Ca0 = 42;
Cb0 = 28;
Cc0 = 4;
fx = @(x) (Cc0+x)-K*((Ca0-2*x)^2)*(Cb0-x);
dfx = @(x) (2*(8*x - 168)*(x - 28))/125 + (2*(2*x - 42)^2)/125 + 1;

x0 = 5;
[x,y] = fnewton(fx,dfx,x0);


ต่อ

**MATLAB**

clear
clc
fx = @(x) x^6 - x -1;
dfx = @(x) 6*x^5 - 1;

x0 = 1.5;
k = 1;
xn = x0 - fx(x0)/dfx(x0);

disp([k x0 xn fx(xn)]);

ep = 0.001;
maxiter = 100;
while abs(xn-x0) > ep & k <= maxiter
%while k <= maxiter
x0 = xn;
k = k+1;
xn = x0 - fx(x0)/dfx(x0);
disp([k x0 xn fx(xn)]);
end

แปะของเก่า ^^

Linear regression
clear, clc
data = [ 10  20  30  40  50  60   70  80
         25  70  380 550 610 1220 830 1450];
x = data(1,:);
y = data(2,:);
plot(x,y,'or')

n = length(x); %number of data point

sum(x.*y) %each of "x" multiply by each of "y"
a1 = ( n*sum(x.*y) - sum(x)*sum(y))/( n*sum(x.*x)-sum(x)^2)
ybar = sum(y)/n;  xbar =sum(x)/n;
a0 = ybar - a1*xbar

yr = a0 +a1*x;
hold on
plot(x,yr,'b-')

st = sum((y-ybar).^2); % square error wrt. ybar
sr = sum((y-yr).^2);   % aquare error wrt. ymodel
rr = (st-sr)/st        % r square

เมื่อ RUN แล้วในหน้า command window  จะแสดง

ans =

      312850


a1 =

   19.4702


a0 =

 -234.2857


rr =

    0.8805

ทำให้เป็น function

function [a0,a1,rr] = fregress_lin(x,y)

n = length(x); %number of data point

a1 = ( n*sum(x.*y) - sum(x)*sum(y))/( n*sum(x.*x)-sum(x)^2);
ybar = sum(y)/n;  xbar =sum(x)/n;
a0 = ybar - a1*xbar;
yr = a0 +a1*x;
st = sum((y-ybar).^2); % square error wrt. ybar
sr = sum((y-yr).^2);   % aquare error wrt. ymodel
rr = (st-sr)/st ;       % r square

ใช้ function ที่ทำไว้มาแก้โจทย์

clear, clc
data = [ 10  20  30  40  50  60   70  80
         25  70  380 550 610 1220 830 1450];
x = data(1,:);
y = data(2,:);
plot(x,y,'or')

[a0,a1,rr] = fregress_lin(x,y);

yr = a0 +a1*x;
hold on
plot(x,yr,'b-')

y10 = a0 + a1*10

เมื่อ RUN แล้วในหน้า command window  จะแสดง

y10 =

  -39.5833

และกราฟ

โจทย์ Enzyme Kinetics
clear,clc

v = [ 1.3    1.8   3     4.5     6      8     9];
s = [ 0.07  0.136  0.22  0.275  0.335  0.35  0.36];

[a0,a1,rr] = fregress_lin(s,v)
vr = a0 + a1*s;
plot (s,v,'or',s,vr,'b-')


เมื่อ RUN แล้วในหน้า command window  จะแสดง
a0 =

   -1.4240


a1 =

   24.9529


rr =

    0.8808


และแสดงกราฟ

แก้โจทย์เดิม โดยใช้ Michaelis Menten model

clear,clc

v = [ 1.3    1.8   3     4.5     6      8     9];
s = [ 0.07  0.136  0.22  0.275  0.335  0.35  0.36];

[a0,a1,rr] = fregress_lin(s,v)
vr = a0 + a1*s;
plot (s,v,'or',s,vr,'b-')

%Michaelis Menten model
y = 1./v;
x = 1./s;
[a0M, a1M, rr] = fregress_lin(x,y)
vm = 1/a0M; ks = vm*a1M;
v_Mic = vm.*s./(ks+s)
plot(s,v,'or',s,vr,'b-',s,v_Mic,'g--')


เมื่อ RUN แล้วในหน้า command window  จะแสดง


a0 =

   -1.4240


a1 =

   24.9529


rr =

    0.8808


a0M =

    0.0133


a1M =

    0.0570


rr =

    0.9224


v_Mic =

    1.2090    2.3138    3.6729    4.5355    5.4532    5.6789    5.8285
และกราฟ

ในโปรแกรม Matlab มี ฟังก์ชั่นการใช้ curve fitting อยู่แล้ว โดยการ พิมพ์ cftool ในหน้า command window จะแสดงหน้าต่าง

Interpolation
Polynomial interpolation
Example : Determining polynomial coefficient


clear,clc

%data
T = [ 300 400 500];
p = [ 0.616  0.525  0.457];

A = [T(1)   1
     T(2)   2];
b = [p(1)
     p(2)];

 a = inv(A)*b;

 TT = 350;  %if T = 350 degree
 pp = a(1)*TT + a(2)


เมื่อ RUN แล้วในหน้า command window  จะแสดง

pp =

    0.5705

เปรียบเทียบ 1st กับ 2nd order


clear,clc

%data
T = [ 300 400 500];
p = [ 0.616  0.525  0.457];

% 1st order
A = [T(1) 1
     T(2) 1];
b = [p(1)
     p(2)];

a = inv(A)*b;

TT = 350;
pp = a(1)*TT + a(2)

%compare to 2nd order
A = [T(1)^2  T(1)  1
     T(2)^2  T(2)  1
     T(3)^2  T(3)  1];
b = [p(1)
     p(2)
     p(3)];
a = inv(A)*b;
TT=350;
pp = a(1)*TT^2 + a(2)*TT + a(3)

เมื่อ RUN แล้วในหน้า command window  จะแสดง


pp =

    0.5705 (1 order)


pp =

    0.5676 (2 order)

นำไปทำเป็น function

function [a1,a0,yi] = finterp_lin(x,y,xi)
A = [x(1)   1
     x(2)   1];
b = [y(1)
      y(2)];
a = inv(A)*b;
a1 = a(1); a0=a(2);

yi = a(1)*xi + a(2);

ใช้ function ที่ทำไว้มาแก้โจทย์


clear,clc

%data
T = [ 300 400 500];
p = [ 0.616  0.525  0.457];
%find P at T=350
[a0,a1,p] = finterp_lin(T,p,350)

เมื่อ RUN แล้วในหน้า command window  จะแสดง


a0 =

 -9.1000e-004


a1 =

    0.8890


p =

    0.5705

การใช้ spline (ดูเอาเอง)

clear,clc

%data
T = [ 300 400];
p = [ 0.616  0.525];
%find P at T=350
[a0,a1,pr] = finterp_lin(T,p,350)
T = [300 400 500];
p = [0.616  0.525  0.457];
pr = interp1(T,p,350,'spline')
   



เมื่อ RUN แล้วในหน้า command window  จะแสดง
a0 =

 -9.1000e-004


a1 =

    0.8890


pr =

    0.5705


pr =

    0.5676

เสร็จจ้าาาาา!!!!

Ordinary differential 
1. Euler's method if  * = f(t,y)
so, y(i+1) = y(i) + f(ti,yi)h ; h is step size

Problem ::   dx/dt = -x +10 = f(ti,xi)
                  x(i+1)  = x(i) + f(ti,xi)*h

พิมสูตรใส่ m-file


clear, clc
%problem
dx = @(x) - x +10;

t(1) = 0;   %t ที่เวลาต่าง
x(1)= 1;   %initial value
h = 1;      %step size

time = 10;
n = time/h;
% make for loop
for i = 1:1:n   %start from 1: increase 1 : to n
    t(i+1) = t(i) + h;
    x(i+1) = x(i) + dx(x(i))*h;
end

plot(t,x)   %วาดกราฟ

เมื่อ RUN แล้ว กราฟจะได้กราฟขึ้นมา
จะได้ คำตอบ x ในหน้า command  window
>> x

x =

     1    10    10    10    10    10    10    10    10    10    10

การพล็อต 2 กราฟ

clear, clc
%problem
dx = @(x) - x +10;

t(1) = 0; %t ที่เวลาต่าง
x(1)= 1;  %initial value
h = 1;    %step size

time = 10;
n = time/h;
% make for loop
for i = 1:1:n %start from 1: increase 1 : to n
    t(i+1) = t(i) + h;
    x(i+1) = x(i) + dx(x(i))*h;
end

plot(t,x, 'b')   %ให้แสดงกราฟเป็นเส้นทึบ


h = 0.01;  
n = time/h;
% make for loop
for i = 1:1:n %start from 1: increase 1 : to n
    t(i+1) = t(i) + h;
    x(i+1) = x(i) + dx(x(i))*h;
end

hold on     %ให้กราฟเดิมยังอยู่
plot(t,x, 'r--')   %r-- ให้แสดงกราฟเป็นเส้นประสีแดง

การทำ ให้เป็น function ไว้เรียกใช้
function Euler

function [t,x] = fEuler(dx,x0,time,h)

t(1) = 0;
x(1)= x0;
n = time/h;
% make for loop
for i = 1:1:n %start from 1: increase 1 : to n
    t(i+1) = t(i) + h;
    x(i+1) = x(i) + dx(x(i))*h;
end

การเรียกใช้ function Euler ที่เราทำไว้

เปิดหน้า M-file ขึ้นมาใหม่ แล้วใส่ปัญหาของเราไป

clear, clc

dx = @(x) - x +10;

x0= 1;  %initial value

h = 1;    %step size

time = 10;

[t,x] = fEuler(dx,x0,time,h);

plot(t,x)

กด RUN แล้วลอง หาX ในหน้า command window

จะได้ คำตอบ x ในหน้า command  window
>> x

x =

     1    10    10    10    10    10    10    10    10    10    10

สร้างสมการจาก system

problem A
สร้างใน M-file

clear,clc

w1 = 500;
w2 = 200;
x1 = 0.4;
x2 =0.75;

model_ss = @(x) w1*(x1-x)+ w2*(x2-x);
[x, fval,exitflag] = fzero(model_ss,1);
%1 equation, 1 variable :: must be use fzero to solve
%check exitflag
***note 1 สมการ 1ตัวแปร ใช้ fzero ในการคำนวณ
แต่ถ้าเป็น หลายสมการ หลายตัวแปร จะใช้ fsolve แทนนะจ๊ะ ***

เชคคำตอบที่หน้าคอมมานวินโดว์

>> x

x =

    0.5000

>> [x, fval,exitflag]

ans =

    0.5000    0.0000    1.0000

Problem B
ใน M-File

%problem b : w1 changes from 500 to 400
w1 = 400;
w2 = 200;
x1 = 0.4;
x2 =0.75;

model_ss = @(x) w1*(x1-x)+ w2*(x2-x);
x0 = 0.5;
[x, fval,exitflag] = fzero(model_ss,x0);

คำตอบใหม่คือ x ลู่เข้าสู่ ...เช็คได้จากหน้า คอมมานวินโดว์

>> x

x =

    0.5167

Pb.A และ Pb.B เป็นการคิดที่ steady state
ถ้าไม่เป็น SS ใช้ ODE คำนวณ

การสร้าง ODE

clear, clc
w1 = 400;
w2 = 200;
x1 = 0.4;
x2 =0.75;

model_ss = @(x) w1*(x1-x)+ w2*(x2-x);
x0 = 0.5;
[x, fval,exitflag] = fzero(model_ss,x0);

v=2;
p=900;
model_dyn = @(x) w1/v/p*(x1-x) + w2/v/p*(x2-x);
time = 30;
h=0.01;
[t,x] = fEuler(model_dyn,0.5,time,h);
plot(t,x)
** คำตอบดูจากกราฟ ว่า มีแนวโน้มถูกต้องหรือไม่ **

Pb. C ใน M-file

%Pb C w2 chenges 200 to 100
clear, clc
w1 = 500;
w2 = 100;
x1 = 0.4;
x2 =0.75;

model_ss = @(x) w1*(x1-x)+ w2*(x2-x);
x0 = 0.5;
[x, fval,exitflag] = fzero(model_ss,x0);

v=2;
p=900;
model_dyn = @(x) w1/v/p*(x1-x) + w2/v/p*(x2-x);
time = 30;
h=0.01;
[t,x] = fEuler(model_dyn,0.5,time,h);
plot(t,x)


จะได้กราฟคำตอบ

Pb.D ใน M-file

%Pb D x1 chenges 0.4 to 0.6
clear, clc
w1 = 500;
w2 = 200;
x1 = 0.6;
x2 =0.75;

model_ss = @(x) w1*(x1-x)+ w2*(x2-x);
x0 = 0.5;
[x, fval,exitflag] = fzero(model_ss,x0);

v=2;
p=900;
model_dyn = @(x) w1/v/p*(x1-x) + w2/v/p*(x2-x);
time = 30;
h=0.01;
[t,x] = fEuler(model_dyn,0.5,time,h);
plot(t,x)

จะได้กราฟคำตอบเป็น

การใช้ model "ode45" หรือ SK model
ใช้โจทย์ Pb.B
สร้าง function ของ dynamic model ขึ้นมาก่อน


function df = fblending_tank(t,x)
%input
w1 = 400;
w2 = 200;
x1 = 0.4;
x2 =0.75;
%parameter
v = 2;
p = 900;
%dynamic model
df = w1/v/p*(x1-x) + w2/v/p*(x2-x);

เปรียบเทียบ ระหว่าง คำตอบที่ได้จาก Euler และ ode45

clear, clc

w1 = 400;
w2 = 200;
x1 = 0.4;
x2 =0.75;

model_ss = @(x) w1*(x1-x)+ w2*(x2-x);
x0 = 0.5;
[xss, fval,exitflag] = fzero(model_ss,x0);

v=2;
p=900;
model_dyn = @(x) w1/v/p*(x1-x) + w2/v/p*(x2-x);
time = 30;
h=0.01;

[t,x] = fEuler(model_dyn,x0,time,h);
plot(t,x,'b')

[t1,x1] = ode45(@fblending_tank,[0 time],0.5);
hold on
plot (t1,x1,'r--')

จะได้กราฟคำตอบ

จะได้ว่ากราฟสองเส้นทับกันพอดี!
check คำตอบที่หน้า คอมมานวินโดว์

>> x

x =.
.
.
.

Column 3001

    0.5167
สุดท้าย x ลู่เข้าสู่ค่า 0.5167

Example System 2

*มีหลายสมการ ใช้ fsolve*
At  Steady state
สร้าง Function ขึ้นมาก่อน ใน M-file

function df = fTank_heater_ss(x)
h= x(1);
T= x(2);

 %inputs
 Fi = 300; %cm3/s
 Ti = 25; %Oc
 Q = 5000; %J
 %parameter
 Cv = 50;
 %A = 100; %cm2
 %p = 1; %g/cm3
 Cp = 4.184; %J/g.C

 df = [Fi - Cv*sqrt(h)
     Fi/h*(Ti-T) + Q/h/Cp]; %steady state Model

แก้สมการ โดยใช้ fsolve ใน M-file

clear, clc
x0 = [1 1];

[x, fval,exitflag] = fsolve(@fTank_heater_ss,x0);
%fval is function value at x

จะได้คำตอบในคอมมาานวินโดว์

>> x

x =

   36.0000   28.9834

Un_SS
สร้างสมการ ที่ Unsteady state in M-file

function df = fTank_heater_dyn(t,x)
h= x(1);
T= x(2);

 %inputs
 Fi = 300; %cm3/s
 Ti = 25; %Oc
 Q = 5000; %J
 %parameter
 Cv = 50;
 A = 100; %cm2
 p = 1; %g/cm3
 Cp = 4.184; %J/g.C

 df = [Fi - Cv*sqrt(h)
     Fi/A/p/h*(Ti-T) + Q/A/p/h/Cp];   %Unsteady state Model

แก้สมการ โดยใช้ ode45

clear, clc
x0 = [1 1];

[xss, fval,exitflag] = fsolve(@fTank_heater_ss,x0);
%fval is function value at x

x0 = xss;
[t,x] = ode45(@fTank_heater_dyn,[0 20],x0);

จะได้คำตอบในคอมมานวินโดว์

>> x

x =

   36.0191   28.9834
   36.0163   28.9834
   36.0140   28.9834
   36.0120   28.9834 ......

พล็อตกราฟเพื่อดูค่า ใน M-File

clear, clc

x0 = [1 1];

[xss, fval,exitflag] = fsolve(@fTank_heater_ss,x0);

%fval is function value at x

x0 = xss;

[t,x] = ode45(@fTank_heater_dyn,[0 20],x0);

plot(t,x(:,1))

figure

plot(t,x(:,2))

จะได้กราฟ

รูปที่ 1 และ 2 แสดงค่าที่ SS ค่าทีไ่ด้จะเป็นเส้นตรง (ที่เหวียงๆข้างหลัง คือ numerical error)

ลองเปลี่ยนค่า 

 Fi = 200 ดู

จะได้กราฟ 

แสดงลักษณะที่ได้ ^^

เสร็จจ้าาาาา!!!!!

check logic ของ Matlab
>> 1>0

ans =

     1

>> 1<0

ans =

     0

>> 1 < 0 & 2==1

ans =

     0

>> 2 ~= 1

ans =

     1

>> 1 <0

ans =

     0

>>  1<0 | 2~=1

ans =

     1

การแก้สมการ 3 ตัวแปร ของ matlab (สมการจาโคบี)

clear, clc

%linear system of equation
%A is 3*3 metrix
%Ax = b
A= [3  -0.1  -0.2
    0.1   7  -0.3
    0.3  -0.2  10];
%A = { 3  -0.1  -0.2; 0.1   7  -0.3; 0.3  -0.2  10];
%b is 3*1 metrix
b = [7.85; -19.3; 71.4];

%initial guess
x10 = 0; x20 = 0;  x30 = 0;
x1 = ( b(1) - A(1,2)*x20 - A(1,3)*x30 ) / A(1,1);
x2 = ( b(2) - A(2,1)*x10 - A(2,3)*x30 ) / A(2,2);
x3 = ( b(3) - A(3,1)*x10 - A(3,2)*x20 ) / A(3,3);

ep1 = 0.01; ep2 = 0.01; ep3 = 0.01;
OK = abs(x1-x10) < ep1 & abs(x2-x20) <ep2 & abs(x3-x30) < ep3; %check eror
while ~OK
    x10 = x1; x20=x2; x30=x3;
    x1 = ( b(1) - A(1,2)*x20 - A(1,3)*x30 ) / A(1,1);
    x2 = ( b(2) - A(2,1)*x10 - A(2,3)*x30 ) / A(2,2);
    x3 = ( b(3) - A(3,1)*x10 - A(3,2)*x20 ) / A(3,3);
    OK1 = abs(x1-x10) < ep1; %check eror of x1
    OK2 = abs(x2-x20) < ep2; %check eror of x2
    OK3 = abs(x3-x30) < ep3; %check eror of x3
    OK = OK1 & OK2 & OK3; %check eror of all2
end

งาน Reactor system

clear, clc
%ใส่ตัวแปร
q01=5;
q03=8;
q12=3;
q15=3;
q23=1;
q24=1;
q25=1;
q34=8;
q31=1;
q44=11;
q54=2;
q55=2;

%สมการ balance ในถัง
%-(q12+q15)*c1  +0                 +q31*c3        +0       +0            =-5*10
%q12*c1        -(q23+q24+q25)*c2  +0             +0       +0            =0
%0             +(q23)*c2          -(q31+q34)*c3  +0       +0            =-160
%0             +q24*c2            +q34*c3        -q44*c4  +q54*c5       =0
%q15*c1        +q25*c2            +0             +0       -(q54+q55)*c5 =0


A = zeros(5,5);
A(1,1) = -(q12+q15);
A(1,3) = q31;
A(2,1) = q12;
A(2,2) = -(q23+q24+q25);
A(3,2) = q23;
A(3,3) = -(q31+q34);
A(4,2) = q24;
A(4,3) = q34;
A(4,4) = -q44;
A(4,5) = q54;
A(5,1) = q15;
A(5,2) = q25;
A(5,5) = -(q54+q55);

b = zeros(5,1);
b(1)=-50;
b(3)=-160;
เมื่อ RUN แล้ว ในหน้าต่าง command window ให้ใส่

>> x = inv(A)*b
เพื่อหาค่า x
คำตอบจะได้
x =

   11.5094
   11.5094
   19.0566
   16.9983
   11.5094

งาน การใช้ fex ,fsolve
ในหน้าต่าง m file
ให้ใส่สมการเข้าไป


function y = fex_nonlinear(x)

y = [ x(1)^2 + x(1)*x(2) - 10
     x(2) + 3*x(1)*x(2)^2 - 57 ];



เมื่อ RUN แล้ววไปที่หน้าต่าง command window ให้ใส่

>> fex_nonlinear([1 1])

ans =

    -8
   -53

>> x = fsolve(@fex_nonlinear,[1 1])
Optimization terminated: first-order optimality is less than options.TolFun.

x =

    2.0000    3.0000

จากนั้นพิม
>>[x,fval,exitflag] = fsolve(@fex_nonlinear,[1 1])
Optimization terminated: first-order optimality is less than options.TolFun.

x =

    2.0000    3.0000


fval =

  1.0e-013 *

         0
   -0.2842


exitflag =

     1
****exitflag = 1 
1,2,3,4 คำตอบเชื่อถือได้
0,-1,-2,-3,-4 ไม่น่าเชื่อถือ

การให้ matlab ช่วยเรา (help)
ที่หน้าต่าง command window ให้ใส่ help fsolve
รอพักนึกจากนัน้กด doc fsolve ที่มันขึ้นมาให้
จากนั้นจะปรากฎหน้าต่าง แสดงว่า   fsolve ใช้อย่างไร มีความหมายอย่างไร

งาน ตัวอย่างการคำนวณการสกัด

(ใน m file )
function y = fsolve_ext(x)
%จากโจทย์
%1st-stage: cR = y2*S -y1*S -x1*R = 0
%2nd-stage: x1R + 0 - y2S - x2R =0
%y1 = mx1
%y2 = mx2
x1 = x(1);
x2 = x(2);
y1 = x(3);
y2 = x(4);

c = 0.1;
R = 1;
S= 12;
m = 0.125;

y = [c*R + y2*S - y1*S - x1*R
    x1*R - y2*S - x2*R
    y1 - m*x1
    y2 - m*x2];

เมื่อ RUN แล้ววไปที่หน้าต่าง command window ให้ใส่


>> x = fsolve(@fsolve_ext,[1 1 1 1])
จะขี้นว่า
Optimization terminated: first-order optimality is less than options.TolFun.

x =

    0.0526    0.0211    0.0066    0.0026
จากนั้น check exitflag เพื่อ check คำตอบ


>> [x,fval,exitflag] = fsolve(@fsolve_ext,[1 1 1 1])
จะขี้นว่า
Optimization terminated: first-order optimality is less than options.TolFun.

x =

    0.0526    0.0211    0.0066    0.0026


fval =

  1.0e-014 *

   -0.1082
    0.0555
   -0.0153
   -0.0167


exitflag =

     1

เสร็จจ้าาาาา!

i lost my assignment of work. i waned someone to talk to, but he was busy about things. i was just a last thing

Dog day of me

MY DOG DAY ^^

you never ever know me at all. when i need you. you don't stand beside me, don't think about me, don't ever give me your hands. i feel like i live my life alone but i'm weak and tried. i want you to know that i'm not too strong too stand in this world alone since you come into my life. so, don't leave me alone. i give you everything that i have. i surrender you at all.
i love you much my dear.


Pls, want to know me more.
Pls, stay with me. 

all of my post i just want someone special for me see this at some day.^^

my first blog

my first blog ^^