Dismiss Notice
Join Physics Forums Today!
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

  1. Nov 10, 2008 #1
    1. The problem statement, all variables and given/known data

    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:


    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

    2. Relevant equations


    3. The attempt at a solution

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





    And wrote this program to represent them,

    function [dx]=secondorder(t,v)



    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

    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)


    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.:cry:
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?