05_Curve_Fitting
ME 3255 Homework 5
Problem 2
##Least Squares Function
function [a,fx,r2] = least_squares_2(Z,y)
%least squares solution for input matrix Z
%and measurements y
%outputs are:
%a: constants for best fit function
%fx: function evaluations at xi
%r2: coefficient of determination
a = Z\y;
fx = Z*a;
e=y-fx;
St=std(y);
Sr=std(e);
r2 = 1-Sr/St;
end
set defaults
Part A
xa=[1 2 3 4 5]'
ya=[2.2 2.8 3.6 4.5 5.5]';
Z=[ones(size(xa)) xa xa.^-1];
[a_a,fx_a,r2_a] = least_squares_2 (Z,ya)
xa_fcn = linspace(min(xa),max(xa));
a = figure(1);
plot(xa,ya,'o',xa_fcn,a_a(1)+a_a(2)*xa_fcn+a_a(3)*xa_fcn.^-1)
saveas(a, 'figure_01.png')
Part B
xb=[0 2 4 6 8 10 12 14 16 18]';
yb=[21.5 20.84 23.19 22.69 30.27 40.11 43.31 54.79 70.88 89.48]';
Z=[ones(size(xb)) xb xb.^2 xb.^3];
[a_b,fx_b,r2_b] = least_squares_2 (Z,yb)
xb_fcn = linspace(min(xb),max(xb));
b = figure(2);
plot(xb,yb,'o',xb_fcn,a_b(1)+a_b(2)*xb_fcn.^1+a_b(3)*xb_fcn.^2+a_b(4)*xb_fcn.^3)
saveas(b, 'figure_02.png')
##Part C
xc=[0.5 1 2 3 4 5 6 7 9]';
yc=[6 4.4 3.2 2.7 2.2 1.9 1.7 1.4 1.1]';
Z=[exp(-1.5*xc), exp(-0.3*xc), exp(-0.05*xc)];
[a_c,fx_c,r2_c] = least_squares_2 (Z,yc)
xc_fcn = linspace(min(xc),max(xc));
c = figure(3);
plot(xc,yc,'o',xc_fcn,a_c(1)*exp(-1.5.*xc_fcn) + a_c(2)*exp(-0.3.*xc_fcn) + a_c(3)*exp(-0.05.*xc_fcn));
saveas(c, 'figure_03.png')
Part D
xd=[0.00000000e+00 1.26933037e-01 2.53866073e-01 3.80799110e-01 5.07732146e-01 6.34665183e-01 7.61598219e-01 8.88531256e-01 1.01546429e+00 1.14239733e+00 1.26933037e+00 1.39626340e+00 1.52319644e+00 1.65012947e+00 1.77706251e+00 1.90399555e+00 2.03092858e+00 2.15786162e+00 2.28479466e+00 2.41172769e+00 2.53866073e+00 2.66559377e+00 2.79252680e+00 2.91945984e+00 3.04639288e+00 3.17332591e+00 3.30025895e+00 3.42719199e+00 3.55412502e+00 3.68105806e+00 3.80799110e+00 3.93492413e+00 4.06185717e+00 4.18879020e+00 4.31572324e+00 4.44265628e+00 4.56958931e+00 4.69652235e+00 4.82345539e+00 4.95038842e+00 5.07732146e+00 5.20425450e+00 5.33118753e+00 5.45812057e+00 5.58505361e+00 5.71198664e+00 5.83891968e+00 5.96585272e+00 6.09278575e+00 6.21971879e+00 6.34665183e+00 6.47358486e+00 6.60051790e+00 6.72745093e+00 6.85438397e+00 6.98131701e+00 7.10825004e+00 7.23518308e+00 7.36211612e+00 7.48904915e+00 7.61598219e+00 7.74291523e+00 7.86984826e+00 7.99678130e+00 8.12371434e+00 8.25064737e+00 8.37758041e+00 8.50451345e+00 8.63144648e+00 8.75837952e+00 8.88531256e+00 9.01224559e+00 9.13917863e+00 9.26611167e+00 9.39304470e+00 9.51997774e+00 9.64691077e+00 9.77384381e+00 9.90077685e+00 1.00277099e+01 1.01546429e+01 1.02815760e+01 1.04085090e+01 1.05354420e+01 1.06623751e+01 1.07893081e+01 1.09162411e+01 1.10431742e+01 1.11701072e+01 1.12970402e+01 1.14239733e+01 1.15509063e+01 1.16778394e+01 1.18047724e+01 1.19317054e+01 1.20586385e+01 1.21855715e+01 1.23125045e+01 1.24394376e+01 1.25663706e+01]';
yd=[9.15756288e-02 3.39393873e-01 6.28875306e-01 7.67713096e-01 1.05094584e+00 9.70887288e-01 9.84265740e-01 1.02589034e+00 8.53218113e-01 6.90197665e-01 5.51277193e-01 5.01564914e-01 5.25455797e-01 5.87052838e-01 5.41394658e-01 7.12365594e-01 8.14839678e-01 9.80181855e-01 9.44430709e-01 1.06728057e+00 1.15166322e+00 8.99464065e-01 7.77225453e-01 5.92618124e-01 3.08822183e-01 -1.07884730e-03 -3.46563271e-01 -5.64836023e-01 -8.11931510e-01 -1.05925186e+00 -1.13323611e+00 -1.11986890e+00 -8.88336727e-01 -9.54113139e-01 -6.81378679e-01 -6.02369117e-01 -4.78684439e-01 -5.88160325e-01 -4.93580777e-01 -5.68747320e-01 -7.51641934e-01 -8.14672884e-01 -9.53191554e-01 -9.55337518e-01 -9.85995556e-01 -9.63373597e-01 -1.01511061e+00 -7.56467517e-01 -4.17379564e-01 -1.22340361e-01 2.16273929e-01 5.16909714e-01 7.77031694e-01 1.00653798e+00 9.35718089e-01 1.00660116e+00 1.11177057e+00 9.85485116e-01 8.54344900e-01 6.26444042e-01 6.28124048e-01 4.27764254e-01 5.93991751e-01 4.79248018e-01 7.17522492e-01 7.35927848e-01 9.08802925e-01 9.38646871e-01 1.13125860e+00 1.07247935e+00 1.05198782e+00 9.41647332e-01 6.98801244e-01 4.03193543e-01 1.37009682e-01 -1.43203880e-01 -4.64369445e-01 -6.94978252e-01 -1.03483196e+00 -1.10261288e+00 -1.12892727e+00 -1.03902484e+00 -8.53573083e-01 -7.01815315e-01 -6.84745997e-01 -6.14189417e-01 -4.70090797e-01 -5.95052432e-01 -5.96497000e-01 -5.66861911e-01 -7.18239679e-01 -9.52873043e-01 -9.37512847e-01 -1.15782985e+00 -1.03858206e+00 -1.03182712e+00 -8.45121554e-01 -5.61821980e-01 -2.83427014e-01 -8.27056140e-02]';
Z = [sin(xd), sin(3*xd)];
[a_d,fx_d,r2_d] = least_squares_2 (Z,yd)
xd_fcn = linspace(min(xd),max(xd));
d = figure(4);
plot(xd,yd,'o',xd_fcn,a_d(1)*sin(xd_fcn)+a_d(2)*sin(3*xd_fcn));
saveas(d, 'figure_04.png')
##Output:
Part A: r2 = 0.9801
Part B: r2 = 0.9232
Part C: r2 = 0.9462
Part D: r2 = 0.9219
Part A
Part B
Part C
Part D
Problem 3
darts = dlmread('compiled_data.csv',',',1,0);
x_darts =darts(:,2).*cosd(darts(:,3));
y_darts =darts(:,2).*sind(darts(:,3));
ur = darts(:,1);
i = 1;
for r = 0:32
i_interest = find(ur==r);
mx_ur(i)= mean(x_darts(i_interest));
my_ur(i)= mean(y_darts(i_interest));
i = i +1;
end
accuracy = mx_ur + my_ur;
val = 0;
abs_accuracy = abs(accuracy-val);
[~, index] = min(abs_accuracy)
closest_value = accuracy(index)
darts = dlmread('compiled_data.csv',',',1,0);
x_darts =darts(:,2).*cosd(darts(:,3));
y_darts =darts(:,2).*sind(darts(:,3));
ur = darts(:,1);
i = 1;
for r = 0:32
i_interest = find(ur==r);
stdx_ur(i)= std(x_darts(i_interest));
stdy_ur(i)= std(y_darts(i_interest));
i = i +1;
end
precision = stdx_ur + stdy_ur;
val = 0;
abs_accuracy = abs(precision-val);
[~, index] = min(abs_accuracy)
closest_value = precision(index)
Part A
User number 32 was closest with a value of 0.0044.
Part B
User number 33 was the most precise with a standard deviation of 3.4073Problem 4
Part A
```matlab
function [mean_buckle_load,std_buckle_load]=buckle_monte_carlo(E,r_mean,r_std,L_mean,L_std);
N=100;
r_mean=0.01;
r_std=.001;
L_mean=5;
L_std=L_mean*0.01;
rrand=normrnd(r_mean,r_std,[N,1]);
Lrand=normrnd(L_mean,L_std,[N,1]);
E = 200*10^9;
P_cr=(pi^3*E.*rrand.^4)./(16.*Lrand.^2);
x1 = mean(P_cr)
x2 = std(P_cr)
P_test=1000*9.81/100;
sum(P_cr<P_test)
end
```matlab
Mean Buckle Load = 164.4154 N
Standard Deviation = 63.6238 N
Part B
```matlab
function [mean_buckle_load,std_buckle_load]=buckle_monte_carlo_2(E,r_mean,r_std,P_cr);
N=100;
r_mean=0.01;
r_std=.001;
P_test=1000*9.81/100;
sum(2.5<P_test)
P_cr = 164.4154;
rrand=normrnd(r_mean,r_std,[N,1]);
E = 200*10^9;
L = ((pi^3*E.*rrand.^4)./(16*P_cr)).^0.5;
mean(L)
end
Output
```matlab
L = 5.0108 m
#Problem 5
Part A
```matlab
function [Cd_out] = sphere_drag(Re_in,spline_type)
data = [2 0.52
5.8 0.52
16.8 0.52
27.2 0.5
29.9 0.49
33.9 0.44
36.3 0.18
40 0.074
46 0.067
60 0.08
100 0.12
200 0.16
400 0.19];
Cd_out=interp1(data(:,1),data(:,2),Re_in,spline_type);
rho=1.3;V=linspace(4,40);D=23.5e-2;mu=1.78e-5;
Re=rho*V*D/mu*1e-4;
end
```matlab
Cd_out1 = sphere_drag(150,'linear')
Cd_out2 = sphere_drag(150,'spline')
Cd_out3 = sphere_drag(150,'pchip')
Output
```matlab
Linear = 0.1400
Spline = 0.1428
Pchip = 0.1447
Part B Code
```matlab
function [Cd_out] = sphere_drag(Re_in,spline_type)
data = [2 0.52
5.8 0.52
16.8 0.52
27.2 0.5
29.9 0.49
33.9 0.44
36.3 0.18
40 0.074
46 0.067
60 0.08
100 0.12
200 0.16
400 0.19];
Cd_out=interp1(data(:,1),data(:,2),Re_in,spline_type)
Re=linspace(data(1,1),data(end,1),1000);
rho=1.3;V=linspace(4,40);D=23.5e-2;mu=1.78e-5;
Re=rho*V*D/mu*1e-4;
end
Part B Script
```matlab
Re=linspace(data(1,1),data(end,1),1000);
rho=1.3;V=linspace(4,40);D=23.5e-2;mu=1.78e-5;
Re=rho*V*D/mu*1e-4
a=figure(1);
cd_interp=sphere_drag(linspace(data(1,1),data(end,1),100),'linear');
plot(data(:,1),data(:,2),'o',Re,cd_interp,'b-','Linewidth',4,'Markersize',10)
hold on;
cd_interp=sphere_drag(linspace(data(1,1),data(end,1),100),'spline');
plot(Re,cd_interp,'r-', 'Linewidth',4,'Markersize',10)
cd_interp=sphere_drag(linspace(data(1,1),data(end,1),100),'pchip');
plot(Re,cd_interp,'g-', 'Linewidth',4,'Markersize',10)
legend('data','linear','spline','pchip')
saveas(a,'figure_14.png')
Script
```matlab
Cd_out1 = sphere_drag(10,'linear')
Cd_out2 = sphere_drag(10,'pchip')
Cd_out3 = sphere_drag(10,'spline')
plot(data(:,1),data(:,2),'o',Re,Cd_out,'b-','Linewidth',4,'Markersize',10)
plot(data(:,1),data(:,2),Re,Cd_out2,'r-', 'Linewidth',4,'Markersize',10)
plot(data(:,1),data(:,2),Re,Cd_out3,'g-', 'Linewidth',4,'Markersize',10)
legend('data','linear','pchip','spline')
Problem 6
##Analytical Solution
```matlab
fun = @(x) (1/6).*(x).^3+(1/2).*(x).^2+(x);
q = integral(fun,2,3)
##Gauss Quadtrature
```matlab
%One Point
x = [0];
c = [2];
a = 2;
b = 3;
tildec = (b-a)/2*c;
tildex = (b-a)/2*x + (b+a)/2;
f = (1/6).*(tildex).^3+(1/2).*(tildex).^2+(tildex);
value1 = sum(tildec.*f)
##Two Point
%Two Point
x = [-0.57735,0.57735];
c = [1,1];
a = 2;
b = 3;
tildec = (b-a)/2*c;
tildex = (b-a)/2*x + (b+a)/2;
f = (1/6).*(tildex).^3+(1/2).*(tildex).^2+(tildex);
value2 = sum(tildec.*f)
%Three Point
x = [-0.77459666,0,0.77459666];
c = [0.5555555,0.88888888,0.55555555];
a = 2;
b = 3;
tildec = (b-a)/2*c;
tildex = (b-a)/2*x + (b+a)/2;
f = (1/6).*(tildex).^3+(1/2).*(tildex).^2+(tildex);
value3 = sum(tildec.*f)
matlab
| Method | Value | Error |
| Analytical | 8.3750 | 0% |
| 1 Gauss Point | 8.2292 | 1.74% |
| 2 Gauss Point | 8.3750 | 0% |
| 3 Gauss Point | 8.3750 | 0% |