Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: 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?
Draft saved Draft deleted