Asked by Shubham Gorde
on 8 May 2018 at 21:05

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)

Answer by Walter Roberson
on 8 May 2018 at 21:19

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)

Shubham Gorde
on 18 May 2018 at 0:27

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

Walter Roberson
on 18 May 2018 at 0:33

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

Shubham Gorde
about 18 hours ago

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.