0

I want to calculate derivative of a function of two variables

f(x,y) = exp(sin(sqrt(x^2+y^2)))

[which for 1D case reduces to

f(x) = exp(sin(x))

well covered under this post Finding first derivative using DFT in Python ] using fourier transforms.

But, I am getting wrong result compared to analytical derivative of this function w.r.t to x and y variable. Below is the used code:

import numpy as np

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt


def surface_plot3D(X,Y,Z):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                           linewidth=0, antialiased=False)
    ax.set_zlim(-3., 3.)

    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

    fig.colorbar(surf, shrink = 0.5, aspect = 5)

def funct(x,y):
    return np.exp(np.sin( np.sqrt(x**2 + y**2) ))    


N = 2**4
xmin=0
xmax=2.0*np.pi
step = (xmax-xmin)/(N)
xdata = np.linspace(step, xmax, N)

ydata = np.copy(xdata)
X,Y = np.meshgrid(xdata,ydata)
Z = funct(X,Y)

surface_plot3D(X,Y,Z)
plt.show()

vhat = np.fft.fft2(Z)

what = 1j * np.fft.fftfreq(N,1./N)
whax, whay = np.meshgrid(what, what)

wha_x = whax * vhat  
w_x = np.real(np.fft.ifft2(wha_x)) 
wha_y = whay * vhat  
w_y = np.real(np.fft.ifft2(wha_y)) 
w_tot = np.sqrt( w_x**2 + w_y**2)
surface_plot3D(X, Y, w_tot)
plt.show()

Here is updated code with comparison:

import numpy as np

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt


def surface_plot3D(X,Y,Z):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                           linewidth=0, antialiased=False)
    ax.set_zlim(-3., 3.)

    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

    fig.colorbar(surf, shrink = 0.5, aspect = 5)

def funct(x,y):
    return np.exp(np.sin( np.sqrt(x**2 + y**2) )) 

#derivative wrt x
def dfunct_x(x,y):
    return np.exp(np.sin( np.sqrt(x**2 + y**2) )) * np.cos( np.sqrt(x**2 + y**2) ) * y/np.sqrt(x**2 + y**2)

#derivative wrt to y
def dfunct_y(x,y):
    return np.exp(np.sin( np.sqrt(x**2 + y**2) )) * np.cos( np.sqrt(x**2 + y**2) ) * x/np.sqrt(x**2 + y**2)  


N = 2**4
xmin=0
xmax=2.0*np.pi
step = (xmax-xmin)/(N)
xdata = np.linspace(step, xmax, N)

ydata = np.copy(xdata)
X,Y = np.meshgrid(xdata,ydata)
Z = funct(X,Y)
#calling the derivative function wrt x
df_x = dfunct_x(X,Y)
#calling the derivative function wrt y
df_y = dfunct_y(X,Y)
mag_df = np.sqrt(df_x**2 + df_y**2 )

surface_plot3D(X,Y,Z)
surface_plot3D(X,Y,mag_df)
plt.show()

vhat = np.fft.fft2(Z)

what = 1j * np.fft.fftfreq(N,1./N)
whax, whay = np.meshgrid(what, what)

wha_x = whax * vhat  
w_x = np.real(np.fft.ifft2(wha_x)) 
wha_y = whay * vhat  
w_y = np.real(np.fft.ifft2(wha_y)) 
w_tot = np.sqrt( w_x**2 + w_y**2)
surface_plot3D(X, Y, w_tot)
plt.show()
snano
  • 23
  • 7
  • How is the result wrong? – Jonas Adler Jul 14 '17 at 08:11
  • If we take derivative of the function wrt x and y, we will get w_x and w_y; then comparing with calculated w_tot using the above code does not agree. – snano Jul 14 '17 at 09:12
  • My point is perhaps, in what sense does it not agree? Is it of by a scale factor? By a very small noise? etc. Give us a figure or something. – Jonas Adler Jul 14 '17 at 09:13
  • I mean w_tot = sqrt( w_x^2 + w_y^2) calculated using direct derivative does not agree with that calculated from above script. – snano Jul 14 '17 at 09:23
  • I have added the comparison part now in the script; needs your suggestion. – snano Jul 14 '17 at 09:42

0 Answers0