%Mason Averill
%ME-480 Fall 2020, 11/20
%Optimizer




%%______________________________________________________________________
%input parameters

l=1; %length of the beam in meters
t=0.05; %width of the beam in meters
e=200*10^9; %Young's Modulus of the steel utilized in Pa
density=7900; %density of the steel utilized in kg/m^3
gravity=9.81; %acceleration due to gravity
Wmax=200; %maximum weight of beam in N
p=5000; %Point load applied at end of beam in N
number_of_elements=100; %specify number of elements to consider
step_size_const_solve=0.0001;

delta_x=l/number_of_elements;

x=0:delta_x:l;

%First lets find the deflection for a constant cross section beam

hmax=Wmax/(density*gravity*t*l); %maximum height permissible to still meet max weight criteria

deflection_constant=(p*l^3)/(3*e*t*(1/12)*hmax^3); %deflection for the constant cross-section beam




%_______________________________________________________________________________________________________________
%Now lets find the ideal linear case
syms x1 m b

h_x_linear=m*x1+b; %general equation for linear h(x)

int_h_x_linear=vpa(int(h_x_linear,x1,[0 l])); %integrating h(x) symbolically from 0 to l

m=solve(int_h_x_linear==Wmax/(density*gravity*t),m); %setting the integral from 0 to l for h(x) equal to the max weight to eliminate one coefficient
                                                     

h_x_linear=m*x1+b;%now have a function h(x) of only one unknown coefficient b

%find deflection in terms of b

deflection_linear=(12*p/(e*t))*vpaintegral((x1^2/h_x_linear^3),x1,[0 l]); %this is the deflection of the beam at the free end in terms of b 
                                                                                                                               
height_max=0.023; %set an upper bound for the unknown coefficient
b_value=0.018; %initialize b value variable 
deflection_b=[];
i=1;
b_store=[];
while(b_value<height_max)
    try
deflection_b(1,i)=double(subs(deflection_linear,b,b_value)); %adds the deflection corresponding to a particular value of b to this array


b_store(1,i)=b_value; %stores the value of b for each of the values in the above array


i=i+1;

    end
b_value=b_value+step_size_const_solve;

end 


[min_value,column]=min(deflection_b); %returns the minimum value of deflection at the free end and the location of b in b_store accompanying this value

b_value=b_store(1,column);%selects the optimum b from all considered b values
figure (1)
plot(b_store,deflection_b);%shows a plot of b values vs the deflection at the free end 
xlabel('b values')
ylabel('deflection at free end')
title('Linear Case: deflection at free end vs b')
result=min_value/deflection_constant; 

deflection_linear_store=min_value;


result=num2str(result);

string_out=strcat('The optimum linear case results in a deflection at the free end of the beam of  ',  result , ' times the deflection for the constant cross section case');

disp(string_out)

h_x_linear_sym=subs(h_x_linear,b,b_value);

h_x_linear=subs(h_x_linear,b,b_value);
h_x_linear=subs(h_x_linear,x1,x);
figure (2)
plot(x,h_x_linear)
xlabel('x')
ylabel('h(x)')
title('Linear Case: h(x) vs x')    
    
%_______________________________________________________________________________________________________________

%Now lets find the ideal polynomial order 2 case

syms x1 c1 c2

h_x_polynomial=c1*x1^2+c2; %general equation for polynomial order 2 h(x)

int_h_x_polynomial=vpa(int(h_x_polynomial,x1,[0 l])); %integrating h(x) symbolically from 0 to l

c1=solve(int_h_x_polynomial==Wmax/(density*gravity*t),c1); %setting the integral from 0 to l for h(x) equal to the max weight to eliminate one coefficient
                                                     

h_x_polynomial=c1*x1^2+c2;%now have a function h(x) of only one unknown coefficient c2

%find deflection in terms of c2

deflection_polynomial=12*p/(e*t)*vpaintegral((x1^2/h_x_polynomial^3),x1,[0.001 l]); %this is the deflection of the beam at the free end in terms of b 
                                                                                                                               
height_max=.035; %set an upper bound for the height at the end of the beam 
c2_value=.033; %initialize c2 value variable 
deflection_c2=[];
i=1;
c2_store=[];
while(c2_value<height_max)
try
deflection_c2(1,i)=double(subs(deflection_polynomial,c2,c2_value)); %adds the deflection corresponding to a particular value of c2 to this array


c2_store(1,i)=c2_value; %stores the value of c2 for each of the values in the above array


i=i+1;
end 

c2_value=c2_value+step_size_const_solve;

end 


[min_value,column]=min(deflection_c2); %returns the minimum value of deflection at the free end and the location of c2 in c2_store accompanying this value

