草庐IT

Python有限差分函数？

coder 2023-05-23 原文

我一直在 Numpy/Scipy 中寻找包含有限差分函数的模块。但是，我发现的最接近的东西是 numpy.gradient()，它适用于二阶精度的一阶有限差分，但如果您想要高阶导数，则效果不佳或更准确的方法。我什至还没有为这类事情找到很多特定的模块。大多数人似乎都在根据需要做“自己动手”的事情。所以我的问题是，是否有人知道任何专门用于高阶(精度和导数)有限差分方法的模块(Numpy/Scipy 的一部分或第三方模块)。我有自己的代码正在处理，但它目前有点慢，如果有可用的东西，我不会尝试对其进行优化。

请注意，我说的是有限差分，而不是导数。我见过 scipy.misc.derivative() 和 Numdifftools ，它采用我没有的分析函数的导数。

最佳答案

快速做到这一点的一种方法是使用高斯核的导数进行卷积。简单的情况是你的数组与 [-1, 1] 的卷积，它给出了简单的有限差分公式。除此之外， (f*g)'= f'*g = f*g' 其中 * 是卷积，所以你最终得到你的导数与一个普通的卷积高斯，所以这当然会使你的数据平滑一点，可以通过选择最小的合理内核来最小化。

import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt

#Data:
x = np.linspace(0,2*np.pi,100)
f = np.sin(x) + .02*(np.random.rand(100)-.5)

#Normalization:
dx = x[1] - x[0] # use np.diff(x) if x is not uniform
dxdx = dx**2

#First derivatives:
df = np.diff(f) / dx
cf = np.convolve(f, [1,-1]) / dx
gf = ndimage.gaussian_filter1d(f, sigma=1, order=1, mode='wrap') / dx

#Second derivatives:
ddf = np.diff(f, 2) / dxdx
ccf = np.convolve(f, [1, -2, 1]) / dxdx
ggf = ndimage.gaussian_filter1d(f, sigma=1, order=2, mode='wrap') / dxdx

#Plotting:
plt.figure()
plt.plot(x, f, 'k', lw=2, label='original')
plt.plot(x[:-1], df, 'r.', label='np.diff, 1')
plt.plot(x, cf[:-1], 'r--', label='np.convolve, [1,-1]')
plt.plot(x, gf, 'r', label='gaussian, 1')
plt.plot(x[:-2], ddf, 'g.', label='np.diff, 2')
plt.plot(x, ccf[:-2], 'g--', label='np.convolve, [1,-2,1]')
plt.plot(x, ggf, 'g', label='gaussian, 2')

既然你提到了 np.gradient 我假设你至少有 2d 数组，所以以下适用于:如果你是内置在 scipy.ndimage 包中想为ndarrays做。不过要小心，因为这当然不会给你完整的渐变，但我相信所有方向的产品。希望有更专业的人说出来。

这是一个例子:

from scipy import ndimage

x = np.linspace(0,2*np.pi,100)
sine = np.sin(x)

im = sine * sine[...,None]
d1 = ndimage.gaussian_filter(im, sigma=5, order=1, mode='wrap')
d2 = ndimage.gaussian_filter(im, sigma=5, order=2, mode='wrap')

plt.figure()

plt.subplot(131)
plt.imshow(im)
plt.title('original')

plt.subplot(132)
plt.imshow(d1)
plt.title('first derivative')

plt.subplot(133)
plt.imshow(d2)
plt.title('second derivative')

使用 gaussian_filter1d 允许您沿某个轴进行方向导数:

imx = im * x
d2_0 = ndimage.gaussian_filter1d(imx, axis=0, sigma=5, order=2, mode='wrap')
d2_1 = ndimage.gaussian_filter1d(imx, axis=1, sigma=5, order=2, mode='wrap')

plt.figure()
plt.subplot(131)
plt.imshow(imx)
plt.title('original')
plt.subplot(132)
plt.imshow(d2_0)
plt.title('derivative along axis 0')
plt.subplot(133)
plt.imshow(d2_1)
plt.title('along axis 1')

上面的第一组结果让我有点困惑(当曲率应该指向向下时，原始的峰值显示为二阶导数的峰值)。在没有进一步研究 2d 版本的工作原理的情况下，我只能真正推荐 1d 版本。如果您想要大小，只需执行以下操作:

d2_mag = np.sqrt(d2_0**2 + d2_1**2)

关于Python有限差分函数？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/18991408/

Python有限39codepltnumpyscipy

