Week 12 Challenge- MATLAB Module for solving Linear Equation and Calculating Shear force and Bending Moment Diagram
1) Develop a function in matlab to solve linear equations by the following methods
Ans :-
a) Cholesky method -
% solve_chol - solve linear equations from the Cholesky factorization.
% Solve A*X = B for X, where A is square, symmetric, positive definite. The
% input to the function is R the Cholesky decomposition of A and the matrix B.
% Example: X = solve_chol(chol(A),B);
%
% NOTE: The program code is written in the C language for efficiency and is
% contained in the file solve_chol.c, and should be compiled using matlabs mex
% facility. However, this file also contains a (less efficient) matlab
% implementation, supplied only as a help to people unfamiliar with mex. If
% the C code has been properly compiled and is avaiable, it automatically
% takes precendence over the matlab code in this file.
%
% Copyright (c) 2004, 2005, 2006 by Carl Edward Rasmussen. 2006-02-08.
function x = solve_chol(A, B);
if nargin ~= 2 | nargout > 1
error('Wrong number of arguments.');
end
if size(A,1) ~= size(A,2) | size(A,1) ~= size(B,1)
error('Wrong sizes of matrix arguments.');
end
x = A\(A'\B);b) Jacobi Iteration Method -
clc
clear
close all
x=input('enter the matrix with augmentation: ');
e=input ('enter the percentage error: ');
[r,c]=size(x); % get the number of rows and coloums in the matrix
error(1:r)=100;% create a rx1 matrix whixh contains the percentage relative error , intialized by the value 100 in each cell
maxerr=1000; % intialize the value of the maximum element of error matrix to 1000
xnew(1:r,1)=0; % intialize value of xnew to 0
while maxerr > e %repeat while maximum of all errors is greater than error entered by the user
xold=xnew; % assign new value of x to old value
for i = 1:r %loop the rows ,assign to xnew the value of the last element in the row
xnew(i,1) = x(i,c); %loop the coloums , if row not equalt to the coloums assign to xnew its value and substract xold multiplied by corresponding element in x
for j = 1 : c-1
if i~=j
xnew(i,1)= xnew(i,1)-x(i ,j)*xold(j,1); %else assign to k the elemnt inx with row number equal to coloumn number
else
k=x(i,j);
end
end
xnew(i,1)=xnew(i,1)/k; % divide xnew by k
end
for t= 1:r
error(t)=100*abs((xnew(t ,1)-xold(t,1))/xnew(t,1)); %loop xnew matrix calculate relative error for each x
end
maxerr = max(error); % get the maximum of all errors calculated
end
disp(['soloution is ' num2str(m)])2) Develop a function in matlab to compute the condensed matrix based on pre-selected degrees of freedom and degrees of freedom that desires to be eliminated
Ans :- As original equation is [k][delta]=[P]
now using condensation process we gonna used specified degree of freedom [delta(b)] and degree of freedom to be eliminated [delta(c)].
so delta = [delta(b) delta(c)] ,
[K] = [Kbb Kbc; Kcb Kcc] and
[P] = [Pb Pc]
therefore original equation becomes
[Kbb Kbc; Kcb Kcc] [delta(b) delta(c)] = [Pb Pc]
solving above we get equatoins,
Kbb*delta(b) + Kbc*delta(c) = Pb ....(1)
Kcb*delta(b) + Kcc*delta(c) = Pc ....(2)
reducing equation (1) into value of delta(b) we get,
delta(b) = [(Kbb^-1)*Pb] - [(Kbb^-1)*Kbc*delta(c)]
putting delta(b) value in equation (2),
K(condensed) = [Kcc - Kcb*Kb^-1*Kbc]
so we get condesnsed equation as above.
3) Develop a function in matlab to calculate critical load capacity of column based on their boundary conditions
a) Fixed-fixed
b) Pin-pinned
c) Fixed-Free
Ans :-
clc
clear all
close all
% Material Properties
E=100/pi()^2
%Dimensions
L=10
I=1
% Discretization
n_elem=200;
Le=L/n_elem;
% Matrix assembly
k0=zeros(n_elem+1,n_elem+1);
k1=zeros(n_elem+1,n_elem+1);
for i=1:n_elem
k0(i,i) = k0(i,i)+Le/3;
k0(i,i+1)=Le/6;
k0(i+1,i)=Le/6;
k0(i+1,i+1)=k0(i+1,i+1)+Le/3;
end
for i=1:n_elem
k1(i,i) = k1(i,i)+1/Le;
k1(i,i+1)=-1/Le;
k1(i+1,i)=-1/Le;
k1(i+1,i+1)=k1(i+1,i+1)+1/Le;
end
k0;
k1;
% Theoretical value of First critical buckling load
P_cr1=(pi())^2*E*I/L^2 % This value is 1.
i=0;
j=0;
P=P_cr1/10;
for index=1:1000
i=i+1;
P=P+P_cr1/10;
kdet= k0*P/(E*I)-k1;
% det_val=det(kdet);
P_val(i,1)=P;
det_val(i,1)=det(kdet);
% Alogrithm to find P values which give
% zero determinant value
% There is change of sign when curve
% P_val vs det_val passes through zero.
if i>1
det_prev=det_val(i-1,1);
det_current=det_val(i,1);
sign_change=det_prev*det_current;
if sign_change<0
j=j+1;
P_cr(j,1)=P_val(i-1,1);
end
end
end
plot(P_val,det_val,'-or','LineWidth',0.25)
title('P vs determinant','FontName','Arial Narrow','FontSize',24)
P_cr
4) Develop a function to plot bending moment diagram and shear force diagram for the following structure frames
a)
b)
Ans :-
% This Matlab code can be used for simply supported beam with single point
% load or uniformly distributed to find the
% * Support reaction
% * Maximum Bending Moment
% * Shear force diagram
% * Bending Moment daigram
clc
clear
close all
disp('Simply Supported Beam');
% Data input section
disp(' ');
L = input('Length of beam in meter = ');
disp(' ');
disp('Type 1 for point load, Type 2 for udl')
Type = input('Load case = ');
if Type == 1
disp(' ');
W = input('Load applied in kN = ');
disp(' ');
a = input('Location of Load from left end of the beam in meter = ');
c = L-a;
R1 = W*(L-a)/L; % Left Support Reaction.
R2 = W*a/L; % Right Support Reaction.
else
disp(' ');
W = input('Uniformly distributed load in kN/m = ');
disp(' ');
b = input('Length of udl in meter = ');
disp(' ');
cg = input('C.G of udl from left end of the beam in meter = ');
a = (cg-b/2);
c = L-a-b;
R1 = W*b*(b+2*c)/(2*L); % Left Support Reaction.
R2 = W*b*(b+2*a)/(2*L); % Right Support Reaction.
end
% Discretization of x axis.
n = 1000; % Number of discretization of x axis.
delta_x = L/n; % Increment for discretization of x axis.
x = (0:delta_x:L)'; % Generate column array for x-axis.
V = zeros(size(x, 1), 1); % Shear force function of x.
M = zeros(size(x, 1), 1); % Bending moment function of x.
% Data processing section
if Type == 1
for ii = 1:n+1
% First portion of the beam, 0 < x < b
V(ii) = R1;
M(ii) = R1*x(ii);
% Second portion of the beam, b < x < L
if x(ii) >= a
V(ii) = R1-W;
M(ii) = R1*x(ii)-W*(x(ii)-a);
end
end
x1 = a;
Mmax = W*a*(L-a)/L;
else
for ii = 1:n+1
% First portion of the beam, 0 < x < a
if x(ii) < a
V(ii) = R1;
M(ii) = R1*x(ii);
elseif a <= x(ii) && x(ii)< a+b
% Second portion of the beam, a < x < a+b
V(ii) = R1-W*(x(ii)-a);
M(ii) = R1*x(ii)-W*((x(ii)-a)^2)/2;
elseif x(ii) >= (a+b)
% Second portion of the beam, a+b < x < L
V(ii) = -R2;
M(ii) = R2*(L-x(ii));
end
end
x1 = a+b*(b+2*c)/(2*L);
Mmax = W*b*(b+2*c)*(4*a*L+2*b*c+b^2)/(8*L^2);
end
disp(' ');disp (['Left support Reaction' ' = ' num2str(R1) ' ' 'kN'])
disp(' ');disp (['Left support Reaction' ' = ' num2str(R2) ' ' 'kN'])
disp(' ');disp (['Maximum bending moment' ' = ' num2str(Mmax) ' ' 'kNm'])
figure
subplot(2,1,1);
plot(x, V, 'r','linewidth',1.5); % Grafica de las fuerzas cortantes.
grid
line([x(1) x(end)],[0 0],'Color','k');
line([0 0],[0 V(1)],'Color','r','linewidth',1.5);
line([x(end) x(end)],[0 V(end)],'Color','r','linewidth',1.5);
title('Shear Force Diagram','fontsize',16)
text(a/2,V(1),num2str(V(1)),'HorizontalAlignment','center','FontWeight','bold','fontsize',16)
text((L-c/2),V(end),num2str(V(end)),'HorizontalAlignment','center','FontWeight','bold','fontsize',16)
axis off
subplot(2,1,2);
plot(x, M, 'r','linewidth',1.5); % Grafica de momentos flectores;
grid
line([x(1) x(end)],[0 0],'Color','k');
line([x1 x1],[0 Mmax],'LineStyle','--','Color','b');
title('Bending Moment Diagram','fontsize',16)
text(x1+1/L,Mmax/2,num2str(roundn(Mmax,-2)),'HorizontalAlignment','center','FontWeight','bold','fontsize',16)
text(x1,0,[num2str(roundn(x1,-2)) ' m'],'HorizontalAlignment','center','FontWeight','bold','fontsize',16)
axis off5) Explain how this course helped in improving your concepts in analysis and coding concepts in Matlab.
Ans :- We get to know that how to write programme for solving problems using problem logic and how to apply mathematical solution into words of programme.
6) How easy was it for you to implement algorithms developed in class as a matlab script? How can it be made better for the future classes ?
Ans :- Good. Maybe more live lesson about it.