George Mason University

Conor Philip Nelson

Source: Home > Classes > Math 447 > Project 1 > Step 6 > Code for Pictures 1&2

Math 447: Numerical Analysis

Project 1

Step 1:

Step 2:

Step 3:

Step 4:

Step 5:

Step 6:

Step 7:

This code produces out puts for both images.

n = 1280
j = n - 2;
L = 2; % Length of beam
h = L / n; % interval length for approximate solution
hsquare = h^2;
hsquares(n,1) = h;
w = .3; % width of the beam
d = .03; % thickness of the beam
g = 9.81;
rho = 480; % density of the beam
E = 1.3 * 10^(10); % Young's modulus of the wood
I = w * d^(3) / 12; % area moment of inertia
p = 100;
x = h:h:L;
s = x;

for i = 1:length(x)
if x(1, i)>=1.8
s(1, i)=-g*(70/0.2);
else
s(1, i) = 0;
end
end
fc = -( rho * w * d * g );
f = -( rho * w * d * g ) * x ./ x + s; % weight of the beam on x_i
A = sparse(n, n); % Commencement of construction of Coeff Matrix
A(1,1) = 16;
A(1,2) = -9;
A(1,3) = (8/3);
A(1,4) = (-1/4);
A(2,1) = -4;
A(2,2) = 6;
A(2,3) = -4;
A(2,4) = 1;
for i=3:j
A(i,i-2) = 1;
A(i,i-1) = -4;
A(i,i) = 6;
A(i,i+1) = -4;
A(i,i+2) = 1;
end
A(n-1,n-3) = (16/17);
A(n-1,n-2) = -(60/17);
A(n-1,n-1) = (72/17);
A(n-1,n) = -(28/17);
A(n,n-3) = -(12/17);
A(n,n-2) = (96/17);
A(n,n-1) = -(156/17);
A(n,n) = (72/17); % Completion of construction of Coeff Matrix
A;
sol = ( ( h^(4) / ( E * I ) ) * f' ); % Solution vector
y = A \ sol; % approximate solution
y_init = (fc/(24*E*I)).*x.^2.*(x.^2-4*L*x+6*L^2); %Step 2
y1 = ( fc / ( 24 * E * I ) ) .* x.^2 .* ( x.^2 - 4 * L * x + 6 * L^2 ) - ( ( p * g * L ) / ( E *I * pi ) ) * ( ( ( L^3 ) / ( pi^3 ) ) * sin( ( pi / L ) *x) - ( x.^3 / 6 ) + ( L / 2 ) * x.^2 - ( ( L^2 )*x / ( pi^2 ) ) );
%plot(x,y,'b',x,y_init,'r',x,y1,'g')

figure(1)
plot(x,y,'b')
title('Diver on the end of the board (Zoomed In)')
xlabel('distance(m)')
ylabel('height(m)')
axis([0 2 -0.3 0])

figure(2)
plot(x,y,'b')
title('Diver on the end of the board (Real Scale)')
xlabel('distance(m)')
ylabel('height(m)')
axis([0 2 -2 2])
y2 = y1';
p1 = y2(n,1);
p2 = y(n,1);
err = p1 - p2;
errors(n,1) = err;
deflection = p2;
deflection

Back to Step 6