PitchUp Labs

A good Symbolic toolbox example

2005-09-28T18:37:00.000-07:00

Problem No 3.23, Engineering Optimization Reklaitis

clc;
close all;
syms x1 x2 
a='x1*x2+(2000/(x1*x2))*(x2+x1)';
m=diff(a,x1);
n=diff(a,x2);
pretty(m)
pretty(n)
l=solve(m,x1)
l1=simplify(subs(n,x1,l(1)));
l2=simplify(subs(n,x1,l(2)));
pretty(l1)
pretty(l2)
simplify(solve(l1,x2))
simplify(solve(l2,x2))
simplify(subs(l(1),x2,solve(l1,x2)))
h=1000/(simplify(subs(l(1),x2,solve(l1,x2)))*simplify(subs(l(1),x2,solve(l1,x2))))

How to use subs in Matlab

2005-09-28T18:29:00.000-07:00

clc;
close all;

syms x1 x2 x3 
a='(x1^2+(x2+1)^2)*(x1^2+(x2-1)^2)';

l1=diff(a,x1);
l2=diff(a,x2);
l1=subs(l1,x1,1);
l1=subs(l1,x2,1)
l2=subs(l2,x1,1);
l2=subs(l2,x2,1)

Symbolic Toolbox

2005-09-28T18:28:00.000-07:00

clc;
close all;

syms x1 x2 x3
a='x1^2+2*x2^2-3*x3^2-6*x1*x2+8*x1*x3-4*x2*x3'

A=[diff(diff(a,x1),x1),diff(diff(a,x1),x2),diff(diff(a,x1),x3);...
    diff(diff(a,x2),x1),diff(diff(a,x2),x2),diff(diff(a,x2),x3);...
    diff(diff(a,x3),x1),diff(diff(a,x3),x2),diff(diff(a,x3),x3)]

delta=[];
for i=1:3
    disp '  ----';
    i
    A(1:i,1:i)
    delta=det(A(1:i,1:i))
end

Golden Section

2005-09-28T17:59:00.000-07:00

clc;

clear;

f=@(x) x.^2+2*x;

a=-3;

b=3;

error=0.8;

r=(sqrt(5)-1)/2;

x0=a;

x3=b;

x1=a+(1-r)*(b-a);

x2=a+r*(b-a);

y1=feval(f,x1);

y2=feval(f,x2);

while (x3-x0)*r>error

if y1<=y2

x3=x2;

x2=x1;

x1=x0+(1-r)*(x3-x0);

else

x0=x1;

x1=x2;

x2=x0+r*(x3-x0);

end

y1=feval(f,x1);

y2=feval(f,x2);

end

y_array=[feval(f,x0),y1,y2,feval(f,x3)];

[ymin,index]=min(y_array);

x_array=[x0,x1,x2,x3];

xmin=x_array(index);

disp(sprintf('minimum of f(x) is %5.3f that occurs at x=%5.3f\n',ymin,xmin))

Secant Method

2005-09-28T17:58:00.000-07:00

clc;

clear;

% Shows the points obtained by the Secant Method when

% solving the equation func(x)=0. The two initial

% points are supplied, along with the number of

% required iterations.

func=@(x) 4*(x-2)^3+4*x-12;

current_1=0

current_2=10

iterations=5;

orbit=[current_1,current_2];

for i=1:iterations

disp('itter')

i

disp('slope')

disp('f(L)=')

feval(func,current_1)

disp('f(R)=')

feval(func,current_2)

disp('slope')

slope=(feval(func,current_2)-feval(func,current_1))/(current_2-current_1)

disp('new L')

current_1=current_2

disp('new R')

current_2=current_2-feval(func,current_2)/slope

orbit=[orbit,current_2];

end

orbit

Golden Section Vs Fminbnd Function

2005-09-28T17:54:00.000-07:00

clc;
clear;
close all;

L=1000; %feet

Q=20; %gmp