c2_value=c2_store(1,column);%selects the optimum c2 from all considered c2 values
figure (3)
plot(c2_store,deflection_c2);%shows a plot of c2 values vs the deflection at the free end 
xlabel('c2 values')
ylabel('deflection at free end')
title('Polynomial order 2 Case: deflection at free end vs c2')
result=min_value/deflection_constant; 

result=num2str(result);

string_out=strcat('The optimum polynomial order 2 case results in a deflection at the free end of the beam of  ',  result , ' times the deflection for the constant cross section case');

disp(string_out)

h_x_polynomial_sym=subs(h_x_polynomial,c2,c2_value);

h_x_polynomial=subs(h_x_polynomial,c2,c2_value);

h_x_polynomial=subs(h_x_polynomial,x1,x);
figure (4)
plot(x,h_x_polynomial)
xlabel('x')
ylabel('h(x)')
title('Polynomial order 2 case: h(x) vs x')   

%_______________________________________________________________________________________________________________

%Now lets find the ideal polynomial order 4 case

syms x1 c1 c2

h_x_polynomial_o4=c1*x1^4+c2; %general equation for polynomial order 4 h(x)

int_h_x_polynomial_o4=vpa(int(h_x_polynomial_o4,x1,[0 l])); %integrating h(x) symbolically from 0 to l

c1=solve(int_h_x_polynomial_o4==Wmax/(density*gravity*t),c1); %setting the integral from 0 to l for h(x) equal to the max weight to eliminate one coefficient
                                                     

h_x_polynomial_o4=c1*x1^4+c2;%now have a function h(x) of only one unknown coefficient c2

%find deflection in terms of c2

deflection_polynomial_o4=12*p/(e*t)*vpaintegral((x1^2/h_x_polynomial_o4^3),x1,[0 l]); %this is the deflection of the beam at the free end in terms of b 
                                                                                                                               
height_max=0.05; %set an upper bound for the height at the end of the beam 
c2_value_o4=0.035; %initialize c2 value variable
deflection_c2_o4=[];
i=1;
c2_store_o4=[];
while(c2_value_o4<height_max)
try
deflection_c2_o4(1,i)=double(subs(deflection_polynomial_o4,c2,c2_value_o4)); %adds the deflection corresponding to a particular value of c2 to this array


c2_store_o4(1,i)=c2_value_o4; %stores the value of c2 for each of the values in the above array


i=i+1;
end

c2_value_o4=c2_value_o4+step_size_const_solve;

end 


[min_value,column]=min(deflection_c2_o4); %returns the minimum value of deflection at the free end and the location of c2 in c2_store accompanying this value

c2_value_o4=c2_store_o4(1,column);%selects the optimum c2 from all considered c2 values
figure (5)
plot(c2_store_o4,deflection_c2_o4);%shows a plot of c2 values vs the deflection at the free end 
xlabel('c2 values')
ylabel('deflection at free end')
title('Polynomial order 4 Case: deflection at free end vs c2')
result=min_value/deflection_constant; 

result=num2str(result);

string_out=strcat('The optimum polynomial order 4 case results in a deflection at the free end of the beam of  ',  result , ' times the deflection for the constant cross section case');

disp(string_out)

h_x_polynomial_o4_sym=subs(h_x_polynomial_o4,c2,c2_value_o4);

h_x_polynomial_o4=subs(h_x_polynomial_o4,c2,c2_value_o4);

h_x_polynomial_o4=subs(h_x_polynomial_o4,x1,x);
figure (6)
plot(x,h_x_polynomial_o4)
xlabel('x')
ylabel('h(x)')
title('Polynomial order 4 Case: h(x) vs x')


%_______________________________________________________________________________________________________________

%Now lets find the ideal polynomial order 6 case

syms x1 c1 c2

h_x_polynomial_o6=c1*x1^6+c2; %general equation for polynomial order 6 h(x)

int_h_x_polynomial_o6=vpa(int(h_x_polynomial_o6,x1,[0 l])); %integrating h(x) symbolically from 0 to l

c1=solve(int_h_x_polynomial_o6==Wmax/(density*gravity*t),c1); %setting the integral from 0 to l for h(x) equal to the max weight to eliminate one coefficient
                                                     

h_x_polynomial_o6=c1*x1^6+c2;%now have a function h(x) of only one unknown coefficient c2

%find deflection in terms of c2

deflection_polynomial_o6=12*p/(e*t)*vpaintegral((x1^2/h_x_polynomial_o6^3),x1,[0 l]); %this is the deflection of the beam at the free end in terms of b 
                                                                                                                               
