Error: ode45 Must return a column vector?

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

 Réponse acceptée

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 commentaires

King Hin Wong
King Hin Wong le 16 Mai 2021
Modifié(e) : King Hin Wong le 16 Mai 2021
Yes! Thank you so much!!! If you wouldn't mind, could you please explain what you've done in the function?
dx=ssmodel(~,x,A,gmat,Fmat,Sinv)
What does the tilde do? And for this:
dx=Sinv*(A*x+gmat+Fmat)
Would the order of multiplication of the inverse matrix with the stuff the in brackets not matter?
Alan Stevens
Alan Stevens le 16 Mai 2021
The tilde is just a placeholder, since your function doesn't make explicit use of t. This isn't essential, you could keep the t there if you prefer.
The order in which matrices and vectors are multiplied is very important!
King Hin Wong
King Hin Wong le 16 Mai 2021
If I had to do
(A*x+gmat+Fmat)*Sinv
For the sake of the order of multiplication, why does it return an error?
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);
King Hin Wong
King Hin Wong le 16 Mai 2021
Solved! Thank you so much, have a nice day :)

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur General Applications dans Centre d'aide et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by