# ODE45 function for 3 Variables

37 views (last 30 days)
Robsn on 30 Mar 2021
Edited: darova on 31 Mar 2021
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 on 30 Mar 2021
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);
.

James Tursa on 30 Mar 2021
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.

MG on 30 Mar 2021
Edited: MG on 30 Mar 2021
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]; )