height_max=0.0475; %set an upper bound for the height at the end of the beam
c2_value_o6=0.0445; %initialize c2 value variable 
deflection_c2_o6=[];
i=1;
c2_store_o6=[];
while(c2_value_o6<height_max)
try
deflection_c2_o6(1,i)=double(subs(deflection_polynomial_o6,c2,c2_value_o6)); %adds the deflection corresponding to a particular value of c2 to this array


c2_store_o6(1,i)=c2_value_o6; %stores the value of c2 for each of the values in the above array


i=i+1;

end
c2_value_o6=c2_value_o6+step_size_const_solve;

end 


[min_value,column]=min(deflection_c2_o6); %returns the minimum value of deflection at the free end and the location of c2 in c2_store accompanying this value

c2_value_o6=c2_store_o6(1,column);%selects the optimum c2 from all considered c2 values
figure (7)
plot(c2_store_o6,deflection_c2_o6);%shows a plot of c2 values vs the deflection at the free end 
xlabel('c2 values')
ylabel('deflection at free end')
title('Polynomial order 6 Case: deflection at free end vs c2')
result=min_value/deflection_constant; 

result=num2str(result);

string_out=strcat('The optimum polynomial order 6 case results in a deflection at the free end of the beam of  ',  result , ' times the deflection for the constant cross section case');

disp(string_out)

h_x_polynomial_o6_sym=subs(h_x_polynomial_o6,c2,c2_value_o6);

h_x_polynomial_o6=subs(h_x_polynomial_o6,c2,c2_value_o6);

h_x_polynomial_o6=subs(h_x_polynomial_o6,x1,x);
figure (8)
plot(x,h_x_polynomial_o6)
xlabel('x')
ylabel('h(x)')
title('Polynomial order 6 Case: h(x) vs x')

%_______________________________________________________________________________________________________________



%Now lets find the ideal polynomial order 12 case

syms x1 c1 c2

h_x_polynomial_o12=c1*x1^12+c2; %general equation for polynomial order 12 h(x)

int_h_x_polynomial_o12=vpa(int(h_x_polynomial_o12,x1,[0 l])); %integrating h(x) symbolically from 0 to l

c1=solve(int_h_x_polynomial_o12==Wmax/(density*gravity*t),c1); %setting the integral from 0 to l for h(x) equal to the max weight to eliminate one coefficient
                                                     

h_x_polynomial_o12=c1*x1^12+c2;%now have a function h(x) of only one unknown coefficient c2

%find deflection in terms of c2

deflection_polynomial_o12=12*p/(e*t)*vpaintegral((x1^2/h_x_polynomial_o12^3),x1,[0 l]); %this is the deflection of the beam at the free end in terms of b 
                                                                                                                               
height_max=0.051; %set an upper bound for the height at the end of the beam 
c2_value_o12=0.046; %initialize c2 value variable
deflection_c2_o12=[];
i=1;
c2_store_o12=[];
while(c2_value_o12<height_max)
try
deflection_c2_o12(1,i)=double(subs(deflection_polynomial_o12,c2,c2_value_o12)); %adds the deflection corresponding to a particular value of c2 to this array


c2_store_o12(1,i)=c2_value_o12; %stores the value of c2 for each of the values in the above array


i=i+1;
end

c2_value_o12=c2_value_o12+step_size_const_solve;

end 


[min_value,column]=min(deflection_c2_o12); %returns the minimum value of deflection at the free end and the location of c2 in c2_store accompanying this value

c2_value_o12=c2_store_o12(1,column);%selects the optimum c2 from all considered c2 values
figure (9)
plot(c2_store_o12,deflection_c2_o12);%shows a plot of c2 values vs the deflection at the free end 
xlabel('c2 values')
ylabel('deflection at free end')
title('Polynomial order 12 Case: deflection at free end vs c2')
result=min_value/deflection_constant; 

result=num2str(result);

string_out=strcat('The optimum polynomial order 12 case results in a deflection at the free end of the beam of  ',  result , ' times the deflection for the constant cross section case');

disp(string_out)

h_x_polynomial_o12_sym=subs(h_x_polynomial_o12,c2,c2_value_o12);

h_x_polynomial_o12=subs(h_x_polynomial_o12,c2,c2_value_o12);

h_x_polynomial_o12=subs(h_x_polynomial_o12,x1,x);
figure (10)
plot(x,h_x_polynomial_o12)
xlabel('x')
ylabel('h(x)')
title('Polynomial order 12 Case: h(x) vs x')

%_______________________________________________________________________________________________________________


%Now lets find the ideal exponential case

syms x1 c1 c2

h_x_exp=c1*exp(c2*x1); %general equation for exp h(x)

int_h_x_exp=vpa(int(h_x_exp,x1,[0 l])); %integrating h(x) symbolically from 0 to l

