I have been trying to get the numerical solution to the following system of ordinary differential equations:
Equations for the movement of a body through air in a inclined lunch:
(apparently LaTeX doesn't work on stack overflow)
u'= -F(u, theta, t)*cos(theta)
v'= -F(v, theta, t)*sin(theta)-mg
by the Runge-Kutta-Fehlberg Algorithm, but in the middle of the computation i have to calculate theta, that is calculated by
arccos(u/sqrt(u^2+v^2)) or arcsin(v/sqrt(u^2+v^2)),
but eventually theta gets too small and I need it to solve the function F( v, theta, t) and to find the value V = sqrt(v^2 + u^2) I use V = (v/sin(theta)), but as theta gets small so does sin(theta) and I get a numerical error from a given iteration forward -1.IND00, It is problably because theta is too small, i tried to make theta go from a small positive angle like 0.00001 to a small negative angle like -0.00001 (if(fabs(theta)<0.00001) theta = -0.00001) but it seems that theta gets trapped into this negative value, does anyone have an indication on what to do to resolve this numerical instability ?
It is a bad idea to use the inverse cosine or sine functions to determine the angle of a point. To get
theta = arg ( u + i*v)
use
theta = atan2(v,u).
This still has the problem that it jumps on the negative half axis, that is for v=0, u<0. That can be solved by making theta a third dynamical variable, so that
theta' = Im( (u'+i*v')/(u+i*v) ) = (u*v' - u'*v) / (u^2+v^2)
But really, the equation for the free fall with air friction is easiest implemented as
def friction(vx, vy):
v = hypot(vx, vy)
return k*v
def freefall_ode(t, u):
rx, ry, vx, vy = u
f=friction(vx, vy)
ax = -f*vx
ay = -f*vy - g
return array([ vx, vy, ax, ay ])
so that you do not need any angle or to try to weaken the coupling of the velocity components by reducing it to the angle of the velocity vector. This you can now plug into the integration method of your choice, applied as a method for vector-valued systems.