CODEWITHSUNDEEP
pythonscipyinterpolation
rland
rland

Reputation: 123

Calculate Line Integral over Bivariate Spline

I can construct a spline of some 2-D data and query area integrals quite easily with scipy.interpolate.RectBivariateSpline. But there is not a function for easily/quickly calculating line integrals over such splines.

from scipy.interpolate import RectBivariateSpline
import numpy as np

# example data
X, Y = np.meshgrid(np.arange(0, 100), np.arange(0, 100))
Z = np.sin(X / 10) + np.cos(Y / 10)

# construct the spline
rbs = RectBivariateSpline(X[0, :], Y[:, 0], Z)

# get integral
area_integral = rbs.integral(0, 10, 30, 50)

print(area_integral)

>>>40.9287528271697

What I am looking for is something like rbs.line_integral(0,10, 30,50) which will give me the integral over the line from [0,30] -> [10,50]; not the area of the rectangle formed by those two points.

Is there such a function in fitpack? I can see that under the hood, scipy makes a call to dfitpack.dblint() but I'm not really sure where to go from there.

I should add also: I don't want to make a series of points on a line, query those points with e.g. rbs.ev(x, y), and then sum the result. That is very slow and probably introduces numerical integration errors.

Edit

What RectBivariateSpline.integral() gives me is:

\int_A f(x, y) dxdy

Where f(x, y) is the function approximated by the spline, and A is the area of the rectangle between points p1=(x0, y0) p2=(x1, y1). So, dxdy is a differential area, A is an area, and the evaluated integral is a volume.

What I am looking for is:

\int_C f(x, y) ds

Where f(x, y) is the function approximated by the spline, and C is a straight line between p1=(x0,y0) and p2=(x1, y1). So ds is a differential length, S is a length, and the evaluated integral is an area.

Upvotes: 2

Views: 397

Answers (1)

user1196549
user1196549

Reputation:

A bivariate spline surface is piecewise polynomial. In the case of a bicubic expression, for instance, by substituting the parametric equations X = X0 + t.δX, Y = Y0 + t.δY, you obtain a sextic polynomial in t. The parameters X0, Y0, δX and δY as well as the integration range are found in every grid cell by tracing the line segment and intersecting it with the grid. Then if scipy is kind enough to give you the coefficients of the bivariate polynomials crossed, you can obtain those of the polynomials in t.

If the expression of ds is √dx²+dy², this amounts to √δX²+δY².dt and the antiderivative is elementary. So it "suffices" to enumerate the polynomials in the cells crossed and compute all required coefficients.

If the expression is √dx²+dy²+dz², you are stuck because by the chain rule dz = (δX.P(t) + δY.Q(t)) dt, where P, Q are polynomials, and the integrals do not have a closed-form expression.

Note. Due to the complexity of the polynomials expansions, it is not unthinkable that a straight Simpson integral evaluation could be faster and as accurate, for a reasonable step size.

Upvotes: 1

Related Questions