Error: ode45 Must return a column vector?
1 view (last 30 days)
Show older comments
King Hin Wong
on 16 May 2021
Commented: King Hin Wong
on 16 May 2021
I'm trying to solve this problem involving a 2DOF system involving springs and dampers, and i've coded all of this below, including the state space matrix, this is urgent as this is due in 2 days, any help would be appreciated :)
Many thanks in advance
function batch_timeresp
% the total time domain response
% System parameters
k3 = 15; k1 = 10.0; k2 = 10.0; m1 = 2.0; m2 = 6.0; L = 1.0; c = 1; ct = 0.5; g = 9.81;kequi=5;
%ball
zeta1 = (c/(2*sqrt(m2*k3)));
%block
zeta2 = (ct/(2*sqrt(m1*kequi)));
%damped natural freq of block
om02=sqrt(kequi/m2);
omd2=om02*(sqrt(1-(zeta2)^2));
%damped natural freq of ball
om01=sqrt(k3/m1);
omd1=om01*(sqrt(1-(zeta1)^2));
T1 = 0.01*((2*pi)/om02);
T2 = 1.50*((2*pi)/om01);
Tpulse = T2-T1;
% Excitation parameters
Ft = 30;
Im = ((Ft*Tpulse)/2);
% State space system matrix and coupling
A=[0,0,1,0;0,0,0,1;-(k1+k2)/(m1+m2),0,-c/(m1+m2),0;0,-k3/(m1*(L^2)),0,-ct/(m1*(L^2))]; %state space matrix
S=[1,0,0,0;0,1,0,0;0,0,1,(m1*L)/(m1+m2);0,0,(m1*L)/(m1*(L^2)),1]; %coupling
Sinv=inv(S);
gmat = [0;0;-g;g/L]; %gravity terms
Fmat=[0;0;Im/(m1+m2);(2*Im*L)/(m1*(L^2))]; %impulse force
% Numerical integration
tspan=linspace(0,5,1e3);
X0=[0;0;0;0];
[T,X]=ode45(@ssmodel,tspan,X0);
% Results
figure
subplot(2,1,1), plot(T,X(:,1)), ylabel('x [m]')
subplot(2,1,2), plot(T,X(:,2)), ylabel('dx/dt [m/s]'), xlabel('time [s]')
function dx=ssmodel(t,x)
% State-space model of 1 DOF system
dx=(A*x+gmat+Fmat).*Sinv;
end
end
0 Comments
Accepted Answer
Alan Stevens
on 16 May 2021
Like this?
% the total time domain response
% System parameters
k3 = 15; k1 = 10.0; k2 = 10.0; m1 = 2.0; m2 = 6.0; L = 1.0; c = 1; ct = 0.5; g = 9.81;kequi=5;
%ball
zeta1 = (c/(2*sqrt(m2*k3)));
%block
zeta2 = (ct/(2*sqrt(m1*kequi)));
%damped natural freq of block
om02=sqrt(kequi/m2);
omd2=om02*(sqrt(1-(zeta2)^2));
%damped natural freq of ball
om01=sqrt(k3/m1);
omd1=om01*(sqrt(1-(zeta1)^2));
T1 = 0.01*((2*pi)/om02);
T2 = 1.50*((2*pi)/om01);
Tpulse = T2-T1;
% Excitation parameters
Ft = 30;
Im = ((Ft*Tpulse)/2);
% State space system matrix and coupling
A=[0,0,1,0;0,0,0,1;-(k1+k2)/(m1+m2),0,-c/(m1+m2),0;0,-k3/(m1*(L^2)),0,-ct/(m1*(L^2))]; %state space matrix
S=[1,0,0,0;0,1,0,0;0,0,1,(m1*L)/(m1+m2);0,0,(m1*L)/(m1*(L^2)),1]; %coupling
Sinv=inv(S);
gmat = [0;0;-g;g/L]; %gravity terms
Fmat=[0;0;Im/(m1+m2);(2*Im*L)/(m1*(L^2))]; %impulse force
% Numerical integration
tspan=linspace(0,5,1e3);
X0=[0;0;0;0];
[T,X]=ode45(@(t,x) ssmodel(t,x,A,gmat,Fmat,Sinv),tspan,X0);
% Results
figure
subplot(2,1,1), plot(T,X(:,1)), ylabel('x [m]')
subplot(2,1,2), plot(T,X(:,2)), ylabel('dx/dt [m/s]'), xlabel('time [s]')
function dx=ssmodel(~,x,A,gmat,Fmat,Sinv)
% State-space model of 1 DOF system
dx=Sinv*(A*x+gmat+Fmat); %%%%%%%%%%%%%%%%%%%%%%
end
5 Comments
Alan Stevens
on 16 May 2021
You are solving for four parameters, which means the function must return a 4x1 vector. If you have
dx=(A*x+gmat+Fmat).*Sinv;
this returns a 4x4 matrix.
Your state-space equation is presumably
S*dx = A*x+gmat+Fmat;
which means you need to multiply from the left to get x' on its own.
Sinv*S*dx = Sinv*(A*x+gmat+Fmat);
or
dx = Sinv*(A*x+gmat+Fmat);
More Answers (0)
See Also
Categories
Find more on Chassis Systems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!