Hello,
At first mpmath seemed to give me right results for numerical
differentiation of some transcendental equations.
These functions are related to the natural vibration modes of a fixed-
fixed beam.
Normally I can differentiate analytically to get the slope and
inflection points. While trying to have a quick look into these points
with numerical differentiation I noticed that mpmath is giving me
values quite distant from the ones obtained with analytical formulas.
Am I using mpmath.diff in the wrong way? I tried to modify some
optional diff switches without success...
In the bottom I placed a simplified script to display the issue.
By the way I run it on OSX 10.6.8 with a binary packaged sage-4.7.1.
Any idea is appreciated.
Best regards,
Guilherme
----------
import mpmath
# first root of: cos(k)*cosh(k)=1
k = mpmath.mpf(4.730040744862704026024048100833884819898342)
alpha = ( mpmath.cos(k) - mpmath.cosh(k) ) / ( mpmath.sinh(k) -
mpmath.sin(k) )
# fundamental mode
def phi(x):
return mpmath.cos(k*x) - mpmath.cosh(k*x) + alpha*
( mpmath.sin(k*x) - mpmath.sinh(k*x) )
# slope of fundamental
def phi_dx(x):
return -mpmath.sin(k*x) - mpmath.sinh(k*x) + alpha*
( mpmath.cos(k*x) - mpmath.cosh(k*x) )
# curvature of fundamental
def phi_dx2(x):
return -mpmath.cos(k*x) - mpmath.cosh(k*x) - alpha*
( mpmath.sin(k*x) + mpmath.sinh(k*x) )
### Why the next functions doesn't match the analytical ones on the y-
axis? ###
# numerical slope of fundamental
def phi_dx_mp(x):
return mpmath.diff(phi, x, 1)
# numerical curvature of fundamental
def phi_dx2_mp(x):
return mpmath.diff(phi, x, 2)
# Plot the functions
import numpy as np
xx = np.linspace(0,1,100)
leg=['Fundamental','Slope','Curvature']
a1=line(([x,phi(x)] for x in xx), legend_label=leg[0],
rgbcolor=(0,0,1))
a2=line(([x,phi_dx(x)] for x in xx), legend_label=leg[1],
rgbcolor=(0,1,0))
a3=line(([x,phi_dx2(x)] for x in xx), legend_label=leg[2],
rgbcolor=(1,0,0))
g=plot(a1+a2+a3); g.show()
### See the differences on the y-axis for the numerically
differentiated ones...
m2=line(([x,phi_dx_mp(x)] for x in xx), legend_label=leg[1]+'
mpmath', rgbcolor=(0,1,0)); m2.show()
m3=line(([x,phi_dx2_mp(x)] for x in xx), legend_label=leg[2]+'
mpmath', rgbcolor=(1,0,0)); m3.show()
----------
--
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
