must return a column vector
1 view (last 30 days)
Show older comments
I know this is a dimension issue, but i just can't see it at the moment. could someone help me out.
n=12;
x=zeros(1,2*n);
ct=zeros(n,n);
l=3; % in m
t=6;
%%%%%%%%%mass matrix %%%%%%%%%%%%%%
g=5000;
m=[ g,0,0,0,0,0,0,0,0,0,0,0;
0,g,0,0,0,0,0,0,0,0,0,0;
0,0,g,0,0,0,0,0,0,0,0,0;
0,0,0,g,0,0,0,0,0,0,0,0;
0,0,0,0,g,0,0,0,0,0,0,0;
0,0,0,0,0,g,0,0,0,0,0,0;
0,0,0,0,0,0,g,0,0,0,0,0;
0,0,0,0,0,0,0,g,0,0,0,0;
0,0,0,0,0,0,0,0,g,0,0,0;
0,0,0,0,0,0,0,0,0,g,0,0;
0,0,0,0,0,0,0,0,0,0,g,0;
0,0,0,0,0,0,0,0,0,0,0,g];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%Damping ratio matrix%%%%%%%
h=0.1;
zeta=[ h,0,0,0,0,0,0,0,0,0,0,0;
0,h,0,0,0,0,0,0,0,0,0,0;
0,0,h,0,0,0,0,0,0,0,0,0;
0,0,0,h,0,0,0,0,0,0,0,0;
0,0,0,0,h,0,0,0,0,0,0,0;
0,0,0,0,0,h,0,0,0,0,0,0;
0,0,0,0,0,0,h,0,0,0,0,0;
0,0,0,0,0,0,0,h,0,0,0,0;
0,0,0,0,0,0,0,0,h,0,0,0;
0,0,0,0,0,0,0,0,0,h,0,0;
0,0,0,0,0,0,0,0,0,0,h,0;
0,0,0,0,0,0,0,0,0,0,0,h];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%Stiffness matrix%%%%%%%%%%
b=4*((3*100)/l); %kN/m
k=[ 2*b,-b,0,0,0,0,0,0,0,0,0,0;
-b,2*b,-b,0,0,0,0,0,0,0,0,0;
0,-b,2*b,-b,0,0,0,0,0,0,0,0;
0,0,-b,2*b,-b,0,0,0,0,0,0,0;
0,0,0,-b,2*b,-b,0,0,0,0,0,0;
0,0,0,0,-b,2*b,-b,0,0,0,0,0;
0,0,0,0,0,-b,2*b,-b,0,0,0,0;
0,0,0,0,0,0,-b,2*b,-b,0,0,0;
0,0,0,0,0,0,0,-b,2*b,-b,0,0;
0,0,0,0,0,0,0,0,-b,2*b,-b,0;
0,0,0,0,0,0,0,0,0,-b,2*b,-b;
0,0,0,0,0,0,0,0,0,0,-b,2*b];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=m^(-0.5);
kt=a*k*a;
[p,q]=eig(kt);
omega=sqrt(q);
for zt=1:n
ct(zt,zt)=zeta(zt)*2*omega(zt,zt);
end
af=a^(-1);
c=af*ct*af;
v0=[0;0;0;0;0;0;0;0;0;0;0;0];%initial velocities
x0=[0;0;0;0;0;0;0;0;0;0;0;0];%initial displacement
ti=0; %initial time
tf=100; %final time
init=[v0,x0]; %initial velocity and displacement initialize
[t,ny]=ode45(@yd,[ti,tf],init);
for kl=1:n
plot(t,ny(:,kl+n),'b');
title(['Response x',num2str(kl)]);
xlabel('Time, seconds');
ylabel('Displacement, m');
grid
figure(gcf)
pause
end
here is the file it references
if true
function yd=yd(t,x)
yd(1)=(-((800)*x(13)+(68.1858)*x(1)*(-400)*x(14))+(0))/(5000);
yd(2)=(-((-400)*x(13)*(800)*x(14)*(-400)*x(15))+(0))/(5000);
yd(3)=(-((-400)*x(14)*(800)*x(15)*(-400)*x(16))+(0))/(5000);
yd(4)=(-((-400)*x(15)*(800)*x(16)*(-400)*x(17))+(0))/(5000);
yd(5)=(-((-400)*x(16)*(800)*x(17)*(-400)*x(18))+(0))/(5000);
yd(6)=(-((-400)*x(17)*(800)*x(18)*(-400)*x(19))+(0))/(5000);
yd(7)=(-((-400)*x(18)*(800)*x(19)*(-400)*x(20))+(0))/(5000);
yd(8)=(-((-400)*x(19)*(800)*x(20)*(-400)*x(21))+(0))/(5000);
yd(9)=(-((-400)*x(20)*(800)*x(21)*(-400)*x(22))+(0))/(5000);
yd(10)=(-((-400)*x(21)*(800)*x(22)*(-400)*x(23))+(0))/(5000);
yd(11)=(-((-400)*x(22)*(800)*x(23)*(-400)*x(24))+(0))/(5000);
yd(12)=(-((-400)*x(23)*(800)*x(24))+(5000*sin(3*t)))/(5000);
yd(13)=x(1);
yd(14)=x(2);
yd(15)=x(3);
yd(16)=x(4);
yd(17)=x(5);
yd(18)=x(6);
yd(19)=x(7);
yd(20)=x(8);
yd(21)=x(9);
yd(22)=x(10);
yd(23)=x(11);
yd(24)=x(12);
end
0 Comments
Accepted Answer
Star Strider
on 11 Dec 2014
I didn’t run your code, but see if making the first statement after the function statement a preallocation of ‘yd’ to be a column vector will help:
yd = zeros(24,1);
Also, it’s not good practice to name your function the same as the variables within it or the variable it returns. Change the function name to ‘ydf’ or something that is meaningful in the context of your code. Also change the name of the .m file containing it and of course the reference in ode45.
0 Comments
More Answers (1)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!