GDS-Education-Community-of-Practice / DSECOP

This repository contains data science educational materials developed by DSECOP Fellows.
Creative Commons Zero v1.0 Universal
43 stars 25 forks source link

Improper formulation of the Velocity Verlet algorithm in 01-differential-equations notebook #15

Open marcobn opened 1 year ago

marcobn commented 1 year ago

Dear Developers, First of all, I want to thank you for this excellent collections of modules. They are great tools to educate students to data science. If I may, I would like to point your attention on an issue that was actually raised initially by my students in an advanced computational physics class. We have been working on the first module of the section on solving differential equations with NN and we noticed an inconsistency in the suggested code for the Velocity Verlet algorithm. In cell 4 the VV algorithm is implemented as:

for i in range(n-1):

expression for acceleration

    a[i] = g - (b/m)*v[i]
    # update position
    y_vv[i+1] = y_vv[i] + v[i]*DeltaT + 0.5*a[i]*DeltaT**2
    # updated expression for acceleration
    a[i+1] = g - (b/m)*v[i+1]
    # update velocity
    v[i+1] = v[i] + 0.5*DeltaT*(a[i+1] + a[i])
    # update time to next time step and compute analytical answer
    t[i+1] = t[i] + DeltaT

The issue is that in line 7 we attempt to update a using a velocity that we have not yet calculated (v[i+1] is zero, given the way arrays are initialized). The problem is that in the usual formulation of the VV we assume the acceleration to be constant (trivial case) or to depend on position (as in a potential field) and not velocity. If a depends on v, as in the present case, the update of the velocity should be done considering that a[i+1] = g-b/mv[i+1], plug this in the formula for the updated velocity and solve for v[i+1]. This results in v[i+1] = (v[i] + (g-0.5b/mv[i])DeltaT)/(1+0.5b/mDeltaT)

Running the original notebook the RMSE Between Velocity-Verlet Method and the Exact Solution is 0.3026318710261818. With the implementation above the RMSE is 9.351368107168857e-05, in line with the expected DeltaT**3 scaling.

I hope you will find this useful. Best,

marco

Marco Buongiorno Nardelli, (he/him) Regents Professor, University of North Texas External Professor, Santa Fe Institute CEMI, Center for Experimental Music and Intermedia iARTA, Initiative for Advanced Research in Technology and the Arts ArtSciLab, ATEC @ the University of Texas at Dallas

http://ermes.unt.edu/ http://www.materialssoundmusic.com/

soltaniehha commented 1 year ago

Thank you @marcobn for bringing this to our attention. @juliebutlerhartley, could you please have a look at Marc's comment above?