有关Python有限差分函数？的更多相关文章

  1. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声？ - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题？更新问题，以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点？我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等)，但我确实想创建一个输出文件。

  2. ruby - 在没有 sass 引擎的情况下使用 sass 颜色函数 - 2

    我想在一个没有Sass引擎的类中使用Sass颜色函数。我已经在项目中使用了sassgem，所以我认为搭载会像以下一样简单:classRectangleincludeSass::Script::FunctionsdefcolorSass::Script::Color.new([0x82,0x39,0x06])enddefrender#hamlengineexecutedwithcontextofself#sothatwithintemlateicouldcall#%stop{offset:'0%',stop:{color:lighten(color)}}endend更新:参见上面的#re

  3. ruby-on-rails - 在 ruby​​ 中使用 gsub 函数替换单词 - 2

    我正在尝试用ruby​​中的gsub函数替换字符串中的某些单词，但有时效果很好，在某些情况下会出现此错误？这种格式有什么问题吗NoMethodError(undefinedmethod`gsub!'fornil:NilClass):模型.rbclassTest"replacethisID1",WAY=>"replacethisID2andID3",DELTA=>"replacethisID4"}end另一个模型.rbclassCheck 最佳答案 啊，我找到了!gsub!是一个非常奇怪的方法。首先，它替换了字符串，所以它实际上修改了

  4. ruby - 在 Ruby 中有条件地定义函数 - 2

    我有一些代码在几个不同的位置之一运行:作为具有调试输出的命令行工具，作为不接受任何输出的更大程序的一部分，以及在Rails环境中。有时我需要根据代码的位置对代码进行细微的更改，我意识到以下样式似乎可行:print"Testingnestedfunctionsdefined\n"CLI=trueifCLIdeftest_printprint"CommandLineVersion\n"endelsedeftest_printprint"ReleaseVersion\n"endendtest_print()这导致:TestingnestedfunctionsdefinedCommandLin

  5. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意，但是谷歌搜索||=并不是很有帮助；)Python中是否有与Ruby和Perl中的||=语句等效的语句？例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外，类似这样的东西的通用术语是什么？条件分配是我的第一个猜测，但Wikipediapage跟我想的不太一样。

  6. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi？ - 2

    什么是ruby​​的rack或python的Java的wsgi？还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

  7. ruby - 在 Ruby 中按名称传递函数 - 2

    如何在Ruby中按名称传递函数？(我使用Ruby才几个小时，所以我还在想办法。)nums=[1,2,3,4]#Thisworks,butismoreverbosethanI'dlikenums.eachdo|i|putsiend#InJS,Icouldjustdosomethinglike:#nums.forEach(console.log)#InF#,itwouldbesomethinglike:#List.iternums(printf"%A")#InRuby,IwishIcoulddosomethinglike:nums.eachputs在Ruby中能不能做到类似的简洁？我可以只

  8. 华为OD机试用Python实现 -【明明的随机数】 2023Q1A - 2

    华为OD机试题本篇题目：明明的随机数题目输入描述输出描述：示例1输入输出说明代码编写思路最近更新的博客华为od2023|什么是华为od，od薪资待遇，od机试题清单华为OD机试真题大全，用Python解华为机试题|机试宝典【华为OD机试】全流程解析+经验分享,题型分享,防作弊指南华为o

  9. python - 如何读取 MIDI 文件、更改其乐器并将其写回？ - 2

    我想解析一个已经存在的.mid文件，改变它的乐器，例如从“acousticgrandpiano”到“violin”，然后将它保存回去或作为另一个.mid文件。根据我在文档中看到的内容，该乐器通过program_change或patch_change指令进行了更改，但我找不到任何在已经存在的MIDI文件中执行此操作的库.他们似乎都只支持从头开始创建的MIDI文件。 最佳答案 MIDIpackage会为您完成此操作，但具体方法取决于midi文件的原始内容。一个MIDI文件由一个或多个音轨组成，每个音轨是十六个channel中任何一个上的

  10. 「Python｜Selenium｜场景案例」如何定位iframe中的元素？ - 2

    本文主要介绍在使用Selenium进行自动化测试或者任务时，对于使用了iframe的页面，如何定位iframe中的元素文章目录场景描述解决方案具体代码场景描述当我们在使用Selenium进行自动化测试的时候，可能会遇到一些界面或者窗体是使用HTML的iframe标签进行承载的。对于iframe中的标签，如果直接查找是无法找到的，会抛出没有找到元素的异常。比如近在咫尺的例子就是，CSDN的登录窗体就是使用的iframe，大家可以尝试通过F12开发者模式查看到的tag_name,class_name,id或者xpath来定位中的页面元素，会抛出NoSuchElementException异常。解决

随机推荐