python函数求导_python – 来自scipy.special的fadeeva函数的二阶导数

我想计算Fadeeva函数special.wofz的二阶导数. Fadeeva函数与误差函数密切相关.因此,如果有人更熟悉erf,那么答案是值得赞赏的.

这是找到wofz的二阶导数的代码：

import numpy as np

import matplotlib.pyplot as plt

from scipy.special import wofz

def Z(x):

return wofz(x)

## first derivative of wofz (analytically)

def Zp(x):

return -2/1j/np.pi**0.5 - 2*x*Z(x)

##second derivative (analytically)

def Zpp(x):

return (Z(x)+x*Zp(x))*x

x = np.float64(np.linspace(1e4,14e4,1000))

plt.plot(x, Zpp(x).imag,"-")

Zpp_num=np.diff(Zp(x))/np.diff(x) ##calc numerically the second derivative

plt.plot(x[:-1],Zpp_num.imag)

代码生成下一个数字：

显然,分析计算存在严重问题.我一直在使用的公式是正确的.它必须是一些数字问题.

问：有人能告诉我这种行为的原因是什么吗？是否由于wofz功能的精确性？有谁知道计算wofz的算法？可以产生可靠结果的论点有多大？我找不到任何关于它的信息.另外,我知道我可以使用wofz的渐近逼近来找到二阶导数但是如果可能的话我想使用scipy.

最佳答案 按照@Andras Deak的回答,您可以分析地找出高x扩展,然后使用一些简单的平滑在它和scipy函数之间进行插值.实际上有两个术语在高x扩展中取消,所以你必须要小心一点.

这是我得到的答案：

import numpy as np

import matplotlib.pyplot as plt

from scipy.special import wofz

def Z(x):

return wofz(x)

## first derivative of wofz (analytically)

def Zp(x):

return -2/1j/np.pi**0.5 - 2*x*Z(x)

def dawsn_expansion(x):

# Accurate to order x^-9, or, relative to the first term x^-8

# So when x > 100, this will be as accurate as you can get with

# double floating point precision.

y = 0.5 * x**-2

return 1/(2*x) * (1 + y * (1 + 3*y * (1 + 5*y * (1 + 7*y))))

def dawsn_expansion_drop_first(x):

y = 0.5 * x**-2

return 1/(2*x) * (0 + y * (1 + 3*y * (1 + 5*y * (1 + 7*y))))

def dawsn_expansion_drop_first_two(x):

y = 0.5 * x**-2

return 1/(2*x) * (0 + y * (0 + 3*y * (1 + 5*y * (1 + 7*y))))

def blend(x, a, b):

# Smoothly blend x from 0 at a to 1 at b

y = (x - a) / (b - a)

y *= (y > 0)

y = y * (y <= 1) + 1 * (y > 1)

return y*y * (3 - 2*y)

def g(x):

"""Calculate `x + (1-2x^2) D(x)`, where D(x) is the dawson function"""

# For x < 50, use dawsn from scipy

# For x > 100, use dawsn expansion

b = blend(x, 50, 100)

y1 = x + (1 - 2*x**2) * special.dawsn(x)

y2 = dawsn_expansion_drop_first(x) - dawsn_expansion_drop_first_two(x) * 2*x**2

return b*y2 + (1-b)*y1

def Zpp(x):

# only return the imaginary component

return -4j/np.pi**0.5 * g(x)

x = np.logspace(0, 5, 2000)

dx = 1e-3

plt.plot(x, (Zp(x+dx) - Zp(x-dx)).imag/(2*dx))

plt.plot(x, Zpp(x).imag)

ax = plt.gca()

ax.set_xscale('log')

ax.set_yscale('log')

产生以下图：

蓝线是数值导数,绿线是使用扩展的导数.后者实际上在大x时具有更好的行为.