George Mason University

Conor Philip Nelson

Source: Home > Classes > Math 447 > Project 1 > Step 5 > 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 both outputs for both pictures:

format long
hsquares = zeros(11,1);
errors = zeros(11,1);
conditionA = zeros(11,1);
nlist = zeros(11,1);
scnd = zeros(11,1);
fnl = zeros(11,1);
for k = 1:11
for n = 10 * 2^k
n
nlist(k,1) = n;
j = n - 2;
L = 2; % Length of beam
h = L / n; % interval length for approximate solution
hsquare = h^2;
hsquares(k,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
x = h:h:L;
p = 100;
s = - p * g * sin( pi * x / L );
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;
conditionA(k,1) = condest(A);
sol = ( ( h^(4) / ( E * I ) ) * f' ); % Solution vector
y = A \ sol; % approximate solution
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 ) ) );
%figure(k+1)
if k == 7
figure(1)
plot(x,y,'ro',x,y1,'b')
title('Plot of Clamped-Free Beam with Sinusoidal Load for n = 1280')
xlabel('Distance (m)')
ylabel('Height (m)')
end
y2 = y1';
p1 = y2(n,1);
p2 = y(n,1);
err = p1 - p2;
errors(k,1) = err;
scnd(k,1) = y(n-1,1);
fnl(k,1) = y(n,1);
end
end
r = [nlist errors]
%q = [scnd fnl]
%conditionA
figure(2)
loglog(hsquares,errors,'bo');
title('Plot of Error against h for a Clamped-Free Diving Board with Sinusoidal Load')
xlabel('Size of h')
ylabel('Size of error')

Back to Step 5