i have this code in finite element method to find basis but i have some errores

Question

i have this code in finite element method to find basis but i have some errores and the plot is missing,could you please fix it for me?

clear;

syms x; %Creates x as a variable

f = 0; %Function which is the u"

of = 0; %Original function

n = 8; %Number of Nodes

a = 0; %Lower Bound

b = 1; %Upper Bound

boundary = [1,n]; %Which nodes are 0;

alpha = 0; %Lower boundary value

beta = 1; %Upper boundary value

orderpairs = zeros(2*n,2);

localstiffness = zeros(2,2); %local stiffness matrix for individual elements

xnodevalues = zeros(n,1); %Create a column Vector of the Values at Xsubi

globalstiffness = zeros(n,n); %Create the Local Stiffness Matrices

globalload = zeros(n,1); %Create the global load matrix

localload = zeros(2,1); %Create a Local load vector

basis = zeros((n-1),2,2); %Creates the Shape Functions for integration as

%First column slope and second column y-int

%Builds the Nodal X Values

for i=1:n

xnodevalues(i,1) = a + ((i-1)*(b-a)/(n-1));

end

%Build (x,y) ordered pairs for y-int and computation of approximations

for i=1:2*n

orderpairs(i,1) = xnodevalues(int8(i/2),1);

if mod(i,2) == 1

orderpairs(i,2) = 1;

else

orderpairs(i,2) = 0;

end

end

%Build Basis Vectors Slope and Y-int

count = 1;

for k=1:n %Creates basis vectors first column is slope, second is y-int

for i=1:n-1 %Basis vector row number

if i == k

basis(i,1,k) = -(n-1)/(b-a);

basis(i,2,k) = -(basis(i,1,k)*orderpairs(count,1))+orderpairs(count,2);

count=count+1;

elseif i == (k-1)

basis(i,1,k) = (n-1)/(b-a);

basis(i,2,k) = -(basis(i,1,k)*orderpairs(count,1))+orderpairs(count,2);

count=count+1;

else

basis(i,1,k) = 0;

basis(i,2,k) = 0;

end

end

end

%Builds the Local Stiffness Matrix

for i=1:n-1 %Basis_1 vector number

for j=1:2 %Basis_2 vector number

for k=1:2

total = 0;

fun = basis(i,1,j + (i-1))*basis(i,1,k + (i-1));

if fun ~= 0

total = total + fun*(xnodevalues(i+1,1)- xnodevalues(i,1));

end

localstiffness(j,k)=total;

end

end

globalstiffness(i,i) = globalstiffness(i,i) + localstiffness(1,1);

globalstiffness(i+1,i) = globalstiffness(i+1,i) + localstiffness(2,1);

globalstiffness(i,i+1) = globalstiffness(i,i+1) + localstiffness(1,2);

globalstiffness(i+1, i+1) = globalstiffness(i+1,i+1) + localstiffness(2,2);

end

% Builds the local Load Vector

for i=1:n %Basis vector number

total = 0;

for j=1:n-1

total = total + (xnodevalues(j+1,1)-xnodevalues(j,1))*((neg_u_double_prime(xnodevalues(j+1,1))*(basis(j,1,i)*...

xnodevalues(j+1,1)+basis(j,2,i))+(neg_u_double_prime(xnodevalues(j,1))*(basis(j,1,i)*xnodevalues(j,1)+basis(j,2,i)))))/2;

%Trapezoidal Rule (b-a)(f(a)+f(b))/2 finds (f, Phi_i)

end

total = total - alpha*globalstiffness(i,1) - beta*globalstiffness(i,n);

globalload(i,1) = globalload(i,1) + total;

end

adjustedglobalstiffness = zeros(n-2, n-2);

adjustedload = zeros(n-2,1);

for i = 1:n-2

for j = 1:n-2

adjustedglobalstiffness(i,j) = globalstiffness (i+1, j+1);

end

adjustedload(i,1) = globalload (i+1,1);

end

coeff = adjustedglobalstiffness adjustedload;

finalcoeff = zeros(n,1);

for i = 1:n

if i==1

finalcoeff(i) = alpha;

elseif i == n

finalcoeff(i) = beta;

else

finalcoeff(i) = coeff(i-1);

end

end

coeff(1) = alpha;

coeff(n) = beta; %Setting boundary values

plot(xnodevalues(1:n),finalcoeff(1:n),'b-x'); %Graphs finite solution

hold on; %Call to graph on the same plot

fplot(@original_function,(a:b),'color','r'); %Graphs the original function

Explanation / Answer

total = total - alpha*globalstiffness(i,1) - beta*globalstiffness(i,n);

globalload(i,1) = globalload(i,1) + total;

end

adjustedglobalstiffness = zeros(n-2, n-2);

adjustedload = zeros(n-2,1);

for i = 1:n-2

for j = 1:n-2

adjustedglobalstiffness(i,j) = globalstiffness (i+1, j+1);

end

adjustedload(i,1) = globalload (i+1,1);

end

coeff = adjustedglobalstiffness adjustedload;

finalcoeff = zeros(n,1);

for i = 1:n

if i==1

finalcoeff(i) = alpha;

elseif i == n

finalcoeff(i) = beta;

else

finalcoeff(i) = coeff(i-1);

end

end

coeff(1) = alpha;

coeff(n) = beta; %Setting boundary values

plot(xnodevalues(1:n),finalcoeff(1:n),'b-x'); %Graphs finite solution

hold on; %Call to graph on the same plot

fplot(@original_function,(a:b),'color','r'); %Graphs the original function

Related Questions

Navigate

• Browse (All)
• Browse I
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.