(adsbygoogle = window.adsbygoogle || []).push({}); 1. The problem statement, all variables and given/known data

A projectile is ﬁred with an initial speed of 1000 m/s at an angle α from the horizontal. A

good model for the path uses an air-resistance force proportional to the square of the speed.

This leads to the following equations for the acceleration in the x and y directions:

[tex]\ddot{x}=-c*\dot{x}*(\dot{x}^{2}+\dot{y}^{2})^{1/2}[/tex]

[tex]\ddot{y}=-c*\dot{y}*(\dot{x}^{2}+\dot{y}^{2})^{1/2}-g[/tex]

where an over dot signiﬁes a derivative with respect to time, and the constant c can be taken to be c = 0.005. The initial conditions are:

[tex] x(0) =y(0) =0\ and\ (\dot{x(0)},\dot{y(0)}) =1000(cos\alpha,sin\alpha)[/tex]

(b) Convert the above system of equations to 4 ﬁrst-order equations.

(c) Solve the system of equations over the time interval [0,30] (with e.g. the RK4 method)

for each launch angle α =10:10:80 degrees (convert to radians!). Plot the resulting

trajectories.

2. Relevant equations

None

3. The attempt at a solution

I transformed the second order ode's to four first order;

[tex]\frac{dx}{dt}=x_{1}[/tex]

[tex]\frac{dx_{1}}{dt}=-c*x_{1}*(x_{1}^{2}+y_{1}^{2})^{1/2}[/tex]

[tex]\frac{dy}{dt}=y_{1}[/tex]

[tex]\frac{dy_{1}}{dt}=-c*y_{1}*(x_{1}^{2}+y_{1}^{2})^{1/2}-g[/tex]

And wrote this program to represent them,

function [dx]=secondorder(t,v)

c=0.005;

g=9.8;

x=v(1);

y=v(2);

dx(1)=x;

dx(2)=-c.*x.*sqrt(x.^2+y.^2);

dx(3)=y;

dx(4)=-c.*y.*sqrt(x.^2+y.^2)-g;

Then fed this program into this fourth order Runge-Kutte program,

function [sol]=rk4_syst(fcn,a,b,y0,N) % fcn now outputs the vector of derivs.

h=(b-a)/N; % and y0 is the vector if initial values

y(1,:) = y0; % copy into the first vector-point

x=a+(0:N)*h;

for k =1:N % loop over steps as before

k1 = feval(fcn,x(k),y(k,; % for each step, compute

k2 = feval(fcn,x(k)+h/2,y(k,:)+h*(k1)/2); % the slopes k1,...,k4

k3 = feval(fcn,x(k)+h/2,y(k,:)+h*k2/2); % for all the equations i

k4 = feval(fcn,x(k)+h,y(k,:)+h*k3); % by using vector-notation

y(k+1,:)=y(k,:)+h*(k1+2*(k2+k3)+k4)/6; % then take vector step

end % (simultaneously for all i)

sol=[x',y];

But everytime I try to solve my system I get silly results, for instance the y part goes way to high or the x part stays fixed! I've spent ages trying to puzzle it out to no avail. Please could someone give me a hand? I know it's an involved question, but I've run of palces to turn.

**Physics Forums | Science Articles, Homework Help, Discussion**

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Using Matlab to solve second order ODE's

Can you offer guidance or do you also need help?

**Physics Forums | Science Articles, Homework Help, Discussion**