Using Matlab to solve second order ODE's

  • Thread starter Vuldoraq
  • Start date
  • #1
272
0

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

Answers and Replies

Related Threads on Using Matlab to solve second order ODE's

  • Last Post
Replies
2
Views
910
Replies
3
Views
5K
Replies
5
Views
4K
  • Last Post
Replies
0
Views
952
Replies
2
Views
5K
Replies
16
Views
2K
Replies
1
Views
386
Replies
4
Views
2K
  • Last Post
Replies
4
Views
751
Top