python中插值函数的加速积分

2024-04-26 17:54:28 发布

您现在位置:Python中文网/ 问答频道 /正文

我有一个函数f(x),它没有解析形式,但是来自一组数据点的插值(依赖于许多其他参数),我需要对这些数据点进行一百万次的积分。问题是被积函数在两端发散,因此我似乎被quad积分困住了,这使得它非常慢。f(x)的形状示例如下。所以我必须插值f(x),然后积分:

import numpy as np
from scipy.integrate import quad
from scipy.interpolate import UnivariateSpline

x = np.array([6.3096e-04, 8.8499e-04, 1.2413e-03, 1.7411e-03, 2.4421e-03, 3.4253e-03,
 4.8043e-03, 6.7386e-03, 9.4517e-03, 1.3257e-02, 1.8595e-02, 2.6081e-02,
 3.6582e-02, 5.1310e-02, 7.1969e-02, 1.0094e-01, 1.4159e-01, 1.9859e-01,
 2.7855e-01, 3.9069e-01, 5.4799e-01, 7.6862e-01, 1.0781e+00, 1.5121e+00,
 2.1210e+00, 2.9749e+00, 4.1726e+00, 5.8526e+00, 8.2089e+00, 1.1514e+01,
 1.6150e+01, 2.2652e+01, 3.1772e+01, 4.4564e+01, 6.2506e+01, 8.7671e+01,
 1.2297e+02, 1.7248e+02, 2.4192e+02, 3.3932e+02, 4.7594e+02, 6.6756e+02,
 9.3633e+02, 1.3133e+03, 1.8421e+03, 2.5837e+03, 3.6240e+03, 5.0830e+03,
 7.1295e+03, 1.0000e+04])
f = np.array([ 3.3914e+06,  2.3765e+06,  1.6964e+06,  1.1964e+06,  8.3927e+05,
  5.9232e+05,  4.1550e+05,  2.8922e+05,  1.9973e+05,  1.3602e+05,
  9.1304e+04,  5.9997e+04,  3.8402e+04,  2.3764e+04,  1.4142e+04,
  8.0330e+03,  4.3891e+03,  2.2771e+03,  1.1380e+03,  5.5416e+02,
  2.4666e+02,  7.1272e+01,  3.6225e+01,  2.6173e+01,  1.8561e+01,
  1.2739e+01,  8.4143e+00,  5.3091e+00,  3.1677e+00,  1.7646e+00,
  9.0203e-01,  4.1415e-01,  1.6602e-01,  5.5657e-02,  1.4487e-02,
  2.9594e-03,  1.5416e-03, -9.8210e-05,  1.5371e-04,  2.2051e-04,
  1.5187e-04,  9.9361e-05,  6.0271e-05,  3.0741e-05,  1.3920e-05,
  6.0695e-06,  2.7415e-06,  1.3585e-06,  7.5640e-07,  4.6422e-07])

rx = np.array([3.1623e-03, 4.5099e-03, 6.4318e-03, 9.1728e-03, 1.3082e-02, 1.8657e-02,
 2.6607e-02, 3.7946e-02, 5.4117e-02, 7.7179e-02, 1.1007e-01, 1.5698e-01,
 2.2387e-01, 3.1928e-01, 4.5534e-01, 6.4938e-01, 9.2612e-01, 1.3208e+00,
 1.8836e+00, 2.6864e+00, 3.8312e+00, 5.4639e+00, 7.7923e+00, 1.1113e+01,
 1.5849e+01])

finterp = UnivariateSpline(x, f, s=0, ext=0)
integrand = lambda x, x0: finterp(x0/x) / (x**2 * (1-x**2)**0.5)
out = np.array([quad(integrand, 0, 1, args=(x0,), limit=100)
                for x0 in rx])

scipy.integrate中的其他函数都没有给出合理的值,特别是对于较大的rx值。上面的数组f每次都是不同的(这取决于许多参数),所以存储f不是很实用(或者至少不是以我能想到的方式)。此外,上面的代码发生在一个已经与multiprocessing合并的更大的程序中,因此我无法并行化上面的代码。你知道吗

任何帮助都将不胜感激!你知道吗


Tags: 数据函数fromimport参数npscipyrx