- #1
Vuldoraq
- 272
- 1
Homework Statement
A projectile is fired 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 signifies 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 first-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.