f=@(D) 0.45*L+(0.245*L*(D.^1.5))+325*((44*10^-8*(L*Q.^3/(D.^5))+1.92*(10^-9)*(L*Q^2.68)/(D^4.68)).^(1/2))+61.6*(44*10^-8*(L*Q.^3/(D.^5))+1.92*(10^-9)*(L*Q^2.68)/(D^4.68)).^0925+102

a=.25;
b=6;
error=0.0001; %Stoping Criteria for Golden Section


%=====================

%Start Golden Section

r=(sqrt(5)-1)/2;
x0=a;
x3=b;
x1=a+(1-r)*(b-a);
x2=a+r*(b-a);
y1=feval(f,x1);
y2=feval(f,x2);
i=0
itteration_goldensection=[];
error_goldensection=[];
while (x3-x0)*r>error
    if y1<=y2
        x3=x2;
        x2=x1;
        x1=x0+(1-r)*(x3-x0);
    else
        x0=x1;
        x1=x2;
        x2=x0+r*(x3-x0);
    end
    y1=feval(f,x1);
    y2=feval(f,x2);
    i=i+1;
    itteration_goldensection=[itteration_goldensection i];
    error_goldensection=[error_goldensection abs(feval(f,(x0+x3)/2)-0.425516)];
end
plot(itteration_goldensection,error_goldensection,'r-.');
y_array=[feval(f,x0),y1,y2,feval(f,x3)];
[ymin,index]=min(y_array);
x_array=[x0,x1,x2,x3];
xmin=x_array(index);
disp(sprintf('Golden Search: minimum of f(x) is %10.6f that occurs at x=%10.6f\n',ymin,xmin))
%====================

%MATLAB FUNCTION

[x, feval]=fminbnd(f,.25,6)

Bisection Search

2005-09-28T17:52:00.000-07:00

clc;
clear;

% Shows the points obtained by the Bisection Method when

% solving the equation func(x)=0. The two initial points

% are supplied, along with the number of required iterations.

% It is assumed that the initial interval is "good".

% Of course, in practice one should not set beforehand the

% number of iterations, but rather stop when reaching a

% sufficiently small neighbourhood of the root.

iterations=1000;
left=0;
right=10;
func=@(x) 4*(x-2)^3+4*x-12;

orbit=[];
for i=1:iterations
disp('======')
disp('itteration')
    i
    disp('left');
    left
    disp('right');
    right
    disp('feval left')
   
feval(func,left)
disp('feval right')
feval(func,right)

disp('middle')
    middle=(left+right)/2
    disp('feval middle')
feval(func,middle)
if feval(func,left)*feval(func,middle)>0 left=middle; disp('left*midd>0');
else right=middle; disp('left*midd<=0')
end
orbit=[orbit,middle];
end

orbit

Golden Section Minimization

2005-09-28T17:51:00.000-07:00

clc;
clear;

f=@(x) x.^2+2*x;
a=-3;
b=3;
error=0.8;


r=(sqrt(5)-1)/2;

x0=a;
x3=b;
x1=a+(1-r)*(b-a);
x2=a+r*(b-a);

y1=feval(f,x1);
y2=feval(f,x2);


while (x3-x0)*r>error
    if y1<=y2
        x3=x2;
        x2=x1;
        x1=x0+(1-r)*(x3-x0);
      
    else
        x0=x1;
        x1=x2;
        x2=x0+r*(x3-x0);
    end
   
y1=feval(f,x1);
y2=feval(f,x2);
end
y_array=[feval(f,x0),y1,y2,feval(f,x3)];
[ymin,index]=min(y_array);
x_array=[x0,x1,x2,x3];
xmin=x_array(index);

disp(sprintf('minimum of f(x) is %5.3f that occurs at x=%5.3f\n',ymin,xmin))

How to find minimum value of a function

2005-09-28T17:45:00.000-07:00

y=0:.001:6;
f=(sqrt(25+(6-y).^2)/3)+(y/5);
plot(y,f);
min(f)
a=y(find(f==min(f))) %returns the optimum value of y
x=sqrt(25+(6-a)^2)   %returns the optimum value of x
heading_angle=acos(5/x)*180/pi   %returns the optimum value of headding angle