Pour être précis, ce genre de solution consiste par approcher la dérivée analytique, exacte donc, par une dérivée numérique. Si tu veux augmenter la précision numérique de ce genre de dérivée, tu peux utiliser la technique du Romberg qui permet d'éliminer les contaminations des ordres supérieurs:
f'(i,j) = (4**j*f'(i,j-1)-f'(i+1,j-1))/(4**j-1)
avec
f'(i,0)=(f(x+dx)-f(x-dx))/(2*dx) et epsilon=2**i*epsilon
En pratique, des premières estimations de la dérivée numérique sont calculée pour des "dx" croissant (i=0,1,2,...): f'(0,0),f'(1,0),f'(2,0). Ces valeurs sont ensuite combinées entre elles suivant la première équation: f'(0,1)=(4*f'(0,0)-f'(1,0))/3, f'(1,1)=(4*f'(1,0)-f'(2,0))/3, f'(0,2)=(16*f'(0,1)-f'(1,1))/15, ...
Une analyse du triangle de valeurs ainsi obtenu te permet en outre de savoir si l'epsilon choisi était trop petit (problème de précision numérique due à la machine) ou trop grand (encore trop de contamination des ordres supérieurs).
Une solution procurant la valeur exacte est d'employer les nombres de Rall. D'un point de vue informatique, la valeur de la fonction en un point, f(x0,y0,...), est stockée ainsi que l'ensemble des dérivées premières par rapport aux différents paramètres, (df/dx(x0),df/dy(y0),...). La dérivation d'un fonction quelconque est effectuée par récurrence en suivant les règles mathématiques de dérivations pour l'addition, la multiplication, ... (voir par exemple: http://www.american.edu/academic.d [...] todiff.pdf, ou faire des recherches sur google avec les mots clés: Rall number differentiation). Si jamais cela t'intéresse, j'ai réalisé le code C++ de ce genre de nombre.