George Mason University

Conor Philip Nelson

Source: Home > Classes > Math 447 > Project 1 > Step 2 > Code for Picture 1

Math 447: Numerical Analysis

Project 1

Step 1:

Step 2:

Step 3:

Step 4:

Step 5:

Step 6:

Step 7:

This is the code for Step 2.

n = 10; % number of intervals
k = n - 2;
L = 2; % Length of beam
h = L / n; % interval length for approximate solution
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
fc = -( rho * w * d * g );
f = -( rho * w * d * g ) * x ./ x; % weight of the beam on x_i
A = zeros(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:k
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
%x_2 = x.^2 - (4 * L * x) + (6 * L^(2))
%x_1 = ( fc / ( 24 * E * I ) ) * x.^(2)
y1 = ( fc / ( 24 * E * I ) ) * x.^(2) .* ( x.^(2) - (4 * L * x) + (6 * L^(2)) )
figure(1)
plot(x,y,'ro',x,y1);
err = y1' - y
figure(2)
plot(x, err);

Back to Step 2