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

Finding derivatives of vector function for numerical integration

  1. Jul 24, 2008 #1

    I have a math problem that I think I've worked out properly, but I'm not entirely sure. The explanation is a bit lengthy, but I don't want to miss anything that might be pertinent.

    Essentially, I have a force equation F(t) that describes the acceleration of a body in two dimensions (F=ma, assuming unit mass).

    The body position, velocity, and acceleration are then calculated using Gear's predictor-corrector method (a method that's popular in molecular dynamics circles, as I understand it).

    The Gear predictor portion is a fifth order Taylor expansion on Newton's equations for curvlinear motion:

    x[n+1] = x[n] + h*v[n] + ((h^2)/2)a[n] + ((h^3)/6)B[n] + ((h^4)/24)C[n]
    v[n+1] = v[n] + h*a[n] + ((h^2)/2)B[n] + ((h^3)/6)C[n]
    a[n+1] = a[n] + h*B[n] + ((h^2)/2)C[n]
    B[n+1] = B[n] + h*C[n]


    B[n] and C[n] are the first and second derivatives of the acceleration (F(t)).
    [n] indicates the value calculated in the nth timestep
    [n+1] indicates the predictor value (calculated before the next timestep).
    h is the timestep width (in my case, 0.1 seconds).

    My force equation:

    F(t) = (v * e(t)) / tau


    v = velocity scalar constant
    tau = time constant

    e(t) vector indicating direction of travel from a current position vector (x) to a desired position vector (p).

    e(t) = (p - x) / ||p - x||

    The trick is to convert an instantaneous vector function to one that is defined with respect to time (t) or in this case, timestep (h).

    The current position (x) is the only part of e(h) that varies with time. Since it describes displacement, I chose to describe it using Newton's equations thusly:

    x(h) = x[n] + hv[n] + ((h^2)/2)a[n]

    Substituting x(h) in for x in e(h) results in a vector function that looks like this (expanded):

    e(h) (numerator):
    ( p[x] - x[n] - hv[x][n] - ((h^2)/2)a[x][n] ) i +
    ( p[y] - y[n] - hv[y][n] - ((h^2)/2a[y][n] ) j

    e(h) (denominator):
    sqrt( (( p[x] - x[n] - hv[x][n] - ((h^2)/2)a[x][n] ) ^2 ) i +
    ((p[y] - y[n] - hv[y][n] - ((h^2)/2)a[y][n] ) ^2 ) j )

    the first derivative, (dF(h) / dh) = B[h]:

    B[h] (numerator):
    [ ( -v[x][n] - h*a[x][n] )i + ( -v[y][n] - h*a[y][n] ) j ] * (v / tau)

    B[h] (denominator):
    sqrt ( (( -v[x][n] - h* a[x][n] ) ^2 ) i +
    (( -v[y][n] - h*a[y][n] ) ^2 ) j )

    the second derivative, (d2F(h) / dh2) = C[h]:

    C[h] (numerator):
    [ (-a[x][n]) i + (-a[y][n]) j ] * (v / tau)

    C[h] (denominator):
    sqrt ( (-a[x][n]^2) i + (-a[y][n]^2) j )

    The big trick was figuring out how to convert an instantaneously-derived vector function to one that is defined with respect to time, which is why x became x(t) (or x(h)). Having done that, I simply want to make sure that finding the first and second derivatives was as straight-forward as it appeared to be.

    I'm pretty sure this is how to do it - I don't see really any other way to get what I need to implement a fifth-order Gear solution. However, I must admit Calculus and Physics were learned in excess of ten years ago for me, and they're hardly necessary disciplines in my current line of work, so I'm a bit rusty.

    Any comments or thoughts? Any help is greatly appreciated.


  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?