Finding derivatives of vector function for numerical integration

[graffy76]

Hello,

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]

where:

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

where:

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.

Thanks.

Joel

 

[jvicens]

Hello Joel,

Thank you for sharing your math problem and explanation. It seems like you have put a lot of thought and effort into this problem and your solution looks sound to me. However, I would recommend double checking your calculations and equations to make sure there are no errors. Also, if possible, it would be beneficial to run some simulations to test your solution and make sure it produces accurate results. It's always good to have some real-world data to validate your mathematical model.

Additionally, I would suggest seeking feedback from other experts in the field of molecular dynamics or numerical methods. They may be able to provide some insights or suggestions on your approach.

Overall, it seems like you have a good understanding of the concepts and equations involved in solving this problem. Keep up the good work and don't hesitate to reach out for help or clarification if needed. Good luck!