这是找到wofz的二阶导数的代码:
import numpy as npimport matplotlib.pyplot as pltfrom scipy.special import wofzdef 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))*xx = 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 derivativeplt.plot(x[:-1],Zpp_num.imag)
代码生成下一个数字:
显然,分析计算存在严重问题.我一直在使用的公式是正确的.它必须是一些数字问题.
问:有人能告诉我这种行为的原因是什么吗?是否由于wofz功能的精确性?有谁知道计算wofz的算法?可以产生可靠结果的论点有多大?我找不到任何关于它的信息.另外,我知道我可以使用wofz的渐近逼近来找到二阶导数但是如果可能的话我想使用scipy.
解决方法 按照@Andras Deak的回答,您可以分析地找出高x扩展,然后使用一些简单的平滑在它和scipy函数之间进行插值.实际上有两个术语在高x扩展中取消,所以你必须要小心一点.这是我得到的答案:
import numpy as npimport matplotlib.pyplot as pltfrom scipy.special import wofzdef 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)*y1def Zpp(x): # only return the imaginary component return -4j/np.pi**0.5 * g(x)x = np.logspace(0,5,2000)dx = 1e-3plt.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时具有更好的行为.
总结以上是内存溢出为你收集整理的python – 来自scipy.special的fadeeva函数的二阶导数全部内容,希望文章能够帮你解决python – 来自scipy.special的fadeeva函数的二阶导数所遇到的程序开发问题。
如果觉得内存溢出网站内容还不错,欢迎将内存溢出网站推荐给程序员好友。
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)