05_curve_fitting
Problem #2
Given:
a. y=0.3745+0.98644x+0.84564/x
b. y= 22.47-1.36x+0.28x^2
c. y=4.0046e^(-1.5x)+2.9213e^(-0.3x)+1.5647e^(-0.05x)
d. y=0.99sin(t)+0.5sin(3t)
The following vectors of best fit constants approach the coefficients of the given functions therefore showing that the aforementioned functions are the best fit lines. A Plot of each function is shown after the descriptive data
a. vector = [0.374496644295301;0.986442953020134;0.845637583892619] coefficient of determination: 0.980050796115443
b. vector = [22.466909090909102;1.362689393939397;0.276344696969697] coefficient of determination: 0.911190560385360
c. vector = [4.004638968589944;2.921297480316984;1.564688485802708] coefficient of determination: 0.946217462294966
d. vector = [0.997274968105860;0.502505402280354] coefficient of determination: 0.921860759509490
% This function fits given discrete data points to that curve which depends on the points themselves
% Part a
xa=[1 2 3 4 5]';
ya=[2.2 2.8 3.6 4.5 5.5]';
Z=[xa.^0,xa.^1,1./xa];
[aa,fxa,r2a]=least_squares(Z,ya);
ha=figure(1);
plot(xa,fxa,xa,ya,'o')
saveas(ha,'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=[xb.^0,-xb.^1,xb.^2];
[ab,fxb,r2b]=least_squares(Z,yb);
hb=figure(2);
plot(xb,fxb,xb,yb,'o')
% 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)];
[ac,fxc,r2c]=least_squares(Z,yc);
hc=figure(3);
plot(xc,fxc,xc,yc,'o')
% 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)];
[ad,fxd,r2d]=least_squares(Z,yd);
hd=figure(4);
plot(xd,fxd,xd,yd,'o')
function [a,fx,r2]=least_squares(Z,y)
% 'least_squares.m` accepts a Z-matrix and dependent variable y as
% input and returns the vector of best-fit constants, a, the
% best-fit function evaluated at each point $f(x_{i})$, and the
% coefficient of determination, r2.
a=Z\y; %or Z'*Z\Z'*ya
fx=Z*a;
St=std(y); %=sum(y)-(mean(y))^2
Sr=std(y-fx);
r2=1-Sr/St;
end
Problem #3
a. Most accurate dart thrower: 31
b. Most precise dart thrower: 32
% This function determines both the most accurate and the most precise of the dart throwers
data=[0 1.49 3.2
0 6.21 45.4
0 11.56 340.1
0 4.12 273.4
0 7.48 205.7
0 9.43 133.2
0 8.11 159.5
0 5.59 21.4
0 4.91 82.3
1 10.74 115.8
1 3.7 295
1 2.4 300
1 4.2 134
1 0.6 236
1 0.5 18
1 3.3 16
1 7.6 224
1 6.1 290
1 2.3 224
2 4.9 8
2 7.5 74
2 4.5 56
2 5.25 102
2 6 92
2 2 134
2 4.5 16
2 4 66
2 13.25 284
2 3 84
3 7 290
3 1.5 18
3 6 63
3 5 115
3 4.5 122
3 5.5 137
3 0.5 205
3 7 76
3 2.5 117
3 4 101
4 2.5 281
4 3 340
4 5.5 185
4 8 67
4 3.5 7
4 2.5 255
4 3 290
4 3.5 295
4 4.5 250
4 8.5 275
5 10.5 304
5 1.2 120
5 16.6 90
5 14.3 152
5 14.4 150
5 4.6 210
5 5.9 100
5 3.7 86
5 14.2 10
5 7.3 65
6 10.1 88
6 1 30
6 4.5 60
6 6 101
6 6.1 110
6 1.4 160
6 0.8 200
6 4 250
6 3.7 273
6 2.3 320
7 1.4 340
7 11 260
7 12 53
7 12.5 88
7 8.2 138
7 3.6 130
7 2 15
7 6.2 175
7 9.5 225
7 12.5 302
8 4 268
8 10 98
8 4.6 83
8 2.7 102
8 3.4 42
8 2.8 24
8 3.4 -6
8 3.5 -42
8 2.9 -46
8 3.8 -47
9 5.6 -84
9 4 310
9 4 140
9 8 155
9 9 75
9 12 85
9 13 65
9 15 110
9 11 200
9 14 190
10 8 290
10 8 140
10 2.2 15
10 9.7 97
10 6 122
10 7.5 330
10 5.9 100
10 3.1 190
10 9.8 150
10 10.2 45
11 4.8 22
11 10 92
11 5.1 92
11 5.4 122
11 3.1 164
11 2.6 338
11 3.4 234
11 5.7 295
11 7.9 263
11 8.9 259
12 12.4 242
12 8.2 85
12 5.3 92
12 3.1 105
12 2.6 11
12 6.4 65
12 1.1 112
12 3.5 121
12 4.8 92
12 3.2 76
13 2 90
13 6.6 216
13 2.4 272
13 3 320
13 0.5 16
13 4.3 44
13 2.6 60
13 1.1 124
13 2.4 358
13 6.1 86
14 3.4 174
14 14 128
14 5 220
14 12.5 338
14 15.5 290
14 13.5 284
14 15 236
14 4 290
14 3 25
14 6.5 85
15 16.5 105
15 2 50
15 2.7 164
15 4 194
15 5 330
15 5.2 262
15 6 90
15 8.1 305
15 8.5 48
15 8.5 136
16 10 240
16 3 275
16 7 61
16 2 51
16 6 255
16 12 64
16 4 77
16 4 184
16 5 111
16 8 119
17 11 15
17 3.2 0
17 4.1 12
17 2 20
17 2.2 38
17 3.2 39
17 1.5 50
17 2.9 52
17 9.4 88
17 3.9 180
18 10.2 241
18 5.5 20
18 2 200
18 7 170
18 4.5 70
18 6.5 330
18 4 290
18 1 240
18 1.5 90
18 6 120
19 4.5 280
19 2.4 15
19 6.7 45
19 3.5 121
19 6.9 126
19 4.5 129
19 3.8 134
19 4.3 221
19 5 257
19 2.4 322
20 1.8 337
20 7.75 9
20 7.6 82.5
20 5.4 122
20 5.8 154.5
20 7.2 211
20 4.5 220
20 3.6 260
20 5 278
20 3.8 316
21 6.5 339
21 10 86
21 6.5 102
21 8 124
21 13 78
21 3.5 270
21 4 180
21 17 251
21 9 93
21 4.5 90
22 6 270
22 5 8.5
23 2.5 317
23 7 35
23 9 166
23 2 112
23 14 330
23 13 240
23 5 105
23 11 60
23 4 52
23 8 262
24 13 155
24 8.9 9
24 5 82
24 5.4 290
24 5.8 20
24 7.2 220
24 4.5 180
24 3.6 270
24 5 278
24 3.8 316
25 6.5 339
25 4 95
25 12 90
25 11.5 67
25 12.5 35
25 11 30
25 6 357
25 13 300
25 18 308
25 10 182
26 13 235
26 3 155
26 1.5 50
26 1.6 172
26 0.75 80
26 2.5 12.5
26 1.5 145
26 1.75 95
26 2.25 170
26 1.75 110
27 0.75 140
27 3.5 108
27 8 252
27 12.5 240
27 14 194
27 1.5 86
27 8 124
27 10.5 64
27 8 46
27 13 0
28 12.5 62
28 4.8 37
28 2.2 44
28 2.5 96
28 4.6 104
28 4.1 115
28 2.3 135
28 3.8 140
28 2.9 197
28 5.1 244
29 2.3 287
29 7.2 276
29 13.5 284
29 10.2 272
29 2.5 72
29 12.8 258
29 6.9 236
29 7.2 70
29 8.5 112
29 9 294
30 13.5 256
30 4 24
30 8 257
30 2 75
30 5 268
30 8 142
30 6 200
30 3 340
30 5 16
30 6 39
31 4 127
31 7 150
31 1 30
31 2.5 45
31 3 62
31 2.8 110
31 4 246
31 1.7 155
31 3.6 176
31 0.5 18
32 3.3 335
32 4.7 286];
mean_x=zeros(1,33);
mean_y=zeros(1,33);
mean_d=zeros(1,33);
std_x=zeros(1,33);
std_y=zeros(1,33);
std_d=zeros(1,33);
for i=1:33
j=find(data(:,1)==i-1);
mean_x(i)=mean(data(j,2).*sind(data(j,3)));
mean_y(i)=mean(data(j,2).*cosd(data(j,3)));
mean_d(i)=mean_x(i)+mean_y(i);
std_x(i)=std(data(j,2).*sind(data(j,3)));
std_y(i)=std(data(j,2).*cosd(data(j,3)));
std_d(i)=std_x(i)+std_y(i);
end
% The following line ascertains the identity of the most accurate thrower
most_accurate_thrower=find(mean_d==min(abs((mean_d))))-1;
% The following line ascertains the identity of the most precise thrower
most_precise_thrower=find(std_d==min(abs((std_d))))-1;
Problem #4
a. Mean buckle load: 171.2036, Standard buckle load: 72.5282
b. Length so that beams are only 2.5% will reach the critical buckling load: 4.8923
n=100; %% number of samples
mass=1000; % (kg)
E=200e9; % Elactic modulus (Pa)
r=0.01; % radius of rod (m)
L=5; % length of rod (m)
r_std=0.001; % standard deviation of radius
L_std=L*0.01; % standard deviation of length
[mean_buckle_load,std_buckle_load]=monte_carlo_buckle(E,r,L,r_std,L_std)
function [mean_buckle_load,std_buckle_load]=monte_carlo_buckle(E,r,L,r_std,L_std)
n=100; %%number of samples
mass=1000; %(kg)
r_rand=normrnd(r,r_std,[n,1]); % create normal distribution of r values
L_rand=normrnd(L,L_std,[n,1]); % create normal distribution of L values
mean_buckle_load=mean((pi^3*E.*r_rand.^4)./(16.*L_rand.^2));
std_buckle_load=std((pi^3*E.*r_rand.^4)./(16.*L_rand.^2));
L=((pi^3*E.*r_rand.^4)./(16*mean_buckle_load)).^0.5;
mean(L)
end
Problem #5
a. Linear: [Cd_out] = 0.0900 Piecewise Cubic: [Cd_out] = 0.0905 Continuous Cubic Spline: [Cd_out] = 0.0960
b.
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];
data(:,1)=(1e-4)*data(:,1);
Re_in=70*1e-4; %Input Reynold's Number
prompt='Input Spline Type: ';
spline_type=input(prompt,'s');
[Cd_out]=sphere_drag(Re_in,spline_type)
mu=1.78e-5;
V=linspace(4,40,1000);
D=23.5e-2; %Diameter in meters
rho=1.3; %density (kg/m^3)
Re=(rho*V*D/mu)*1e-8;
Cd_in=interp1(data(:,1),data(:,2),Re,'linear');
Cd_in_pc=interp1(data(:,1),data(:,2),Re,'pchip');
Cd_in_spline=interp1(data(:,1),data(:,2),Re,'spline');
plot(Re,Cd_in,data(1:end-4,1),data(1:end-4,2))
hold on
plot(Re,Cd_in_pc,data(1:end-4,1),data(1:end-4,2))
hold on
plot(Re,Cd_in_spline,data(1:end-4,1),data(1:end-4,2))
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];
data(:,1)=(1e-4)*data(:,1);
% x=(data(:,1),data(:,2),linspace(min(data(:,1)));
% Cd_in=interp1(data(:,1),data(:,2),linspace(min(data(:,1)),max(data(:,1))),spline_type);
%plot(linspace(min(data(:,1)),max(data(:,1))),Cd_in)
Cd_out=interp1(data(:,1),data(:,2),Re_in,spline_type);
end
Problem #6
The following are the approximate function values for the integrand which was given in the description for the sixth part of the fifth homework assignment. Their errors will be listed in a table which will append the following list:
a. In an analytical way: 8.375000000000000
b. By the means of the 1-point Gauss quadrature: 8.229166666666666
c. By the means of the 2-point Gauss quadrature: 8.375000000000000
d. By the means of the 3-point Gauss quadrature: 8.374999999999998
| method | value | error |
|---|---|---|
| analytical | 8.375000000000000 | 0% |
| 1 Gauss point | 8.229166666666666 | 1.741293532338315% |
| 2 Gauss point | 8.375000000000000 | 0% |
| 3 Gauss point | 8.374999999999998 | 2.121023091821194e-14% |
f_x=@(x) (1/6).*(x).^3+(1/2).*(x).^2+(x); %integrand
F_x=integral(f_x,2,3); %analytical solution to the integral
% Part a: 1 Point Gauss Quadrature
x=0;
w=2;
a=2;
b=3;
C=(b-a)/2*w;
X=(b-a)/2*x+(b+a)/2;
f=(1/6).*(X).^3+(1/2).*(X).^2+(X);
F_x1=sum(C.*f);
e1=(abs((F_x1-F_x)/F_x))*100;
%Part b: 2 Point Gauss Quadrature
x=[-1/sqrt(3),1/sqrt(3)];
w=[1,1];
C=(b-a)/2*w;
X=(b-a)/2*x+(b+a)/2;
f=(1/6).*(X).^3+(1/2).*(X).^2+(X);
F_x2=sum(C.*f);
e2=(abs((F_x2-F_x)/F_x))*100;
%Part c: 3 Point Gauss Quadrature
x=[-sqrt(3/5),0,sqrt(3/5)];
w=[5/9,8/9,5/9];
C=(b-a)/2*w;
X=(b-a)/2*x+(b+a)/2;
f=(1/6).*(X).^3+(1/2).*(X).^2+(X);
F_x3=sum(C.*f);
e3=(abs((F_x3-F_x)/F_x))*100;