ODE45 function for 3 Variables

I have to solve these differential equations with the ode45 function. I only know how to transform the equation for only one variable, but how can i create a proper StateSpaceForm with 4 equations?

Would this be the right approach:

clear all
clear variables
close all
clc


tspan = 0:.1:3;
y0 = [0;0;0;0];

% zA = x(1) zA' = x(2) zR = x(3) zR' = x(4)

% zA'' = dA/mA * zR' + cA/mA * zR - cA/mA * zA - dA/mA * zA'
% zR'' = dA/mR * (zA' - zR') + cA/mR * (zA - zR) - cR/mR*(zR - zF)

% zA=x1 , zA'=x2 . zR=x3 , zR'=x4


[t,x]=ode45(@StateSpaceForm, tspan, y0);

plot(t,x(:,1))
plot(t,x(:,2))

function dF=StateSpaceForm(t,x)

mA=256; %kg
mR=31; %kg
cA=1000; 
cR=12800;
dA=100; %Ns/m
zF = 0;

dF(1,1) = x(3); 
dF(2,1) = x(4);
dF(3,1) = (dA/mA)* x(4) + (cA/mA) * x(3) - (cA/mA) * x(1) - (dA/mA) * x(2);
dF(4,1) = (dA/mR) * (x(2) - x(4)) + (cA/mR) * (x(1) - x(3)) - (cR/mR)* (x(3) - zF);

end

Now the problem is that i get a wrong plot. its all 0. I dont know why. Can someone pls help?

[Star Strider]

I always use the Symbolic Math Toolbox to check my derivations:

syms za(t) zr(t) ca cr ma da mr
Eqn = [diff(za); diff(zr); diff(za,2); diff(zr,2)] == [0 0 1 0; 0 0 0 1; -ca/ma ca/ma -da/ma da/ma; ca/mr -(ca+cr)/mr da/mr -da/mr] * [za; zr; diff(za); diff(zr)] + [0; 0; 0; cr/mr];
[VF, Subs] = odeToVectorField(Eqn);

.

Answer 1

[James Tursa]

This

y0 = [0;0;0;0];

and this

zF = 0;

means that your derivative function will evaluate to 0 exactly, so the solution is 0 exactly. I.e., there is nothing forcing the solution away from 0. You need to have something driving the solution away from 0 with a non-zero y0 or zF.

Answer 2

[MG]

Hi Robsn,

your StateSpaceForm function doesn't seem to exactly match the ODE in your problem. I believe it should read as follows (I here kept the order of the terms as in the problem, as I don't see a need to rearrange them; and my notation is x(1)=zA,x(2)=zB,x(3)=\dot{zA},x(4)=\dot{zA}):

dF(1,1) = x(3); 
dF(2,1) = x(4);
dF(3,1) = -(cA/mA)*x(1) +(cA/mA)      *x(2) -(dA/mA)*x(3) + (dA/mA)*x(4) +0         ;
dF(4,1) =  (cA/mR)*x(1) -(cA/mR+cR/mR)*x(2) +(dA/mR)*x(3) - (dA/mR)*x(4) +(cR/mR)*zF;

(ps, see also James reply that you would only find the trivial solution if you keep zF = 0, with the intial condition y0 = [0;0;0;0]; )