Shouryya Ray's equation



Hi,
you might have seen a picture of the Indian teenager holding a poster 
with the equation he recently
discovered. If you are curious to know what it is about, here is a 
possible derivation of the equation.

The motion of a projectile on a vertical plane xy is due to two external 
forces, its weight [0, -m*g],
and the air resistance (Fr) which is proportional to the square of the 
speed v (v=sqrt(vx^2+vy^2) and
it is always opposite to the velocity [vx,vy]

    Fr/m = -c*v*[vx,vy]

where m is the mass and c is a constant which depends on geometric shape 
and size, and the air density.
We will also use the slope of the trajectory, s=tan(vy/vx).

Let us first declare the dependence of vx, vy and s on time t

(%i2) depends([vx,vy,s],t)$

The total force divided by the mass gives the acceleration vector and 
thus we have two differential
equations for vx and vy

(%o2) [vx(t),vy(t),s(t)]
(%i3) eq1: diff(vx,t) = -c*vx*sqrt(vx^2+vy^2);

(%o3) 'diff(vx,t,1) = -c*vx*sqrt(vy^2+vx^2)
(%i4) eq2: diff(vy,t) = -c*vy*sqrt(vx^2+vy^2)-g;

replacing vy by s in the equations, they become

(%i5) eq3: ev(eq1,vy=s*vx,ratsimp);

(%o5) 'diff(vx,t,1) = -c*sqrt(s^2+1)*vx*abs(vx)

(%i6) eq4: ev(eq2,vy=s*vx,diff,ratsimp);

(%o6) s*'diff(vx,t,1)+'diff(s,t,1)*vx = -c*s*sqrt(s^2+1)*vx*abs(vx)

and we can write them in a simpler form

(%i7) sol: first (solve([eq3,eq4],[diff(vx,t),diff(s,t)]));

(%o7) ['diff(vx,t,1) = -c*sqrt(s^2+1)*vx*abs(vx),
        'diff(s,t,1) = -g/vx]

The only horizontal force is the component of the air resistance, which 
is proportional to v^2 and,
therefore, vx will either be positive or negative, but not both. Let us 
assume that the projectile
was launched in the x-positive direction

(%i8) assume(vx>=0)$

We can know now obtain an ordinary differential equation (ODE)  for s and vx

(%i9) eq5: 'diff(s,vx) = ratsimp(rhs(sol[2])/rhs(eq3));

(%o9) 'diff(s,vx,1) = g/(c*sqrt(s^2+1)*vx^3)

which is a very simple separable ODE and its solution is

(%i10) eq6: ode2(eq5,s,vx);

(%o10) (c*asinh(s)+c*s*sqrt(s^2+1))/(2*g) = %c-1/(2*vx^2)
(%i11) eq7: expand(first(solve(eq6,%c)));

(%o11) %c = 1/(2*vx^2)+c*asinh(s)/(2*g)+c*s*sqrt(s^2+1)/(2*g)
(%i12) eq8: subst(sqrt(s^2+1)=sqrt(vx^2+vy^2)/vx,eq7);

(%o12) %c = c*s*sqrt(vy^2+vx^2)/(2*g*vx)
           +1/(2*vx^2)+c*asinh(s)/(2*g)
(%i13) eq9: subst(s=vy/vx,eq8);

(%o13) %c = c*asinh(vy/vx)/(2*g)+c*vy*sqrt(vy^2+vx^2)/(2*g*vx^2)+1/(2*vx^2)

this last result is the famous equation. If the initial speed was v0 
with an angle a0 above the horizontal,
we have

(%i14) ic1(eq9,vx=v0*cos(a0),s=tan(a0));

(%o14) (cos(a0)^2*c*v0^2*asinh(vy/(cos(a0)*v0))
  +c*vy*sqrt(vy^2+cos(a0)^2*v0^2)+g) /(2*cos(a0)^2*g*v0^2)
   = (c*vx^2*asinh(vy/vx)+c*vy*sqrt(vy^2+vx^2)+g)/(2*g*vx^2)

Regards,
Jaime