l=5; B=zeros(1,4); syms x; N1=1-(3*(x^2)/l^2)+(2*(x^3)/l^3); N2=x-(2*(x^2)/l)+((x^3)/l^2); N3=(3*(x^2)/l^2)-(2*(x^3)/l^3); N4=-((x^2)/l)+((x^3)/l^2); N=[N1 N2 N3 N4]; for i=1:4; a=diff(N(i,1),x); B(i,1)=a; end
MATLAB: The following error occurred converting from sym to double: DOUBLE cannot convert the input expression into a double array. Error in bending (line 12) B(i,1)=a;
differentiationdoublesymbolic
Related Solutions
%subdivisions space /time
ht=0.01; Tmax=1.2; nx=33; hx=1/(nx-1); x=[0:hx:1]';%matrices
K=stiff2(1/pi^2,hx,nx); M=mass(1/ht,hx,nx); A=M+K;%boundary conditions by deleting rows
A(nx,:)=[]; A(:,nx)=[]; A(1,:)=[]; A(:,1)=[];%creation
R=chol(A);%initial condition
u=cos(pi*(x-1/2));%time step loop
k=0; while (k*ht<Tmax)k=k+1;b=M*u;b(nx)=[];b(1)=[];%solve
u=zeros(nx,1);u(2:nx-1)=R\(R'\b);end% Set the matrices to zero
k=zeros(6,6,2); K=zeros(6,6,2); Gamma=zeros(6,6,2);% Enter parameter values, in units: mm^2, mm^4, and MPa(10^6 N/m)
A=6000; II=200*10^6; EE=200000;% Convert units into meters and kN
A=A/10^6; II = II/10^12; EE =EE*1000;% Element 1
i=[0,0]; j=[7.416,3];[k(:,:,1),K(:,:,1),Gamma(:,:,1)]=stiff(EE,II,A,i,j);% Element 2
i=j; j=[15.416,3];[k(:,:,2),K(:,:,2),Gamma(:,:,2)]=stiff(EE,II,A,i,j);% Define the elementary stiffness matrix
ID=[-4 1 -7;-5 2 -8;-6 3 -9];% Define the connection matrix
LM=[-4 -5 -6 1 2 3; 1 2 3 -7 -8 -9];% Assemble the augmented stiffness matrix
Kaug = zeros(9); for elem=1:2for r=1:6lr=abs(LM(elem,r));for c=1:6lc=abs(LM(elem,c));Kaug(lr,lc)=Kaug(lr,lc)+K(r,c,elem);endendend% Extract the stiffness matrix
Ktt=Kaug(1:3,1:3);% Determine the reactions at the nodes in local coordinates
fea(1:6,1)=0; fea(1:6,2)=[0,8*4/2,4*8^2/12,0,8*4/2,-4*8^2/12]';% Determine the reactions in the global coordinate system
FEA(1:6,1)=Gamma(:,:,1)'*fea(1:6,1);FEA(1:6,2)=Gamma(:,:,2)'*fea(1:6,2);% FEA_Rest for constrained nodes
FEA_Rest=[0,0,0,FEA(4:6,2)'];% Assemble the right-hand side for non-constrained nodes
P(1)=50*3/8; P(2)=-50*7.416/8-FEA(2,2); P(3)=-FEA(3,2);% Solve to find the displacements in meters and in radians
Displacements=inv(Ktt)*P'% Extract Kut
Kut = Kaug(4:9,1:3);% Compute the reactions and introduce boundary conditions
Reactions=Kut*Displacements+FEA_Rest'% Solve to find the internal forces, excluding fixed points
dis_global(:,:,1)=[0,0,0,Displacements(1:3)'];dis_global(:,:,2)=[Displacements(1:3)',0,0,0]; for elem=1:2dis_local= Gamma(:,:,elem)*dis_global(:,:,elem)';int_forces= k(:,:,elem)*dis_local+fea(1:6,elem)end%The above script calls the stiff function, which can be implemented as follows:
function [k,K,Gamma] = stiff( EE,II,A,i,j )% Find the length
L=sqrt((j(2)-i(2))^2+(j(1)-i(1))^2);% Compute the angle theta
if(j(1)-i(1))~=2alpha=atan((j(2)-i(2))/(j(1)-i(1)))elsealpha=-pi/2;end% Form the rotation matrix Gamma
Gamma =[cos(alpha) sin(alpha) 0 0 0 0;-sin(alpha) cos(alpha) 0 0 0 0;0 0 1 0 0 0;0 0 0 cos(alpha) sin(alpha) 0;0 0 0 -sin(alpha) cos(alpha) 0;0 0 0 0 0 1];% Form the elementary stiffness matrix in local coordinates
EI=EE*II; EA=EE*A; k=[EA/L, 0, 0, -EA/L, 0, 0; 0, 12*EI/L^3, 6*EI/L^2, 0, -12*EI/L^3,6*EI/L^2;0, 6*EI/L^2, 4*EI/L, 0 -6*EI/L^2, 2*EI/L;-EA/L, 0 ,0 , EA/L, 0, 0;0, -12*EI/L^3, -6*EI/L^2, 0, 12*EI/L^3, -6*EI/L^2;0, 6*EI/L^2, 2*EI/L, 0, -6*EI/L^2, 4*EI/L];% Elementary matrix in global coordinates
K=Gamma'*k*Gamma;endfunction R=stiff2(nu,h,n)%
[m1,m2]=size(nu);if (length(nu)==1)ee=nu*ones(n-1,1)/h; e1=[ee;0]; e2=[0;ee]; e=e1+e2;elseee=.5*(nu(1:n-1)+nu(2:n))/h; e1=[ee;0]; e2=[0;ee]; e=e1+e2;endR=spdiags([-e1 e -e2],-1:1,n,n);returnendfunction M=mass(alpha,h,n)[m1,m2]=size(alpha);if(length(alpha)==1)ee=alpha*h*ones(n-1,1)/6;e1=[ee;0]; e2=[0;ee]; e=e1+e2;elseee=h*(alpha(1:n-1)+alpha(2:n))/12;e1=[ee;0]; e2=[0;ee]; e=e1+e2;endM=spdiags([e1 2*e e2],-1:1,n,n);returnendTry symbolic expression
clc,clearalpha=1.5; A=0.0207; gamma=0.25;syms beta epsx = (beta-alpha*(beta-gamma)+sqrt((beta-alpha*(beta-gamma)).^2+4*alpha*beta*gamma*A))./(2*beta);y = x.*(1-x).*beta/gamma;t = -x+x.*y.*(alpha-eps*beta)./((A+x+y).^2);d = eps*beta.*x.*y.*(x.^2+(A+x).*x+A*alpha)./((A+x+y).^3); % convert symbolic expression to function handle
F = matlabFunction(t.^2 - d,'vars',[beta eps]);ezplot(F,[0,1,0,1]);
Best Answer