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:

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


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