How do i transform a second grade ODE into coordinates?

Hello,

I'm new to matlab and i'm quite struggling with it.

I have been assigned a coursework in which i must find the trajectory of a projectile.Within it there is a second order ODE which i must transform into x&y's in order to get a system of first-order equations.

The second order equation is: r''(t) = (1 /m)* (((− 1 /2 )*ρ*cd*A*|r'|*r' )− mg );

Mathematically, this system will be expressed a:

y' = F(t, y(t))=dy/dt,

where: y(t) = [vx(t), vy(t), rx(t), ry(t)]

Is there any way to accomplish that in matlab and could anyone give suggestions with examples on how to do so?

Thanks!

Answer 1

[James Tursa]

Looks like motion is restricted to be on the x-y plane, so r is 2D. You have a 2nd order ODE in a 2D frame, so that means you have a 2 * 2 = 4th order system of 1st order equations. Define the following:

r = [x;y]

From that define a state vector consisting of all the variables and derivatives up to but not including the highest order:

Y = [x;y;xdot;ydot]

Then the derivative of your state vector is

Ydot = d/dt [x;y;xdot;ydot] = [xdot;ydot;xdotdot;ydotdot] = [Y(3);Y(4);xdotdot;ydotdot]

So then you simply need to get the expressions for xdotdot and ydotdot. You do this by rewriting your rdotdot equation above in terms of x and y using the definition r = [x;y]. You do that by replacing the following

r' = d/dt r = d/dt [x;y] = [xdot;ydot] = [Y(3);Y(4)]

Thus

|r'| = norm([Y(3);Y(4)])

That being said, you are missing something on the -mg term. This is a vector equation so you need to attach a vector direction for this. I would assume this should be in the y direction (vertical).

So write a derivative function or function handle that takes in a 4-element Y vector and outputs the Ydot vector.

[Nicolae Lungu]

Thank you very much for helping me

Actually i have contacted the proffesor so could give some insights on how to solve it.

And so he said: i need to take r''(t) = (1/m)*(((− 1/2)* ρcdA|r'|r') − mg ). equation (1)! and project it into x and y, which then is gonna give the equations for y' = F(t, y(t)), in x and y directions. Which then will have y(t) = [vx(t), vy(t), rx(t), ry(t)].

Can you please show me how can this be incorporated into matlab language?, because i don't quite understand the logic how to implement it...

[James Tursa]

Just plugging into what I have already written above:

rdot = [Y(3);Y(4)];
rdotdot = (1/m) * (((-1/2)*p*cd*A*norm(rdot)*rdot)-m*[0;g]);
Ydot = [rdot;rdotdot];

[Nicolae Lungu]

thank you very much

[Nicolae Lungu]

Ok, so i have defined the function ,written below, but when tryng to find out if it works i get the folowing error:

function Ydot = drag_ode(t, y)
m= 0.005;          % mass of the projectile
g= 9.81;          % acceleration due to gravity constant
cd= 0.479;         % Drag coefficient
D= 0.1;            % Diametet m
A=pi*(D/2)^2;     % Area of sphere, A=0.0314 (m)
% equations of motion for projectile 
Y = [t;y;tdot;ydot];
Ydot = [Y(3);Y(4);tdotdot;ydotdot];
rdot = [Y(3);Y(4)];
rdotdot = (1/m)*(((-1/2)*p*cd*A*norm(rdot)*rdot)-m*[0;g]);
Ydot = [rdot;rdotdot];
end

Not enough input arguments.

->Y = [t;y;tdot;ydot];

i dont quite understand what is the problem......

[James Tursa]

One fix would be to replace this

Y = [t;y;tdot;ydot];
Ydot = [Y(3);Y(4);tdotdot;ydotdot];

with this

Y = y;

[Nicolae Lungu]

gives me the same error

[Steven Lord]

I interpreted this section of James's post:

Define the following:

r = [x;y]

From that define a state vector consisting of all the variables and derivatives up to but not including the highest order:

Y = [x;y;xdot;ydot]

as describing the mathematics behind convering the second-order ODE into a system of first-order ODEs, not as describing the code you'll need to write to solve the system. So your function could start off something more like:

function Ydot = drag_ode(t, y)
% y is [x coordinate; y coordinate; dx/dt; dy/dt]

Then when you define the vector of left-hand sides you can fill in the first one and accompany it with a comment:

Ydot = zeros(size(y));
Ydot(1) = y(3); % d (x coordinate)/dt = dx/dt

IMO the comment makes it clearer how Ydot(1) and y(3) are related. Obviously y(3) is dx/dt since that's what the comment at the top of the function says. You'd then continue on to fill in Ydot(2) through Ydot(4) the same way.

[James Tursa]

@Nicolae: Please post all of your current code for the drag_ode function, as well as the complete code that is used to call the drag_ode function. You have to have code that calls your drag_ode function ... you can't simply run this function directly without giving it any inputs or you will get an error. Also show us the complete error message you are getting including the offending line.

[Nicolae Lungu]

ok, so at the moment the code looks like this:

function Ydot = drag_ode(t, y)
% y is [x coordinate; y coordinate; dx/dt; dy/dt]

m= 0.005;          % mass of the projectile
g= 9.81;          % acceleration due to gravity constant
cd= 0.479;         % Drag coefficient
D= 0.1;            % Diametet m
A=pi*(D/2)^2;     % Area of sphere, A=0.0314 (m)
y = [t;y;tdot;ydot];

