quadrature

数值积分

trapz - 使用梯形法则积分采样值

状态

实验性

描述

返回表示函数离散样本的数组 y 的梯形法则积分。积分是假设等距横坐标间距为 dx 或任意横坐标 x 计算的。

语法

result = trapz (y, x)

result = trapz (y, dx)

参数

y: 应为类型为 real 的一维数组。

x: 应为类型为 real 的一维数组,其类型和大小与 y 相同。

dx: 应为类型为 real 的标量,其类型与 y 相同。

返回值

结果是类型为 real 的标量,其类型与 y 相同。

如果 y 的大小为零或一,则结果为零。

示例

program example_trapz
  use stdlib_quadrature, only: trapz
  implicit none
  real, parameter :: x(5) = [0., 1., 2., 3., 4.]
  real :: y(5) = x**2
  print *, trapz(y, x)
! 22.0
  print *, trapz(y, 0.5)
! 11.0
end program example_trapz

trapz_weights - 给定横坐标的梯形法则权重

状态

实验性

描述

给定横坐标数组 x,计算权重数组 w,使得如果 y 代表在 x 处制表的函数值,则 sum(w*y) 产生对积分的梯形法则近似值。

语法

result = trapz_weights (x)

参数

x: 应为类型为 real 的一维数组。

返回值

结果是一个 real 数组,其大小和类型与 x 相同。

如果 x 的大小为一,则结果的唯一元素为零。

示例

program example_trapz_weights
  use stdlib_quadrature, only: trapz_weights
  implicit none
  real, parameter :: x(5) = [0., 1., 2., 3., 4.]
  real :: y(5) = x**2
  real :: w(5)
  w = trapz_weights(x)
  print *, sum(w*y)
! 22.0
end program example_trapz_weights

simps - 使用辛普森法则积分采样值

状态

实验性

描述

返回表示函数离散样本的数组 y 的辛普森法则积分。积分是假设等距横坐标间距为 dx 或任意横坐标 x 计算的。

对于奇数长度的数组,使用辛普森普通(“1/3”)法则。对于偶数长度的数组,还将根据 even 的值以某种方式使用辛普森 3/8 法则。如果 even 为负(正),则在数组的开头(结尾)使用 3/8 法则。如果 even 为零或不存在,则结果与首先在数组开头使用 3/8 法则,然后在数组结尾使用 3/8 法则,并将这两个结果取平均值相同。

语法

result = simps (y, x [, even])

result = simps (y, dx [, even])

参数

y: 应为类型为 real 的一维数组。

x: 应为类型为 real 的一维数组,其类型和大小与 y 相同。

dx: 应为类型为 real 的标量,其类型与 y 相同。

even: (可选)应为默认类型的 integer

返回值

结果是类型为 real 的标量,其类型与 y 相同。

如果 y 的大小为零或一,则结果为零。

如果 y 的大小为二,则结果与调用 trapz 相同。

示例

program example_simps
  use stdlib_quadrature, only: simps
  implicit none
  real, parameter :: x(5) = [0., 1., 2., 3., 4.]
  real :: y(5) = 3.*x**2
  print *, simps(y, x)
! 64.0
  print *, simps(y, 0.5)
! 32.0
end program example_simps

simps_weights - 给定横坐标的辛普森法则权重

状态

实验性

描述

给定横坐标数组 x,计算权重数组 w,使得如果 y 代表在 x 处制表的函数值,则 sum(w*y) 产生对积分的辛普森法则近似值。

对于奇数长度的数组,使用辛普森普通(“1/3”)法则。对于偶数长度的数组,还将根据 even 的值以某种方式使用辛普森 3/8 法则。如果 even 为负(正),则在数组的开头(结尾)使用 3/8 法则,并在其他地方使用 1/3 法则。如果 even 为零或不存在,则结果与首先在数组开头使用 3/8 法则,然后在数组结尾使用 3/8 法则,最后将这两个结果取平均值相同。

语法

result = simps_weights (x [, even])

参数

x: 应为类型为 real 的一维数组。

even: (可选)应为默认类型的 integer

返回值

结果是一个 real 数组,其大小和类型与 x 相同。

如果 x 的大小为一,则结果的唯一元素为零。

如果 x 的大小为二,则结果与调用 trapz_weights 相同。

示例

program example_simps_weights
  use stdlib_quadrature, only: simps_weights
  implicit none
  real, parameter :: x(5) = [0., 1., 2., 3., 4.]
  real :: y(5) = 3.*x**2
  real :: w(5)
  w = simps_weights(x)
  print *, sum(w*y)
! 64.0
end program example_simps_weights

gauss_legendre - 高斯-勒让德求积法(又称高斯求积法)节点和权重

状态

实验性

描述

计算高斯-勒让德求积法(也称为高斯求积法)节点和权重,适用于任何 N(节点数)。使用节点 x 和权重 w,您可以计算某个函数 f 的积分,如下所示:integral = sum(f(x) * w)

只支持双精度 - 如果需要较低的精度,则必须自行执行适当的转换。精度已通过将计算结果与已知精确到机器精度的表格值进行比较来验证,直到 N=64(与这些值的最大差异为 2 epsilon)。

语法

subroutine gauss_legendre (x, w[, interval])

参数

x: 应为类型为 real(real64) 的一维数组。它是一个 *输出* 参数,表示求积节点。

w: 应为类型为 real(real64) 的一维数组,其维度与 x 相同。它是一个 *输出* 参数,表示求积权重。

interval: (可选)应为类型为 real(real64) 的双元素数组。如果存在,则从 interval(1)interval(2) 计算节点和权重。如果未指定,则默认积分范围为 -1 到 1。

示例

program example_gauss_legendre
  use iso_fortran_env, dp => real64
  use stdlib_quadrature, only: gauss_legendre
  implicit none

  integer, parameter :: N = 6
  real(dp), dimension(N) :: x, w
  call gauss_legendre(x, w)
  print *, "integral of x**2 from -1 to 1 is", sum(x**2*w)
end program example_gauss_legendre

gauss_legendre_lobatto - 高斯-勒让德-罗巴托求积法节点和权重

状态

实验性

描述

计算高斯-勒让德-罗巴托求积法节点和权重,适用于任何 N(节点数)。使用节点 x 和权重 w,您可以计算某个函数 f 的积分,如下所示:integral = sum(f(x) * w)

只支持双精度 - 如果需要较低的精度,则必须自行执行适当的转换。精度已通过将计算结果与已知精确到机器精度的表格值进行比较来验证,直到 N=64(与这些值的最大差异为 2 epsilon)。

语法

subroutine gauss_legendre_lobatto (x, w[, interval])

参数

x: 应为类型为 real(real64) 的一维数组。它是一个 *输出* 参数,表示求积节点。

w: 应为类型为 real(real64) 的一维数组,其维度与 x 相同。它是一个 *输出* 参数,表示求积权重。

interval: (可选)应为类型为 real(real64) 的双元素数组。如果存在,则从 interval(1)interval(2) 计算节点和权重。如果未指定,则默认积分范围为 -1 到 1。

示例

program example_gauss_legendre_lobatto
  use iso_fortran_env, dp => real64
  use stdlib_quadrature, only: gauss_legendre_lobatto
  implicit none

  integer, parameter :: N = 6
  real(dp), dimension(N) :: x, w
  call gauss_legendre_lobatto(x, w)
  print *, "integral of x**2 from -1 to 1 is", sum(x**2*w)
end program example_gauss_legendre_lobatto