04_linear_algebra
Homework #4
Code is known as HilbertHWP2.m The following are my answers:
ans 9.4366e+05
fnormH4 1.5097
fnormH5 1.5809
fnormiH4 1.0342e+04
fnormiH5 3.0416e+05
H4 4x4 double
H5 5x5 double
i 5
iH4 4x4 double
iH5 5x5 double
infnorm_H4 2.0833
infnorm_H5 2.2833
infnorm_iH4 140.0000
infnorm_iH5 630.0000
j 5
n1 4
n2 5
norm1_H4 2.0833
norm1_H5 2.2833
norm1_iH4 140.0000
norm1_iH5 630.0000
norm2_H4 1.5002
norm2_H5 1.5671
norm2_iH4 1.0341e+04
norm2_iH5 3.0414e+05
Problem #2
% HilbertHWP2 creates and modifies 4x4 and 5x5
% Hilbert matrix. These modifications are the calculations of certain
% norms and their condition numbers of both the regular and inverse Hilbert
% matrices.
n1=4;
n2=5;
%% Hilbert matricies
H4=zeros(n1,n1); %Create first Hilbert matrix
for i=1:n1
for j=1:n1
H4(i,j)=1/(i+j-1);
end
end
H5=zeros(n2,n2); %Create second Hilbert matrix
for i=1:n2
for j=1:n2
H5(i,j)=1/(i+j-1);
end
end
% Calculate the Frobenius Norms and Their Conditions
for i=1:n1
for j=1:n1
fnormH4(i,j)=H4(i,j)^2;
end
end
fnormH4=sqrt(sum(sum(fnormH4)));
cond(H4,'fro');
for i=1:n2
for j=1:n2
fnormH5(i,j)=H5(i,j)^2;
end
end
fnormH5=sqrt(sum(sum(fnormH5)));
cond(H5,'fro');
% Calculate the 1-Norms and Their Conditions
norm1_H4=[];
for j=1:n1
norm1_H4=max([norm1_H4,sum(H4(:,j))]);
end
cond(H4,1);
norm1_H5=[];
for j=1:n2
norm1_H5=max([norm1_H5,sum(H5(:,j))]);
end
cond(H5,1);
% Compute 2-Norms and Conditions
norm2_H4=max(sqrt(eig(H4'*H4)));
cond(H4,2);
norm2_H5=max(sqrt(eig(H5'*H5)));
cond(H5,2);
% Compute Infinty Norms and Conditions
infnorm_H4=[];
for i=1:n1
infnorm_H4=max([infnorm_H4,sum(H4(i,:))]);
end
cond(H4,inf);
infnorm_H5=[];
for i=1:n2
infnorm_H5=max([infnorm_H5,sum(H5(i,:))]);
end
cond(H5,inf);
%% Part 2b
% Compute for Inverse Hilbert Matricies
% Compute Frobenius Norms and Conditions
iH4=inv(H4);
iH5=inv(H5);
for i=1:n1
for j=1:n1
fnormiH4(i,j)=iH4(i,j)^2;
end
end
fnormiH4=sqrt(sum(sum(fnormiH4)));
cond(iH4,'fro');
for i=1:n2
for j=1:n2
fnormiH5(i,j)=iH5(i,j)^2;
end
end
fnormiH5=sqrt(sum(sum(fnormiH5)));
cond(iH5,'fro');
%% Compute 1-Norms and Conditions
norm1_iH4=[];
for j=1:n1
norm1_iH4=max([norm1_iH4,sum(iH4(:,j))]);
end
cond(iH4,1);
norm1_iH5=[];
for j=1:n2
norm1_iH5=max([norm1_iH5,sum(iH5(:,j))]);
end
cond(iH5,1);
%% Compute the 2-Norms and Their Conditions
norm2_iH4=max(sqrt(eig(iH4'*iH4)));
cond(iH4,2);
norm2_iH5=max(sqrt(eig(iH5'*iH5)));
cond(iH5,2);
%% Compute Infinty Norms and Conditions
infnorm_iH4=[];
for i=1:n1
infnorm_iH4=max([infnorm_iH4,sum(iH4(i,:))]);
end
cond(iH4,inf);
infnorm_iH5=[];
for i=1:n2
infnorm_iH5=max([infnorm_iH5,sum(iH5(i,:))]);
end
cond(iH5,inf);
%{
fprintf('Frobenius norm of 4x4 Hilbert matrix = %.3f\nFrobenius norm of 4x4 Hilbert matrix = %.3f\n',fnormH4, fnormH5)
%}
Problem #3
function [d,u]=chol_tridiag(e,f);
%This function chol.tridiag.m calculates the Cholevsky factorization of a
%tri-diagonal matrix
A=diag(e)+diag(f,1)+diag(f,-1); %initialize matrix to be factorized
[m,n]=size(A);
C=zeros(m,n); %create Cholesky matrix
%do Cholesky factorization
C(1,1)=sqrt(A(1,1));
for j=2:n
C(1,j)=A(1,j)/C(1,1);
end
for i=2:n
C(i,i)=sqrt(A(i,i)-sum(C(:,i).^2));
for j=i+1:n
C(i,j)=(A(i,j)-sum(C(:,i).*C(:,j)))/(C(i,i));
end
end
d=diag(C)
u=diag(C,1)
end
Problem #4
function x = solve_tridiag(d,u,b)
%This function solves Cx = b when given inputs of the diagonal (d) and the
%1 offset from upper diagonal (u) of the Cholesky matrix
d=diag(C); %change the diagonal of the Cholesky matrix to vector form
u=diag(C,1); %put the 1 upper diagonal of the Cholesky matrix into vector form
n=length(d);
b=ones(n,1); %create b matrix of Cx=b
v=zeros(n,1);
x=zeros(n,1);
x(1)=b(1)/d(1);
for i=2:n %start forward substitution
v(i-1)=u(i-1)/d(1);
d(1)=d(i);
x(i)=(b(i))/d(1);
end
for j=n-1:-1:1 %start back substitution
x(j)=x(j)-v(j)*x(j+1);
end
end
Problem #5
%stiffness.m calculates the stiffness matrix for four masses of pre-defined
%spring constants
%give values for spring constants
K1=1000;
K2=1000;
K3=1000;
K4=1000;
K=[K1+K2,-K2,0,0;-K2,K2+K3,-K3,0;0,-K3,K3+K4,-K4;0,0,-K4,K4]; %create stiffness matrix
cond(K) %compute the condition of the stiffness matrix
Results: 29.2841 is the condition of this matrix
Problem #6
%assign values for spring constants
K1=1000;
K2=1000;
K3=1000;
K4=1000;
%assign values for masses
m1=1;
m2=1;
m3=1;
m4=1;
g=9.81;
e=[K1+K2,K2+K3,K3+K4,K4];
f=[-K2,-K3,-K4];
b=g*[m1;m2;m3;m4];
[d,u]=chol_tridiag(e,f);
x=solve_tridiag(d,u,b)
d =
44.7214 38.7298 36.5148 15.8114
u =
-22.3607 -25.8199 -27.3861
x
x =
0.5907
0.7426
0.7340
0.6204
Problem #7
function [freq, modes] = mass_spring_vibrate(K1,K2,K3,K4)
K1=10;
K2=20;
K3=20;
K4=10;
m1=1; %input the masses
m2=2;
m3=4;
M=[m1,0,0;0,m2,0;0,0,m3]; %create mass matrix
K=[K1+K2,-K2,0;-K2,K2+K3,-K3;0,-K3,K3+K4]; %create stiffness matrix
[V,D]=eig(K,M); %calculate modes values and eigenvalues using eig command
for i=1:length(M)
freq(i)=sqrt(D(i,i)); %take sqrt of eigenvalues to get frequencies
end
modes=V;
end
ans =
1.6028 3.7959 6.3657
Problem #7
%This computes different modes of a beam given a particular nummber of
%segments and finds both the minimum and maximum loads in this beam
P=10; %force in N
E=76*10^9; %elastic modulus in Pa
I=4*10^-8; %moment of inertia in m^4
L=5; %length in m
%psqrd=P/(E*I);
N=5; %# of segmentsw
dx=L/N; %step size
A=zeros(N-1);
u=2*ones(1,N-1);
v=ones(1,N-2);
A=A+diag(u)+diag(v,-1)+diag(v,1); %create the coefficient matrix
[V,D]=eig(A); %find the eigenvalues and eigenvectors
plot(V(:,end)) %plot mode shapes (end=mode 1, end-1=mode 2, etc.)
P=E*I*(diag(D))./(dx^2); %solve for load vector
max_load=max(P); %find maximum load
min_load=min(P); %find minimum load
number_of_eigenvalues=length(diag(D));
%if the segment length approached zero there would be infinite eigenvalues
of segments | largest load (N) | smallest load (N) | # of eigenvalues |
| ------------ | ---------------- | ----------------- | ---------------- | | 5 | 1.099e4 | 1.1612e3 | 4 | | 6 | 1.6337e4 | 1.1730e3 | 5 | | 10 | 4.7450e4 |1.1903e3 | 9 |