Ydot = zeros(4,1);
Ydot(1) = y(3); % d (x coordinate)/dt = dx/dt

Ydot = [y(3);y(4);tdotdot;ydotdot];
rdot = [y(3);y(4)];
rdotdot = (1/m)*(((-1/2)*p*cd*A*norm(rdot)*rdot)-m*[0;g]);
Ydot = [rdot;rdotdot];

end

And this is the function calling the drag_ode:

function [rx,ry,vx,vy]= solve_ode45(v0,h,dt,theta)
g = 9.81  ; % acceleration due to gravity at Earth surface(m/s^2)
cD = .479  ; % coefficient of drag
m  = 5.5  ; % mass(kg)
v0 = 150;                 % the initial velocity
t(1) = 0;                 % the initial time
tF = 20;                  % the end time
dt = .001                 %time step
theta=45
tspan = (0:dt:25);

IC=[0; v0*cos(theta); h; v0*sin(theta)]; %IC - initial coordinates
[t, oput] = ode45(@drag_ode, tspan, IC);

rx = oput(:,1);%Position coordinates for x axis
vx = oput(:,2);%Velocity coordinates for x axis
ry = oput(:,3);%Position coordinates for y axis
vy = oput(:,4);%Velocity coordinates for y axis

figure(1); clf;
plot(x,y);
xlabel('Distance (m) ');
ylabel('Height (m) ');
title('Motion of projectile');
end

When trying to run te program ,this is the message stating the error in the command window:

solve_ode45(150, 5, 25, 45)

dt =

0.001

theta =

45

Undefined function or variable 'tdot'.

Error in drag_ode (line 9)

y = [t;y;tdot;ydot];

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in solve_ode45 (line 13)

[t, oput] = ode45(@drag_ode, tspan, IC);

[James Tursa]

Delete all of this code completely:

y = [t;y;tdot;ydot];
Ydot = zeros(4,1);
Ydot(1) = y(3); % d (x coordinate)/dt = dx/dt
Ydot = [y(3);y(4);tdotdot;ydotdot];

And change this:

C=[0; v0*cos(theta); h; v0*sin(theta)];

to this:

C=[0; h; v0*cosd(theta); v0*sind(theta)];

You will also need to define the p value (i.e., rho) in your derivative function.

[Nicolae Lungu]

ok now the code looks like this drag_ode:

function Ydot = drag_ode(t, y)

% y is [x coordinate; y coordinate; dx/dt; dy/dt]
m= 0.005;          % mass of the projectile
g= 9.81;          % acceleration due to gravity constant
cd= 0.479;         % Drag coefficient
D= 0.1;            % Diametet m
A=pi*(D/2)^2;     % Area of sphere, A=0.0314 (m)
rho = 1.225;      % density of atmosphere near Earth surface(kg/m^3)
%y = [t;y;tdot;ydot];
% Ydot = [y(3);y(4);tdotdot;ydotdot];
%Ydot = zeros(4,1);
ydot(1) = y(3); % d (x coordinate)/dt = dx/dt
ydot(2) = y(4); % d (y coordinate)/dt = dy/dt
ydot(3) = tdotdot;   
ydot(4) = ydotdot;
rdot = [ydot(3);ydot(4)];
rdotdot = (1/m)*(((-1/2)*rho*cd*A*norm(rdot)*rdot)-m*[0;g]);
Ydot = [rdot;rdotdot];

end

and like this for solve_ode45:

function [rx,ry,vx,vy]= solve_ode45(v0,h,dt,theta)
g = 9.81;                 % acceleration due to gravity at Earth surface(m/s^2)
cD = .479;                % coefficient of drag
m  = 5.5;                 % mass(kg)
v0 = 150;                 % the initial velocity
t(1) = 0;                 % the initial time
tF = 20;                  % the end time
dt = 0.001;
theta = 45;
tspan = (0:dt:25);

IC=[0; h; v0*cosd(theta); v0*sind(theta)];% IC-the initial coordinates
[t(1), oput] = ode45(@drag_ode, tspan, IC);

rx = oput(:,1);           %Position coordinates for x axis
vx = oput(:,2);           %Velocity coordinates for x axis
ry = oput(:,3);           %Position coordinates for y axis
vy = oput(:,4);           %Velocity coordinates for y axis

figure(1); clf;
plot(x,y);
xlabel('Distance (m) ');
ylabel('Height (m) ');
title('Motion of projectile');
end

and i still get the error:

solve_ode45(150, 5, 25, 45)

Undefined function or variable 'tdotdot'.

Error in drag_ode (line 15)

ydot(3) = tdotdot;

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in solve_ode45 (line 13)

[t(1), oput] = ode45(@drag_ode, tspan, IC);

[James Tursa]

Why do you keep putting code in where I tell you to delete it? Why do you keep using variables such as tdotdot and ydotdot that don't exist?

Delete this code completely:

ydot(1) = y(3); % d (x coordinate)/dt = dx/dt
ydot(2) = y(4); % d (y coordinate)/dt = dy/dt
ydot(3) = tdotdot;   
ydot(4) = ydotdot;

and put this line back the way it was:

rdot = [y(3);y(4)];

And then the output is in x,y,xdot,ydot order, so change to this:

rx = oput(:,1);           %Position coordinates for x axis
ry = oput(:,2);           %Position coordinates for y axis
vx = oput(:,3);           %Velocity coordinates for x axis
vy = oput(:,4);           %Velocity coordinates for y axis