- #1

- 272

- 0

## Homework Statement

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.

## Homework Equations

None

## 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.