- You will see updates in your activity feed.
- You may receive emails, depending on your notification preferences.

- 3 questions asked
- 0 answers
- 0 accepted answers
- reputation: 0

- 3 questions asked
- 0 answers
- 0 accepted answers
- reputation: 0

Accepted Answer by Walter Roberson

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

4 views (last 30 days)

4 views (last 30 days)

syms pi p q E A l x

n=3;

k=(n:n);

r=(1:n);

a=(1:n);

N=sym(1:n);

pi=eval(pi);

p=q*x/E*A*l;

for i=1:n

syms x

N(i)=(x.^2-2*x).^i;

%N(i)=sin((2*i-1)/2*pi*x);

%N(i)=piecewise(x>=(i-1)/n & x<1/n+(i-1)/n,n*(x-(i-1)/n),x>=i/n & x<=(i+1)/n & x<1,-n*x+i+1,0);

Y1=(N(i));

Y2=diff(N(i),x);

Y3=diff(N(i),x,2);

x=linspace(0,1);

plot(x,subs(Y1))

hold on

plot(x,subs(Y2))

hold on

plot(x,subs(Y3))

hold off

legend('N','N''','N''''');

pause(1)

end

for i=1:n

r(i)=-int(N(i)*p,x=linspace(0,1))*l ; %Galerkin

%r(i)=-int(diff(N(i),x,2)*p,x=linspace(0,1))/l ; %Least squares

for j=1:n

%k(i,j)=1/l*int(N(i)*diff(N(i),x,2),x=linspace(0,1)); %Galerkin

%k(i,j)=1/l^3*int(diff (N(j),x,2)*diff(N(i),x,2),x=linspace(0,1)); %Least Squares

k(i,j)=-1/l*int(diff(N(i),x)*diff(N(j),x),x=linspace(0,1)); %Displacement formulation

end

end

display(k)

display(r)

a=linsolve(k,r);

U=q*l^2*(x/2-x^3/6)/E*A;

u=o;

for i=1:n

u=u+a(i)*N(i);

end

u=simplify(u);

plot((U/(q*l^2/E*A)),(u/(q*l^2/E*A),x=linspace(0,1))

X1=0:0.9999;

Z3=(U/(q*l^2/E/A)-u/(q*l^2/E/A));

plot(X1,Z3)

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

Answer by Walter Roberson
on 8 May 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

Accepted Answer

This is as close as you can get:

syms pi p q E A l x

n=3;

k = sym(n:n);

r = sym(1:n);

a=(1:n);

N=sym(1:n);

pi=eval(pi);

p=q*x/E*A*l;

for i=1:n

syms x

N(i)=(x.^2-2*x).^i;

%N(i)=sin((2*i-1)/2*pi*x);

%N(i)=piecewise(x>=(i-1)/n & x<1/n+(i-1)/n,n*(x-(i-1)/n),x>=i/n & x<=(i+1)/n & x<1,-n*x+i+1,0);

Y1=(N(i));

Y2=diff(N(i),x);

Y3=diff(N(i),x,2);

x=linspace(0,1);

plot(x,subs(Y1))

hold on

plot(x,subs(Y2))

hold on

plot(x,subs(Y3))

hold off

legend('N','N''','N''''');

pause(1)

end

syms x for i=1:n

r(i)=-int(N(i)*p,x,0, 1)*l ; %Galerkin

%r(i)=-int(diff(N(i),x,2)*p,x=linspace(0,1))/l ; %Least squares

for j=1:n

%k(i,j)=1/l*int(N(i)*diff(N(i),x,2),x=linspace(0,1)); %Galerkin

%k(i,j)=1/l^3*int(diff (N(j),x,2)*diff(N(i),x,2),x=linspace(0,1)); %Least Squares

k(i,j)=-1/l*int(diff(N(i),x)*diff(N(j),x),x, 0, 1); %Displacement formulation

end

end

display(k)

display(r)

a=linsolve(k,r.');

U=q*l^2*(x/2-x^3/6)/E*A;

u = sym(0);

for i=1:n

u=u+a(i)*N(i);

end

u=simplify(u);

X = linspace(0,1); expr1 = subs((U/(q*l^2/E*A)), x, X); expr2 = subs((u/(q*l^2/E*A)), x, X); plot(expr1, expr2)

X1=0:0.9999;

Z3 = subs((U/(q*l^2/E/A)-u/(q*l^2/E/A)), x, X1);

plot(X1,Z3)

The two plot() at the end will fail. The first one will fail because expr2 involves the unresolved symbolic variable `l` (lower-case L). The second one will fail because Z3 involves the unresolved symbolic variables `A` and `l` (lower-case L)

It shows the error at

r(i)=-int(N(i)*p,x,0, 1)*l ; %Galerkin

as its maple version is

r(i):=-int(N[i]*p,x=0..1)*l:

and for this part:

X = linspace(0,1);

expr1 = subs((U/(q*l^2/E*A)), x, X);

expr2 = subs((u/(q*l^2/E*A)), x, X);

plot(expr1, expr2)

X1=0:0.9999;

Z3 = subs((U/(q*l^2/E/A)-u/(q*l^2/E/A)), x, X1);

plot(X1,Z3)

Its maple version is:

plot([U/(q*l^2/E/A),(u/(q*l^2/E/A)],x=0..1);

plot([U/(q*l^2/E/A)-(u/(q*l^2/E/A)],x=0..0.9999999);

Walter Roberson
on 8 May 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

What error is it showing? Did you miss my change

r = sym(1:n);

You still have the problems that A and l (lower-case L) are undefined. I noticed that you assigned to `a` (lower-case A) but that is a different variable than A is.

However, you should change X1 to

X1 = linspace(0,0.9999);

By the way, you would not use

r(i):=-int(N[i]*p,x=0..1)*l:

for Maple. The Maple version would be

r[i]:=-int(N[i]*p,x=0..1)*l:

It is giving value for force vector r =

[ (5*A*l^2*q)/(12*E), 2, 3]

The following error occurred converting from sym to double: Error using mupadmex Error in MuPAD command: DOUBLE cannot convert the input expression into a double array. If the input expression contains a symbolic variable, use the VPA function instead.

Error in Galerkin (line 34) k(i,j)=-1/l*int(diff(N(i),x)*diff(N(j),x),x,0,1); %Displacement formulation

Actual values for force vectors are:

r1=(5*q*l)/(12*E*A)

r2=-(11*q*l)/(30*E*A)

r3=(93*q*l)/(280*E*A)

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

Yeah, I have missed some of the changes you made.Thank you for the help.

U=q*l^2*(x/2-x^3/6)/E*A;

If divide the expression U/(q*l^2/E*A) it wont have A,E l(lower case L) and same with u/(q*l^2/E*A),then this expression will contain only x which varies from 0 to 1.Plot is not getting done as the error is: Error using plot Conversion to double from sym is not possible.

Walter Roberson
on 9 May 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

Are you doing subs() to bring in the numeric x values? I need to see your revised code lines.

In revised code in maple, after a:=linsolve(k,r);

U=q*l^2*(x/2-x^3/6)/E*A; #analytical solution for linear load

u:=0;

for i from 1 by 1 to n do

u:=u+a[i]*N[i]:

end do:

u:=simplify(u);

plot([U/(q*l^2/E*A),u/(q*l^2/E*A)],x=0..1);

plot(U/(q*l^2/E*A)-u/(q*l^2/E*A),x=0..0.99999);

so for this U that is displacement is a function of x (varies from 0 to 1) and will get only numerical values which helped me for plotting the graph in maple.

for last plot it plots the difference between analytical solution with approximation method to help to realize the accuracy.

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

U involves l^2 so when you divide by l^2 the l cancel out, leaving the first expression in the plot with only variables that have been assigned values, so the U portion of the expression can be plotted.

u involves l^3, so when you divide by l^2 there is still an unresolved l left over.

You can attach your maple code as text; I have maple and can go through and compare with the MATLAB code.

Right now I have access to the pdf so have attached for your reference with graphs.

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

syms p q E A l x

n=3;

k = zeros(n,n,'sym');

r = zeros(1,n,'sym');

N = zeros(1,n,'sym');

Pi = sym(pi);

p = q/E/A*x; %linear load

X = linspace(0,1) .';

for i=1:n N(i) = (x.^2-2*x).^i; Y1 = N(i); Y2 = diff(N(i),x); Y3 = diff(N(i),x,2);

plot(X,subs([Y1, Y2, Y3], x, X)); ylim([-6 6]); title(char(N(i))); legend('N','N''','N''''');

pause(1)

end

for i=1:n

r(i) = -int(N(i)*p, x, 0, 1)*l ; %Galerkin

for j=1:n

k(i,j) = -1/l*int(diff(N(i),x)*diff(N(j),x), x, 0, 1); %Displacement formulation

end

end

disp(k)

disp(r)

a = linsolve(k,r.');

U = q*l^2/E/A*(x/2-x^3/6);

u = sym(0);

for i=1:n

u = u + a(i)*N(i);

end

u = simplify(u);

expr1 = U/(q*l^2/E/A); expr2 = u/(q*l^2/E/A); plot(X, subs([expr1, expr2], x, X));

X1 = linspace(0, 0.9999);

Z3 = expr1 - expr2;

plot(X1, subs(Z3,x,X));

The above code you did is not taking piecewise shape function.

In the maple file I have attached above, the piecewise shape function is

N(i)=(x>=(i-1)/n & x<1/n+(i-1)/n,n*(x-(i-1)/n),x>=i/n & x<(i+1)/n & x<1,-n*x+i+1,0)

Also what is the keyword for piecewise in matlab?

Walter Roberson
on 16 May 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

You are mistaken. In Maple, # is the comment character. The line that defined N in piecewise form is commented out in the Maple code.

In R2016b and later, piecewise definition is done by calling piecewise().

Yes that was commented out but we want to check for every shape functions two of them are working and even I tried piecewise definition N(i)=piecewise(x>=(i-1)/n & x<1/n+(i-1)/n , n*(x-(i-1)/n), x>=i/n & x<(i+1)/n & x<1, -n*x+i+1, 0); but error is

Undefined function 'piecewise' for input arguments of type 'sym'.

I am in learning mode,so your help is much appreciated.

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

Which MATLAB release are you using?

I am using R2014a

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

I recommend that you upgrade MATLAB releases.

Yeah I will upgrade.I have one more question but this is about gauss integration over 4 node isoparametric element I have tried to extract this code from maple code into matlab but problem is for n=2 the answer I am getting is 52.5201 but actual answer is 83.27.I am not able to detect what is the problem.Will you please look into my code and the maple code attached?

This is my code:

syms xi eta x1 x2 x3 x4 y1 y2 y3 y4

N1=1/4*(1+xi )*(1+eta);

N2=1/4*(1-xi)*(1+eta);

N3=1/4*(1-xi)*(1-eta);

N4=1/4*(1+xi )*(1-eta);

x= N1 * x1 + N2*x2 + N3*x3 + N4*x4;

y= N1*y1 + N2*y2 + N3*y3 + N4*y4;

dN=[diff(N1,xi),diff(N2,xi),diff(N3,xi),diff(N4,xi);diff(N1,eta),diff(N2,eta),diff(N3,eta),diff(N4,eta)];

XY=[x1,y1;x2,y2;x3,y3;x4,y4];

Jac=dN*XY;

x1=2;y1=3;

x2=-1;y2=4;

x3=-4;y3=-2;

x4=3;y4=-1;

disp(XY)

xig(1,1)=0;Hi(1,1)=2;

xig(2,1)=(1/sqrt(3));Hi(2,1)=1;xig(2,2)=-(1/sqrt(3));Hi(2,2)=1;

xig(3,1)=sqrt(0.6);Hi(3,1)=5/9;xig(3,2)=0;Hi(3,2)=8/9;xig(3,3)=-sqrt(0.6);Hi(3,3)=5/9;

xig(4,1)=0.861136311594953;Hi(4,1)=0.347854845137454;xig(4,2)=0.339981043584856;Hi(4,2)=0.652145154862546;

xig(4,3)=-0.861136311594953;Hi(4,3)=0.347854845137454;xig(4,4)=-0.339981043584856;Hi(4,4)=0.652145154862546;

xig(5,1)=0.906179845938664;Hi(5,1)=0.236926885056189;xig(5,2)=0.538469310105683;Hi(5,2)=0.478628670499366;

xig(5,3)=0;Hi(5,3)=0.56888888888888889;xig(5,4)=-0.906179845938664;Hi(5,4)=0.236926885056189;

xig(5,5)=-0.538469310105683;Hi(5,5)=0.478628670499366;

%simplify((x*y)^2*det(Jac));

n=2;

inte=0;

for i=1:n

for j=1:n

xi=xig(n,i);

eta=xig(n,j);

inte=inte+((x*y)^2*det(Jac)* Hi(n,i) * Hi(n,j));

end

end

eval(inte)

Walter Roberson
on 18 May 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

clearvars xi = sym('xi'); eta = sym('eta'); x1 = sym('x1'); x2 = sym('x2'); x3 = sym('x3'); x4 = sym('x4'); y1 = sym('y1'); y2 = sym('y2'); y3 = sym('y3'); y4 = sym('y4'); Q = @(v) sym(v, 'r');

N1=1/4*(1+xi )*(1+eta);

N2=1/4*(1-xi)*(1+eta);

N3=1/4*(1-xi)*(1-eta);

N4=1/4*(1+xi )*(1-eta);

x= N1 * x1 + N2*x2 + N3*x3 + N4*x4;

y= N1*y1 + N2*y2 + N3*y3 + N4*y4;

dN=[diff(N1,xi),diff(N2,xi),diff(N3,xi),diff(N4,xi); ... diff(N1,eta),diff(N2,eta),diff(N3,eta),diff(N4,eta)];

XY=[x1,y1; ... x2,y2; ... x3,y3; ... x4,y4];

Jac=dN*XY;

x1=Q(2);y1=Q(3);

x2=Q(-1);y2=Q(4);

x3=Q(-4);y3=Q(-2);

x4=Q(3);y4=Q(-1);

disp(subs(XY))

xig(1,1)=Q(0);Hi(1,1)=Q(2);

xig(2,1)=(1/sqrt(Q(3)));Hi(2,1)=1;xig(2,2)=-(1/sqrt(Q(3)));Hi(2,2)=1;

xig(3,1)=sqrt(Q(0.6));Hi(3,1)=Q(5)/9;xig(3,2)=0;Hi(3,2)=Q(8)/9;xig(3,3)=-sqrt(Q(0.6));Hi(3,3)=Q(5)/9;

xig(4,1)=Q(0.861136311594953);Hi(4,1)=Q(0.347854845137454);xig(4,2)=Q(0.339981043584856);Hi(4,2)=Q(0.652145154862546);

xig(4,3)=Q(-0.861136311594953);Hi(4,3)=Q(0.347854845137454);xig(4,4)=Q(-0.339981043584856);Hi(4,4)=Q(0.652145154862546);

xig(5,1)=Q(0.906179845938664);Hi(5,1)=Q(0.236926885056189);xig(5,2)=Q(0.538469310105683);Hi(5,2)=Q(0.478628670499366);

xig(5,3)=0;Hi(5,3)=Q(0.56888888888888889);xig(5,4)=Q(-0.906179845938664);Hi(5,4)=Q(0.236926885056189);

xig(5,5)=Q(-0.538469310105683);Hi(5,5)=Q(0.478628670499366);

%simplify((x*y)^2*det(Jac));

n=2;

inte=0; for i=1:n

for j=1:n

xi=xig(n,i);

eta=xig(n,j);

inte = inte + subs(((x*y)^2*det(Jac)* Hi(n,i) * Hi(n,j)));

end

end

vpa(inte)

Yeah thank you so much for the input.Now I am learning what you did in the code.

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

The biggest change was to subs() the expression being added.

I have attached the image where I want to plot the isoparametric mapped element in cartesian coordinate which has curved part.How to plot that in matlab I tried this way:

xi=[-1,0,1,1,1,0,-1,-1,-1];

eta=[-1,-1,-1,0,1,1,1,0,-1];

plot(xi,eta,'-o');

axis([-2,2,-2,2]);

x1=2;y1=3;

x2=-1;y2=3;

x3=-1;y3=-1;

x4=2;y4=-1;

x5=0.5;y5=4.5;

x6=-1;y6=1;

x7=0.5;y7=-3;

x8=2;y8=1;

disp(subs(XY))

%X=[x1 x5 x2 x6 x3 x7 x4 x8 x1];

%Y=[y1 y5 y2 y6 y3 y7 y4 y8 y1];

rectangle('Position',[-2 2 3 3],'Curvature',0.8)

%plot(X,Y,'-o');

axis([-8 8 -8 8]);

Walter Roberson
on 2 Jun 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

I have plotted the curve using quadratic equation.Now I am working on plane stress problem.In this task the procedure is :

1.Calculate the global stiffness matrix and I am stuck at reduction matrix

2.Then using plane stress calculate the stiffness in symbolic form first and then in numerical values

3.Next step is for element 1 as in pdf the coordinates are for i=1,j=2 and k=3 there are x1,x2,x3 and y1,y2,y3 values, calculate area for element 1 = 0.5*det([1 1 1;x1 x2 x3;y1 y2 y3])

I tried this way:

syms x1 x2 x3 y1 y2 y3 z1 z2 z3

x1=0;y1=0;z1=1;

x2=3;y2=1;z2=1;

x3=2;y3=2;z3=1;

A =[z1 z2 z3; x1 x2 x3 ; y1 y2 y3];

Area=det(A);

I want to use syms form because if I chane the values of x1,x2,x3 then it can calculate accordingly

4. After that calculate b1= x2*y3-x3*y2, b2= x3*y1-x1*y3 ,b3= x1*y2-x2*y1 c1= y2-y3, c2=y3-y1,c3=y1-y2,d1=x3-x2,d2=x1-x3,d3=x2-x1 for element 1(in terms of i,j,k the formulae for b(i)=x(j)*y(k)-x(k)*y(j) c(i)=y(j)-y(k) and d(i)=x(k)-x(j) for element 1: i=1,j=2,k=3)

similarly for element 2 we need to calculate b(i),c(i) and d(i) but this time i=1,j=3,d=4 that means b1= x3*y4-x4*y3, b3= x4*y1-x1*y4 ,b4= x1*y3-x3*y1 , c1= y3-y4, c3=y4-y1, c4=y1-y3, d1=x4-x3, d3=x1-x4, d4=x3-x1 for element 2 then we will get stiffness matrices k22(element 1),k23(elt 1),k32(elt 1),k33(elt 1),k33(elt 2) according to plane stress formulae in pdf in terms of c and d for elements 1 and 2

4.Last step is calculate force vector r given in the pdf [100 0 100 0] and calculate displacement using K*a=r a=displacement vector

I really want your help in coding this problem?Attached the pdf

Walter Roberson
on 9 Jun 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

I seem to have missed your question in that?

Note: if you want to use syms for your Area equation, then

syms x1 x2 x3 y1 y2 y3 z1 z2 z3 A =[z1 z2 z3; x1 x2 x3 ; y1 y2 y3]; Area = det(A);

X1=0; Y1=0; Z1=1; X2=3; Y2=1; Z2=1; X3=2; Y3=2; Z3=1;

this_Area = subs(Area, [x1, x2, x3, y1, y2, y3, z1, z2, z3], [X1, X2, X3, Y1, Y2, Y3, Z1, Z2, Z3])

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

"I am stuck" is not an actionable question.

Actually with direct values I can code the problem but with for loop I want to convert full stiffness matrix (4*4)for element 1 and 2 to reduction stiffness matrix (2*2) and then with c and d values (also need to be in i,j form) will calculate reduction stiffness matrix in numerical form as on the page number 47 in the pdf?

I tried this way

syms x1 x2 x3 y1 y2 y3 z1 z2 z3

A =[z1 z2 z3; x1 x2 x3 ; y1 y2 y3];

E=100000; %Modulus of elasticity

h=1; %thickness

V=0.1; %poissons ratio

%for Element 1

Area1 = det(A);

X11=0; Y11=0; Z11=1;

X21=100; Y21=0; Z21=1;

X31=100; Y31=10; Z31=1;

this_Area1 = 0.5*subs(Area1, [x1, x2, x3, y1, y2, y3, z1, z2, z3],

[X11, X21, X31, Y11, Y21,Y31, Z11, Z21, Z31])

b11= X21*Y31-X31*Y21

b21= X31*Y11-X11*Y31

b31= X11*Y21-X21*Y11

c11= Y21-Y31

c21= Y31-Y11

c31= Y11-Y21

d11= X31-X21

d21= X11-X31

d31= X21-X11

%For Element 2

A =[z1 z2 z3; x1 x2 x3 ; y1 y2 y3]; Area2 = det(A);

X12=0; Y12=0; Z12=1;

X32=100; Y32=10; Z32=1;

X42=0; Y42=10; Z42=1;

this_Area2 = subs(Area2, [x1, x2, x3, y1, y2, y3, z1, z2, z3], [X12, X32, X42, Y12, Y32,Y42, Z12, Z32, Z42])

b12= X32*Y42-X42*Y32

b32= X42*Y12-X12*Y42

b42= X12*Y32-X32*Y12

c12= Y32-Y42

c32= Y42-Y12

c42= Y12-Y32

d12= X42-X32

d32= X12-X42

d42= X32-X12

%for element 1 and 2

syms c d V k r

k = zeros(n,n,'sym');

r = zeros(1,n,'sym');

n=2;

for i=1:n

j=1:n

kij=E*h/4*Area1*(1-V^2)*([c(i)*c(j)+d(i)*d(j)*(1- V)/2 ; c(i)*d(j)*V+c(j)*d(i)*(1-V)/2 ;

c(j)*d(i)*V+c(i)*d(j)*(1- V)/2 ; c(i)*c(j)+d(i)*d(j)*(1-V)/2])

end

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

To proceed, I need an error message, or something particular whose output is not what is expected (along with information about what is expected).

I have not done any FEM work in... 30 years? And it was not much then. You should assume that I do not know the theory of FEM. You should also assume that I am not going to go and study the theory of FEM for a few days so that I can return to this question to figure out how your uncommented code does not fit the theory. But I am good at tracking backwards things like why array sizes do not match in concatenation.

This is the code with comments:

syms x1 x2 x3 x4 y1 y2 y3 y4 x4 z1 z2 z3 z4

A =[z1 x1 y1; z2 x2 y2 ;z3 x3 y3];

E=100000; %Modulus of elasticity

h=1; %thickness

V=0.1; %poissons ratio

%for Element 1

Area1 = det(A);

X11=0; Y11=0; Z11=1; %these are the coordinates for element 1

X21=100; Y21=0; Z21=1;

X31=100; Y31=10; Z31=1;

this_Area1 = 0.5*subs(Area1, [x1, x2, x3, y1, y2, y3, z1, z2, z3],

[X11, X21, X31, Y11, Y21,Y31, Z11, Z21, Z31])

b11= X21*Y31-X31*Y21;

b21= X31*Y11-X11*Y31;

b31= X11*Y21-X21*Y11;

c11= Y21-Y31;

c21= Y31-Y11;

c31= Y11-Y21;

d11= X31-X21;

d21= X11-X31;

d31= X21-X11;

%For Element 2

A =[z1 z3 z4; x1 x3 x4 ; y1 y3 y4];

Area2 = det(A);

X12=0; Y12=0; Z12=1; % these are the coordinates for element 2

X32=100; Y32=10; Z32=1;

X42=0; Y42=10; Z42=1;

this_Area2 = subs(Area2, [x1, x3, x4, y1, y3, y4, z1, z3, z4], [X12, X32, X42, Y12, Y32,Y42, Z12, Z32, Z42]);

b12= X32*Y42-X42*Y32; % b12 means b1 value for element 2,likewise all b,c d values for element 2

b32= X42*Y12-X12*Y42;

b42= X12*Y32-X32*Y12;

c12= Y32-Y42;

c32= Y42-Y12;

c42= Y12-Y32;

d12= X42-X32;

d32= X12-X42;

d42= X32-X12;

%for element 1 and 2

k221=(E*h/4*Area1*(1-V^2))*([c21*c21+d21*d21*(1-V)/2 c21*d21*V+c21*d21*(1-V)/2; c21*d21*V+c21*d21*(1-V)/2 c21*c21+d21*d21*(1-V)/2]) % k22 2 by 2 stiffness matrix for element 1 and the output I am getting is

k221 =

[ 113850000*x1*y2*z3 - 113850000*x1*y3*z2 - 113850000*x2*y1*z3 + 113850000*x2*y3*z1 + 113850000*x3*y1*z2 - 113850000*x3*y2*z1, 13612500*x1*y3*z2 - 13612500*x1*y2*z3 + 13612500*x2*y1*z3 - 13612500*x2*y3*z1 - 13612500*x3*y1*z2 + 13612500*x3*y2*z1] [ 13612500*x1*y3*z2 - 13612500*x1*y2*z3 + 13612500*x2*y1*z3 - 13612500*x2*y3*z1 - 13612500*x3*y1*z2 + 13612500*x3*y2*z1, 113850000*x1*y2*z3 - 113850000*x1*y3*z2 - 113850000*x2*y1*z3 + 113850000*x2*y3*z1 + 113850000*x3*y1*z2 - 113850000*x3*y2*z1]

the expected output does not have x and y terms it has only numerical terms like

k221 = 500/99*[4600 -550 ; -550 10045]

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

Your K221 involves Area1, which is an area formula. You have also calculated this_Area1 but you do not use it.

When you calculate b11 and so on you should be using the symbolic variables. Then after you have calculated k211, subs() in the variables exactly the same way as for the areas.

I want to plot curve through 3 points so I found quadratic equation for 3 horizontal points and simply plotted it but when I tried to do with 3 vertical points it is not working.am I missing something here? Here is my code:

a=-1:0.01:2;

b=-0.66667*a.^2+0.66667*a+4.33333;

c=-1:0.01:2;

d=0.66667*c.^2-0.66667*c-2.3333;

e=-1:0.01:3;

f=-0.125*e.^2+0.25*e-0.625; %this is that curve equation passing % through 3 vertical points

X=[2 0.5 -1 -0.5 -1 0.5 2 1.5 2];

Y=[3 4.5 3 1 -1 -2.5 -1 1 3];

plot(X,Y,'-o',a,b,c,d,e,f);

axis([-8 8 -8 8]);

the three vertical points are:(-1,3),(-0.5,1),(-1,-1) and the equation I get for these points are typed in e and f part.

Walter Roberson
on 11 Jun 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

What formula did you try to use to find the equation passing through three points?

Walter Roberson
on 12 Jun 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

plot(X,Y,'-o',a,b,c,d,f,e)

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

To get the curve you created an equation for x in terms of y.

y=-1:0.01:3; x=-0.125*y.^2+0.25*y-0.625;

and you plot(x, y)

You happened to call y as e and you called x as f, so you plot(f, e)

You want to get two y values for each x value, but you also want to vary x linearly. In order to get two y values from a linearly changing x, your equation for y would need to produce two values for each x that was input. That would not look like a polynomial: it would look like a pair of roots of a quadratic equation.

Yeah.Got it.Thank you.

How to do b,c,d values in terms of i,j,k where i,j,k are 1,2,3 for b1 then anticlockwise rotation for b2 i,j,k changes to 2,3,1 then to for b3 i,j,k are 3,1,2 standard formulas are b(i) = x(j)*y(k) − x(k)*y(j) , c(i) = y(j) − y(k), d(i) = x(k) − x(j) where x(i),y(i),x(j),y(j),x(k),y(k) are coordinates of the triangle element

I want this b1,b2,b3,c1,c2,c3,d1,d2,d3 in terms of loop where for these coordinates it can calculate the values automatically.

X1=0; Y1=0; %these are the coordinates for element 1

X2=100; Y2=0;

X3=100; Y3=10;

b1= X2*Y3-X3*Y2;

b2= X3*Y1-X1*Y3;

b3= X1*Y2-X2*Y1;

c1= Y2-Y3;

c2= Y3-Y1;

c3= Y1-Y2;

d1= X3-X2;

d2= X1-X3;

d3= X2-X1;

I tried this way:

X(1)=0; Y(1)=0; Z(1)=1; %these are the coordinates for element 1

X(2)=100; Y(2)=0; Z(2)=1;

X(3)=100; Y(3)=10; Z(3)=1;

n=3;

%c(i)=zeros(1,n);

for i=1:n

c(i)=(Y(i+1)-Y(i+2)); b(i)= X(i+1)*Y(i+2)-X(i+2)*Y(i+1); d(i)= X(i+2)-X(i+1);

end

but not working?

Walter Roberson
on 16 Jun 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

I do not understand what you are doing, but perhaps the following will help:

b1ijk = 1 : 3; b2ijk = circshift(b1ijk, [1 -1]); b3ijk = circshift(b1ijk, [1 -2]);

for idx = 1 : length(b1ijk) [b1(b1ijk(idx)), b2(b2ijk(idx)), b3(b3ijk(idx))] end

Sometimes though it is easier to circshift the data vectors and then do straight indexing of that

b1s = b1; b2s = circshift(b2, [1 -1]); b3s = circshift(b3, [1 -2]); for idx = 1 : length(b1) [b1s(idx), b2s(idx), b3s(idx)] end

Actually I am solving plane stress triangular element problem.I solved this in this manner but want to standardized it for more number of elements.My code is this:

syms x1 x2 x3 x4 y1 y2 y3 y4 x4 z1 z2 z3 z4

A =[z1 x1 y1; z2 x2 y2 ;z3 x3 y3];

E=10000; %Modulus of elasticity

h=1; %thickness

V=0.1; %poissons ratio

%for Element 1

Area1 = det(A);

X11=0; Y11=0; Z11=1; %these are the coordinates for element 1

X21=100; Y21=0; Z21=1;

X31=100; Y31=10; Z31=1;

this_Area1 = 0.5*subs(Area1, [x1, x2, x3, y1, y2, y3, z1, z2, z3], [X11, X21, X31, Y11, Y21,Y31, Z11, Z21, Z31]);

b11= X21*Y31-X31*Y21;

b21= X31*Y11-X11*Y31;

b31= X11*Y21-X21*Y11;

c11= Y21-Y31;

c21= Y31-Y11;

c31= Y11-Y21;

d11= X31-X21;

d21= X11-X31;

d31= X21-X11;

%For Element 2

A =[z1 z3 z4; x1 x3 x4 ; y1 y3 y4];

Area2 = det(A);

X12=0; Y12=0; Z12=1; % these are the coordinates for element 2

X32=100; Y32=10; Z32=1;

X42=0; Y42=10; Z42=1;

this_Area2 =0.5* subs(Area2, [x1, x3, x4, y1, y3, y4, z1, z3, z4], [X12, X32, X42, Y12, Y32,Y42, Z12, Z32, Z42]);

b12= X32*Y42-X42*Y32; % b12 means b1 value for element 2,likewise all b,c d values for element 2

b32= X42*Y12-X12*Y42;

b42= X12*Y32-X32*Y12;

c12= Y32-Y42;

c32= Y42-Y12;

c42= Y12-Y32;

d12= X42-X32;

d32= X12-X42;

d42= X32-X12;

%for element 1 and 2

k221=E*h/4/this_Area1/(1-V.^2)*([c21*c21+d21*d21*(1-V)/2 c21*d21*V+c21*d21*(1-V)/2; c21*d21*V+c21*d21*(1-V)/2 d21*d21+c21*c21*(1-V)/2]); % matrix for k22 matrix for element 1

k231=E*h/4/this_Area1/(1-V.^2)*[c21*c31+d21*d31*(1-V)/2 c21*d31*V+c31*d21*(1-V)/2; c31*d21*V+c21*d31*(1-V)/2 d21*d31+c21*c31*(1-V)/2]; % matrix for k23 matrix for element 1

k321=E*h/4/this_Area1/(1-V.^2)*[c31*c21+d31*d21*(1-V)/2 c31*d21*V+c21*d31*(1-V)/2; c21*d31*V+c31*d21*(1-V)/2 d31*d21+c31*c21*(1-V)/2]; % matrix for k32 matrix for element 1

k331=E*h/4/this_Area1/(1-V.^2)*[c31*c31+d31*d31*(1-V)/2 c31*d31*V+c31*d31*(1-V)/2; c31*d31*V+c31*d31*(1-V)/2 d31*d31+c31*c31*(1-V)/2]; % matrix for k33 matrix for element 1

k332=E*h/4/this_Area2/(1-V.^2)*[c32*c32+d32*d32*(1-V)/2 c32*d32*V+c32*d32*(1-V)/2; c32*d32*V+c32*d32*(1-V)/2 d32*d32+c32*c32*(1-V)/2]; % matrix for k33 matrix for element 2

KG=[k221 k231; k321 k331+k332];

R=[100; 0 ;100; 0];

a=(KG)\R

I am getting the same answer as pdf but b,c,d,stiffness matrix are of manual input and this is limited to this particular problem. so wanted to make standard procedure for many problems.

The solution above u have given its not giving b1,b2,b3,c1,c2,c3,d1,d2,d3 values?

Walter Roberson
on 17 Jun 2018

- 30 questions asked
- 39,270 answers
- 14,804 accepted answers
- reputation: 76,494

I think I am a bit too tired to look at this today, sorry.

Will you able to check this today?

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply TodayUnable to complete the action because of changes made to the page. Reload the page to see its updated state.

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

## 0 Comments

Sign in to comment.