c1=solve(int_h_x_exp==Wmax/(density*gravity*t),c1); %setting the integral from 0 to l for h(x) equal to the max weight to eliminate one coefficient
                                                     

h_x_exp=c1*exp(c2*x1);%now have a function h(x) of only one unknown coefficient c2

%find deflection in terms of c2

deflection_exp=12*p/(e*t)*vpaintegral((x1^2/h_x_exp^3),x1,[0 l]); %this is the deflection of the beam at the free end in terms of b 
                                                                                                                               
height_max=1.14; %set an upper bound for the height at the end of the beam
c2_value_exp=1.115; %initialize c2 value variable
deflection_c2_exp=[];
i=1;
c2_store_exp=[];
while(c2_value_exp<height_max)
try
deflection_c2_exp(1,i)=double(subs(deflection_exp,c2,c2_value_exp)); %adds the deflection corresponding to a particular value of c2 to this array


c2_store_exp(1,i)=c2_value_exp; %stores the value of c2 for each of the values in the above array


i=i+1;
end

c2_value_exp=c2_value_exp+step_size_const_solve;

end 

[min_value,column]=min(deflection_c2_exp); %returns the minimum value of deflection at the free end and the location of c2 in c2_store accompanying this value

c2_value_exp=c2_store_exp(1,column);%selects the optimum c2 from all considered c2 values
figure (11)
plot(c2_store_exp,deflection_c2_exp);%shows a plot of c2 values vs the deflection at the free end 
xlabel('c2 values')
ylabel('deflection at free end')
title('Exponential Case: deflection at free end vs c2')
result=min_value/deflection_constant; 

result=num2str(result);

string_out=strcat('The optimum exponential case results in a deflection at the free end of the beam of  ',  result , ' times the deflection for the constant cross section case');

disp(string_out)

h_x_exp_sym=subs(h_x_exp,c2,c2_value_exp);

h_x_exp=subs(h_x_exp,c2,c2_value_exp);

h_x_exp=subs(h_x_exp,x1,x);
figure (12)
plot(x,h_x_exp)
xlabel('x')
ylabel('h(x)')
title('Exponential Case: h(x) vs x')



%_______________________________________________________________________________________________________________

%Now lets find the ideal cos case

%syms x1 c1 c2

%h_x_cos=c1*cos(c2*x1); %general equation for cos h(x)

%int_h_x_cos=vpa(int(h_x_cos,x1,[0 l])); %integrating h(x) symbolically from 0 to l

%c1=solve(int_h_x_cos==Wmax/(density*gravity*t),c1); %setting the integral from 0 to l for h(x) equal to the max weight to eliminate one coefficient
                                                     

%h_x_cos=c1*cos(c2*x1);%now have a function h(x) of only one unknown coefficient c2

%find deflection in terms of c2

%deflection_cos=12*p/(e*t)*vpaintegral((x1^2/h_x_cos^3),x1,[0 l]); %this is the deflection of the beam at the free end in terms of b 
                                                                                                                               
%height_max=0.04; %set an upper bound for the height at the end of the beam
%c2_value_cos=-0.1; %initialize c2 value variable
%deflection_c2_cos=[];
%i=1;
%c2_store_cos=[];
%while(c2_value_cos<height_max)
%try
%deflection_c2_cos(1,i)=double(subs(deflection_cos,c2,c2_value_cos)); %adds the deflection corresponding to a particular value of c2 to this array


%c2_store_cos(1,i)=c2_value_cos; %stores the value of c2 for each of the values in the above array


%i=i+1;
%end

%c2_value_cos=c2_value_cos+step_size_const_solve;

%end 

%[min_value,column]=min(deflection_c2_cos); %returns the minimum value of deflection at the free end and the location of c2 in c2_store accompanying this value

%c2_value_cos=c2_store_cos(1,column);%selects the optimum c2 from all considered c2 values
%figure (13)
%plot(c2_store_cos,deflection_c2_cos);%shows a plot of c2 values vs the deflection at the free end 
%xlabel('c2 values')
%ylabel('deflection at free end')
%title('Cosine Case: deflection at free end vs c2')
%result=min_value/deflection_constant; 

%result=num2str(result);

%string_out=strcat('The optimum cosine case results in a deflection at the free end of the beam of  ',  result , ' times the deflection for the constant cross section case');

%disp(string_out)

%h_x_cos_sym=subs(h_x_cos,c2,c2_value_cos);

%h_x_cos=subs(h_x_cos,c2,c2_value_cos);

%h_x_cos=subs(h_x_cos,x1,x);
%figure (14)
%plot(x,h_x_cos)
%xlabel('x')
%ylabel('h(x)')
%title('Cosine Case: h(x) vs x')