Help solving using ode45

Hello everyone.

I need help solving a system of 4 first order differential equations shown below, to find the displacements x1 and y1 as well as velocities x1_dot and x2_dot. Time span is 0 to 20 seconds. Initial conditions should be zero.

x1==>displacement, x1_dot==>velocity, x2==>velocity, x2_dot==>acceleration

Where k=8500, m1=2650, m2=1830

Both Fd and Fb are functions of time and have different values for different intervals.

Fd;

Fb;

I think Fd and Fb are going to need seperate functions. Any help to solve this is much appreciated!

Thank you.

[Bruno Luong]

Many common aspects with this recent thread: Link

[Nadun Chandrabahu]

Thanks for your reply,

I'm mostly confused with how to deal with the y1 in the 2nd equation and the x1 in the 4th equation, since these states are not in their proper equations.

I know its going to be long, but if you could help me write this code I would really appreciate it.

[Bruno Luong]

Sorry, I'm not going to write code for you.

However if you make an effort to write something (you might adapt the code in this thread) we can take a look.

[Nadun Chandrabahu]

I'll see about it thanks.

[Nadun Chandrabahu]

Hello, I did some coding. But I'm getting the following error messages: Could you help me on how to go from here?

Below is my code,

%main script
clear all;
clc;

global x y xdot ydot t;
global k m1 m2;

k=8500;
m1=2650;
m2=1830;

%xdot(1)=x(2); %trying to write the 1st & 3rd DE's not sure if this is correct
%ydot(1)=y(2);

[t,x] = ode45('func_4wd', [0, 20], [0, 0]);
[t,y] = ode45('func_car', [0, 20], [0, 0]);

plot(t,x);
%plot(t,x);

%function (2nd DE)
function xdot = func_4wd(t,x,y)
global m1;
m1=2650;

xdot = [x(2); (DrivingF(t)/m1)-Tension(x(1),y(1))-(140/m1)*x(2)*abs(x(2))];

%function (4th DE)
function ydot = func_car(t,x,y)
global m2;
m2=1830;
ydot = [y(2); Tension(x(1),y(1))-(BrakingF(t)/m2)-((925/m2)*y(2)*abs(y(2)))];

%function for tension
function TensionOutput=Tension(X1,Y1)

if (X1-Y1)<3.5
    TensionOutput=0;
else
    TensionOutput=8500*(X1-Y1-3.5);
end

end

%function for driving force
function DrivingOutput = DrivingF(t)

     DrivingOutput=0;

     if (t>0)&&(t<=4);
         DrivingOutput=(7500/197)*(150-150*cos(pi*t/4));
     end
     if (t>=4)&&(t<=12);
         DrivingOutput=(2250000/197);
     end
     if (t>=12)&&(t<=16);
         DrivingOutput=(7500/197)*(150-150*cos(pi*t/4));
     end
  end

%function for braking force
function BrakingOutput = BrakingF(t)

BrakingOutput=0;

     if (t>=6)&&(t<=7);
         BrakingOutput=(4000*t-24000);
     end
     if (t>=7)&&(t<=8);
         BrakingOutput=(32000-4000*t);
     end
     if (t>=13)&&(t<=14);
         BrakingOutput=(4000*t-52000);
     end   
     if (t>=14)&&(t<=15);
         BrakingOutput=(60000-4000*t);
     end
   end

Answer 1

[Bruno Luong]

%main script
close all
clear

k=8500;
m1=2650;
m2=1830;

[t,xy] = ode45(@(t,xy) xyder(t,xy,k,m1,m2), [0, 20], [0, 0, 0, 0]);

x = xy(:,1);
y = xy(:,3);
hold on
plot(t,x);
plot(t,y);
legend('x','y');
xlabel('t')

function xyd = xyder(t,xy,k,m1,m2)
x1 = xy(1);
x2 = xy(2);
y1 = xy(3);
y2 = xy(4);
x2d = DrivingF(t)/m1 - Tension(x1,y1,k,m1) - (140/m1)*x2*abs(x2);
y2d = Tension(x1,y1,k,m2) - BrakingF(t)/m2 -(925/m2)*y2*abs(y2);
xyd = [x2; x2d; y2; y2d];
end

%function for tension
function TensionOutput=Tension(x,y,k,m)

TensionOutput=k*(x-y-3.5)/m;
% if (x-y)<3.5 % I don't see this IF in your note
%     TensionOutput=0;
% else
% TensionOutput=k*(x-y-3.5);
%end
end

%function for driving force
function DrivingOutput = DrivingF(t)
if t<4
    DrivingOutput=(7500/197)*(150-150*cos(pi*t/4));
elseif t<12
    DrivingOutput=(2.25e6/197);
elseif t<=16
    DrivingOutput=(7500/197)*(150-150*cos(pi*t/4));
else
    DrivingOutput=0;
end
end

%function for braking force
function BrakingOutput = BrakingF(t)
if (t>=6)
    if t<7
        BrakingOutput=(4000*t-24000);
    elseif t<8
        BrakingOutput=(32000-4000*t);
    elseif t<14
        BrakingOutput=(4000*t-52000);
    elseif t<15
        BrakingOutput=(60000-4000*t);
    else % t>=15
        BrakingOutput=0;
    end
else % % t<6
    BrakingOutput=0;
end
end

[Nadun Chandrabahu]

Great, this works. Thank you so much!!