CODEWITHSUNDEEP
matlabequationode
Queen Olang
Queen Olang

Reputation: 1

Solving a system of second order differential equations using MATLAB

I am trying to solve a projectile motion problem, to determine the take off velocity under given initial conditions, the problem was reduced to a system of two second order differential equations. My code and question is in the pictures below. The values of constants in the problem equations have been reduced to constants a, b, c and d.

x¨(t)=-1/2m ρAC_d cos⁡(arctan⁡((y˙(t))/(x˙(t) )))(〖x˙(t)〗^2+ 〖y˙(t)〗^2)
y¨(t)=-1/2m(2mg+ρAC_d sin⁡(arctan⁡((y˙(t))/(x˙(t) )))(〖x˙(t)〗^2+ 〖y˙(t)〗^2)

# With the initial conditions:

x˙(0)=cosθ ∙ V_0

y˙(0)=sinθ ∙ V_0

x(0)=0

y(0)=0

My solution code is shown below;

syms x(t) y(t) a b c d u theta
% Equations
% d2x = -a*cos(arctan(dy/dx))*(((dx)^2)+(dy)^2));
% d2y = -b*(c + d*sin(arctan(dy/dx))*(((dx)^2)+(dy)^2));

%Constants
dx=diff(x,t);
dy=diff(y,t);
d2x=diff(x,t,2);
d2y=diff(y,t,2);
a=1;
b=2;
c=3;
d=4;

%Initial Conditions

cond1 = x(0) == 0;
cond2 = y(0) == 0;
cond3 = dx(0) == u*cos(theta);
cond4 = dy(0) == u*sin(theta);

conds = [cond1 cond2 cond3 cond4];

eq1 = -a*cos(atan(dy/dx))*(((dx)^2)+(dy)^2);
eq2 = -b*(c + d*sin(atan(dy/dx))*(((dx)^2)+(dy)^2));

vars = [x(t); y(t)];
V = odeToVectorField([eq1,eq2]);
M = matlabFunction(V,'vars', {'t','Y'});
interval = [0 5];  %time interval    
ySol = ode23(M,interval,conds);

The error messages are shown below;

Error using mupadengine/feval (line 187)
System contains a nonlinear equation in 'diff(y(t), t)'. The system must be quasi-linear:
highest derivatives must enter the differential equations linearly.

Error in odeToVectorField>mupadOdeToVectorField (line 189)
T = feval(symengine,'symobj::odeToVectorField',sys,x,stringInput);

Error in odeToVectorField (line 138)
sol = mupadOdeToVectorField(varargin);

Error in velocity_takeoff (line 29)
V = odeToVectorField([eq1,eq2]);

Why do I get and how can I alleviate these errors?

Upvotes: 0

Views: 570

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25972

That looks sub-optimal. There is no need for trigonometric functions, especially as in the form used they may introduce sign errors, use atan2 to avoid that. The equations are in a simplified form

as vectors dv/dt = -g*e_2 - k*speed*v, where speed = |v| is the norm of the vector
and additionally dxy/dt=v,  xy=[x,y] being the position vector

Implemented in components, this gives

function du_dt = motion(t,u)
  x=u(1); y=u(2); vx=u(3); vy=u(4);
  speed = hypot(vx,vy);
  ax = -k*speed*vx;
  ay = -k*speed*vy - g;
  du_dt = [vx; vy; ax; ay];
end%function

This direct implementation is shorter and better readable than the way via symbolic equations.

You can adapt this for the different constants used, but any change in the constant structure is unphysical, or would need to be justified with some uncommon physical arguments.

Upvotes: 1

Related Questions