数值积分

2024-12-07

数值积分（精选8篇）

数值积分篇1

一、引言

实际问题当中常常需要计算积分, 依据高等数学中的微积分基本定理, 对于积分I (f) =∫ $_{a}^{b}$ f (x) dx只要找到被积函数f (x) 的原函数F (x) , 便可利用牛顿—莱布尼兹公式∫ $_{a}^{b}$ f (x) dx=F (b) -F (a) 进行计算.然而, 在现代科学技术中使用这种求积方法往往有困难, 因为: (1) 当f (x) 是由测量或数值计算给出的一张数据表时, 无法得到函数f (x) 的原函数F (x) 的具体形式, 从而牛顿-莱布尼兹公式就不能直接运用; (2) 大量的被积函数如 $\frac{s i n x}{x}, \sin x^{2}, \sqrt{1 + x^{3}}, e^{- x^{2}}$ 等等, 很难找到用初等函数表示的原函数F (x) , 固然有些函数的原函数F (x) 可以找到, 但在应用牛顿—莱布尼兹公式时, 涉及大量的函数值计算, 还不如应用数值积分的方法来得方便, 既节省工作量, 又满足精度要求, 因此很有必要研究积分的数值计算问题.

在积分的数值计算中, 并不一定所有的积分都是在定区间上积分, 常常遇到在无穷区间上的积分问题, 如I (f) =∫ $_{a}^{+ \infty}$ f (x) dx, I (f) =∫ $_{- \infty}^{a}$ f (x) dx, I (f) =∫ $_{- \infty}^{+ \infty}$ f (x) dx等类似问题.本文便从这一问题出发, 将无穷积分转换为定积分, 对于定积分问题, 我们便可以用相关的复合求积公式进行求解.

二、数值积分

要计算定积分I (f) =∫ $_{a}^{b}$ f (x) dx (y=f (x) 为已知的可积函数) , 我们可以由定积分的定义知, 定积分∫ $_{a}^{b}$ f (x) dx是和的极限, 数值积分就是将定积分∫ $_{a}^{b}$ f (x) dx的计算用有限和的形式近似, 由于多项式的积分非常容易计算, 和建立数值微分公式一样, 我们可用函数f (x) 的拉格朗日插值多项式Ln (x) 的积分Q[f]近似I[f], 即

I[f]=∫ $_{a}^{b}$ f (x) dx≈Q[f]= $\sum_{i = 0}^{n}$ Aif (xi) ,

其中Ai由式子

$\int_{a}^{b} \prod_{\begin{array}{l} j \neq i \\ j = 0 \end{array}}^{n}$

$\frac{x - x_{j}}{x_{i} - x_{j}} d x$ 确定, 称为求积系数, 它仅与节点有关而与被积函数无关.

数值积分的特点是直接用一些离散节点上的函数值f (xi) 的线性组合计算定积分的近似值, 从而将定积分的计算归结为函数值的计算, 这就避开了牛顿—莱布尼兹公式中需要寻求原函数的困难, 并为用计算机求积分提供了可行性.

在实际计算中, 由于高阶求积公式的计算过程可能出现数值不稳定性以及低阶求积公式不能满足精度要求, 通常使用复合求积公式, 下面便是复合抛物线求积公式的介绍:

将积分区间[a, b]分成2n个等长的小区间[xi, xi+1], 然后在每两个相邻小区间[xi, xi+2]上应用抛物线求积公式, 最后相加即得到复合抛物线求积公式.

令 $h = \frac{b - a}{2 n}, x_{i} = a + i h (i = 0, 1, 2, \dots, 2 n) .$

在区间[x2i, x2i+2] (i=0, 1, 2, …, n-1) 上利用抛物线求积公式然后叠加得

I[f]=∫ $_{x_{0}}^{x_{2}}$ f (x) dx+∫ $_{x_{2}}^{x_{4}}$ f (x) dx+…+∫ $_{x_{2 n - 2}}^{x_{2 n}}$ f (x) dx≈ $\frac{h}{3} [f (x_{0}) + 4 f (x_{1}) + f (x_{2})] + \frac{h}{3} [f (x_{2}) + 4 f (x_{3}) + f (x_{4})] + \dots + \frac{h}{3} [f (x_{2 n - 2}) + 4 f (x_{2 n - 1}) + f (x_{2 n})] = \frac{h}{3}$ $f (a) + 4 \sum_{i = 1}^{n} f (x_{2 i - 1}) + 2 \sum_{i = 1}^{n - 1} f (x_{2 i}) + f (b)$ =S (h) .

求积公式S (h) 称为复合抛物线求积公式.

复合抛物线求积公式不但容易编排程序上机计算, 而且精度也比较高, 是一个较好的数值积分方法, 应用较广泛.

三、瑕积分的计算

当被积函数在积分区间上有瑕点时, 积分I (f) =∫baf (x) dx为瑕积分, 为方便起见, 对于瑕积分, 仅考虑x=a有瑕点的瑕积分.

事实上, 若函数在积分区间上有多个瑕点, 则可将其化成多个仅在端点上有一个瑕点的积分.而此时, 又若x=b为函数f (x) 的瑕点, 则可作变量t=-x, 使得被积函数的瑕点落在左端点上.

考虑瑕积分:I (f) =∫ $_{a}^{b}$ f (x) dx, 其中 $f (x) = \frac{g (x)}{(x - a)^{p}}$ , 而0<p<1, g (x) 在区间[a, b]上连续.

下面主要介绍瑕积分计算中的分段积分法:

分段积分法的思想是将积分分成两部分, 一部分为含有瑕点的积分, 此积分可以解析求出来;另一部分为具有某种光滑程度的函数的积分, 此积分可用前面讲过的求积方法计算求解.设在瑕积分I (f) =∫ $_{a}^{b}$ f (x) dx中g (x) ∈C5[a, b], 则由泰勒展开式

$g (x) \approx Ρ_{4} (x) = g (a) + g^{'} (a) (x - a) + \frac{g^{″} (x)}{2!} (x - a)^{2} + \frac{g^{‴} (x)}{3!} (x - a)^{3} + \frac{g^{(4)} (x)}{4!} (x - a)^{4} .$

将原积分分成两部分:

I (f) =∫ $_{a}^{b}$ f (x) dx= $\int_{a}^{b} \frac{Ρ_{4} (x)}{(x - a)^{p}} d x + \int_{a}^{b} \frac{g (x) - Ρ_{4} (x)}{(x - a)^{p}} d x .$

其中一部分为瑕积分

$\int_{a}^{b} \frac{Ρ_{4} (x)}{(x - a)^{p}} d x = \int_{a}^{b} \sum_{k = 0}^{4} \frac{g^{(k)} (a)}{k!} (x - a)^{k - p} d x = \sum_{k = 0}^{4} \frac{g^{(k)} (a)}{k!} \int_{a}^{b} (x - a)^{k - p} d x = \sum_{k = 0}^{4} \frac{g^{(k)} (a)}{k! (k - p + 1)} (x - a)^{k - p + 1} |_{a}^{b} = \sum_{k = 0}^{4} \frac{g^{(k)} (a)}{k! (k - p + 1)} (b - a)^{k - p + 1}$ .

即得到解析解.

下面分析第二部分∫ $_{a}^{b}$ $\frac{g (x) - Ρ_{4} (x)}{(x - a)^{p}} d x$ 的计算, 定义

$G (x) = {\begin{cases} \frac{g (x) - Ρ_{4} (x)}{(x - a)^{p}}, a < x \leq b, \\ 0, x = a . \end{cases}$

由于0<p<1以及P $_{4}^{(k)}$ (a) =g (k) (a) , k=0, 1, 2, 3, 4知函数G (x) ∈C4[a, b], 因此可用复合抛物线求积公式计算积分 $\int_{a}^{b} \frac{g (x) - Ρ_{4} (x)}{(x - a)^{p}} d x = \int_{a}^{b}$ G (x) dx.

四、无穷积分的转换及数值计算

1.无穷积分的转换

对于无穷区间上的可积函数f (x) , 有关它的无穷积分主要有以下5种类型 (a>0) : (1) I (f) =∫ $_{a}^{+ \infty}$ f (x) dx, (2) I (f) =∫a-∞f (x) dx, (3) I (f) =∫ $_{- \infty}^{- a}$ f (x) dx, (4) I (f) =∫ $_{- a}^{+ \infty}$ f (x) dx, (5) I (f) =∫ $_{- \infty}^{+ \infty}$ f (x) dx, 一般情况下, 可对其进行转化, 使其变成有限区间上的定积分, 从而便于计算.事实上, a<0时积分也可以转化为以上类型中的积分.

类似第一种类型的无穷积分I (f) =∫+∞af (x) dx (a>0) , 可令 $x = \frac{1}{t}$ , 此时 $d x = - \frac{1}{t^{2}} d t$ , 因此,

I (f) = $\int_{a}^{+ \infty} f (x) d x = \int_{\frac{1}{a}}^{0} (- \frac{1}{t^{2}}) f (\frac{1}{t}) d t = \int_{0}^{\frac{1}{a}} t^{- 2} f (\frac{1}{t}) d t$ ,

从而化为我们可以处理的定积分.

类似第二种类型的无穷积分I (f) =∫ $_{- \infty}^{a}$ f (x) dx (a>0) , 我们可以对其进行分段积分, 有:

$Ι (f) = \int_{- \infty}^{a} f (x) d x = \int_{- \infty}^{- a} f (x) d x + \int_{- a}^{0} f (x) d x + \int_{0}^{a} f (x) d x = - \int_{+ \infty}^{a} f (- t) d t + (- 1) \int_{a}^{0} f (- t) d t + \int_{0}^{a} f (t) d t = \int_{0}^{\frac{1}{a}} f (t) d t + \int_{0}^{a} [f (- t) + f (t)] d t .$

同理可得:

$\begin{array}{l} Ι (f) = \int_{- \infty}^{- a} f (x) d x = \int_{a}^{+ \infty} f (- x) d x = \int_{0}^{\frac{1}{a}} t^{- 2} f (- \frac{1}{t}) d t, \\ Ι (f) = \int_{- a}^{+ \infty} f (x) d x = - \int_{a}^{- \infty} f (- t) d t = \int_{- \infty}^{a} f (- t) d t = \int_{0}^{\frac{1}{a}} t^{- 2} f (\frac{1}{t}) d t + \int_{0}^{a} [f (- t) + f (t)] d t, \\ Ι (f) = \int_{- \infty}^{+ \infty} f (x) d x = \int_{- \infty}^{a} f (x) d x + \int_{a}^{+ \infty} f (x) d x = \int_{0}^{\frac{1}{a}} t^{- 2} f (- \frac{1}{t}) d t + \int_{0}^{a} [f (- t) + f (t)] d t + \int_{0}^{\frac{1}{a}} t^{- 2} f (\frac{1}{t}) d t = \int_{0}^{\frac{1}{a}} t^{- 2} [f (- \frac{1}{t}) + f (\frac{1}{t})] d t + \int_{0}^{a} [f (- t) + f (t)] d t . \end{array}$

2.数值计算

例计算积分∫ $_{1}^{+ \infty}$ $\frac{\sin (\frac{1}{x})}{\sqrt[3]{x^{5}}} d x .$

解观察积分形式, 可知这种积分属于第一类无穷积分, 从而可对其进行变形 (这里a=1) :

I (f) =∫ $_{1}^{+ \infty}$ $\frac{\sin (\frac{1}{x})}{\sqrt[3]{x^{5}}} d x =$ $\int_{0}^{1} t^{- 2} \frac{s i n t}{\sqrt[3]{(\frac{1}{t})^{5}}} d t = \int_{0}^{1} t^{- \frac{1}{3}} \sin t d t$ .

这样就将无穷积分化成了定积分, 但存在瑕点x=0, 由瑕积分的分段积分法可得:

$g (t) = \sin t \approx g^{'} (0) + g (0) (x - 0) + \frac{g^{″} (0)}{2!} (x - 0)^{2} + \frac{g^{‴} (0)}{3!} (x - 0)^{3} + \frac{g^{(4)} (0)}{4!} (x - 0)^{4} = \sin (0) + \sin ´ (0) x + \frac{s i n^{″} (0)}{2!} x^{2} + \frac{s i n^{‴} (0)}{3!} x^{3} + \frac{\sin^{(4)} (0)}{4!} x^{4} = t - \frac{1}{6} t^{3} = Ρ_{4} (t) .$

可令

$G (t) = {\begin{cases} \frac{s i n t - Ρ_{4} (t)}{t^{\frac{1}{3}}} = \frac{s i n t - t + \frac{1}{6} t^{3}}{t^{\frac{1}{3}}}, 0 < t \leq 1 ‚ \\ 0, t = 0 . \end{cases}$

此时有

I (f) = $\int_{1}^{+ \infty} \frac{\sin (\frac{1}{x})}{\sqrt[3]{x^{5}}} d x = \int_{0}^{1} t^{- \frac{1}{3}} (t - \frac{1}{6} t^{3}) d t + \int_{0}^{1} G (t) d t = \frac{3}{5} - \frac{1}{22} + \int_{0}^{1} G (t) d t = \frac{61}{110} + \int_{0}^{1} G (t) d t .$

用复合抛物线求积公式对∫ $_{0}^{1}$ G (t) dt进行计算, 再加上 $\frac{61}{110}$ 便可得结果:

积分上限:1

积分下限:0

取n=200, 步长 $h = \frac{b - a}{2 n} = \frac{1 - 0}{400} = 0.0025$

算法: (略)

运行结果为:0.555910

所以, 用复合抛物线求积公式 (n=200) 的计算结果为:0.555910

五、结论

在进行积分的数值计算时, 无论碰到什么样的无穷积分, 我们都可以将无穷积分问题化为定积分问题, 考虑到被积函数在积分区间上是否存在瑕点采用分段积分法, 再用复合抛物线求积公式进行求解.从实验结果可以看出, 用这种求积方法算出的积分值与MATLAB计算出来的积分值非常接近甚至更精确.

数值积分篇2

结构拟动力实验数值积分方法的比较研究

在对结构拟动力实验显式数值积分方法的.理论推导的基础上,通过数值计算分别对中央差分法、显式Newmark法、双β参数法的稳定性和计算精度进行了比较研究.计算结果表明:中央差分法和显式Newmark法稳定条件为Ω<2,双β参数法是无条件稳定的;当Ω<0.8,3种方法都具有良好的计算精度.

作者：祝辉杨佑发 ZHU Hui YANG You-fa 作者单位：重庆大学土木工程学院,重庆,400045刊名：实验科学与技术英文刊名：EXPERIMENT SCIENCE AND TECHNOLOGY年，卷(期)：20097(3)分类号：P315.9 O241.84关键词：拟动力实验数值积分方法稳定性计算精度

Cotes数值积分公式的改进篇3

在数值积分计算中, Newton-Cotes求积公式被广泛应用。文献中对梯形求积公式和Simpson求积公式进行了改进[1—3]。Simpson3/8求积公式和Cotes求积公式也是两类重要的积分近似计算公式。

设函数f (x) 在区间[a, b]上有6阶导数, 则

$\begin{array}{l} \int_{a}^{b} f (x) d x = \frac{b - a}{90} [7 f (a) + 32 f (a + h) + \\ 12 f (\frac{a + b}{2}) + 32 f (b - h) + 7 f (b)] - \end{array}$

$\frac{8}{945} (\frac{b - a}{4})^{7} f^{(6)} (ξ)$ 。

其中h= (b-a) /4, ξ∈ (a, b) 。

若不考虑积分余项, 则有

$\begin{array}{l} \int_{a}^{b} f (x) d x \approx \frac{b - a}{90} [7 f (a) + 32 f (a + h) + \\ 12 f (\frac{a + b}{2}) + 32 f (b - h) + 7 f (b)] 。 \end{array}$

上式即为著名的Cotes求积公式。Cotes求积公式具有5次代数精确度。本文首先讨论Cotes数值求积公式余项表达式中中间点的渐进估计, 进一步利用中间点的渐进估计给出Cotes求积公式的改进。

1 Cotes求积公式中间点的渐进性

考虑Cotes求积公式积分余项中间点的渐进性质, 我们有如下结果。

定理1 设函数f (x) 在[a, b]上有7阶导数, 且f (7) (a) ≠0, 则有

$\lim_{b \to a^{+}} \frac{ξ - a}{b - a} = \frac{1}{2}$ (1)

证明由于所考虑的极限为b→a+的情形, 故题设条件和Taylor公式有

$\int_{a}^{b} f (x) d x = f (a) (b - a) + \dots + \frac{f^{(7)} (a)}{8!} (b - a)^{8} + o (b - a)^{8} f (\frac{a + b}{2}) = f (a) + f^{'} (a) (b - a) + \dots + \frac{f^{(7)} (a)}{7!} (\frac{b - a}{2})^{7} + o (b - a)^{7}$ 。

f (6) (ξ) =f (6) (a) +f (7) (a) (ξ-a) +o (ξ-a) 。

类似地, 可以将函数f (b) , f (a+h) , f (b-h) 在a点展开, 并将这些展开式代入Cotes求积公式积分余项表达式, 经过化简可得

$\frac{1}{2} f^{(7)} (a) (b - a) + o (b - a) = f^{(7)} (a) (ξ - a) + o (ξ - a)$ 。

由于b→a+, ξ∈ (a, b) , 所以ξ→a+, 又因为

f (7) (a) ≠0, 从而有 $\lim_{b \to a^{+}} \frac{ξ - a}{b - a} = \frac{1}{2}$ 。

定理1表明Cotes求积公式积分余项表达式中中间点的渐进估计正是区间的中点。

2 Cotes求积公式的改进及误差估计

定理2 若f (x) ∈C $_{[a, b]}^{6}$ , 则求积公式 (为了方便, 称为改进Cotes求积公式)

$\begin{array}{l} \int_{a}^{b} f (x) d x \approx \frac{b - a}{90} [7 f (a) + 32 f (a + h) + \\ 12 f (\frac{a + b}{2}) + 32 f (b - h) + 7 f (b)] - \\ \frac{8}{945} (\frac{b - a}{4})^{7} f^{(6)} (\frac{a + b}{2}) \end{array}$

具有7次代数精确度。

证明若f (x) ∈C $_{[a, b]}^{6}$ , Cotes求积公式有余项估计 $R (f) = - \frac{8}{945} (\frac{b - a}{4})^{7} f^{(6)} (ξ)$ , 利用余项估计, 确定求积公式

$\begin{array}{l} \int_{a}^{b} f (x) d x \approx \frac{b - a}{90} [7 f (a) + 32 f (a + h) + \\ 12 f (\frac{a + b}{2}) + 32 f (b - h) + 7 f (b)] - \end{array}$

$\frac{8}{945} (\frac{b - a}{4})^{7} f^{(6)} (η)$ 。

其中η=a+λ (b-a) , 确定λ∈ (0, 1) 使上式具有尽可能高的代数精确度就能够给出Cotes求积公式的改进。利用1节中的余项的渐进估计, 可知λ=1/2时, 可以提高求积公式的代数精确度。当λ=1/2时, 通过简单计算可以证明定理2中的求积公式具有7次代数精确度。

下面进一步讨论改进Cotes求积公式的误差估计, 有

定理3 若f (x) ∈C $_{[a, b]}^{7}$ , 则改进Cotes求积公式有误差估计

$| R_{c} (f) | \leq \frac{(b - a)^{8}}{3 870 720} Μ$ 。

其中M为 $| f^{(7)} (x) |$ 在区间[a, b]上的最大值。

证明由Cotes求积公式的余项估计有

$\begin{array}{l} R (f) = \int_{a}^{b} f (x) d x - \frac{b - a}{90} [7 f (a) + 32 f (a + h) + \\ 12 f (\frac{a + b}{2}) + 32 f (b - h) + 7 f (b)] = \end{array}$

$- \frac{8}{945} (\frac{b - a}{4})^{7} f^{(6)} (ξ)$ 。

从而改进Cotes求积公式的余项为

$\begin{array}{l} R_{c} (f) = - \frac{8}{945} (\frac{b - a}{4})^{7} (f^{(6)} (ξ) - f^{(6)} (\frac{a + b}{2})) = \\ - \frac{8}{945} (\frac{b - a}{4})^{7} f^{(7)} (η) (ξ - \frac{a + b}{2}) \\ | R_{c} (f) | \leq \frac{8}{945} (\frac{b - a}{4})^{7} | f^{(7)} (η) (ξ - \frac{a + b}{2}) | \leq \frac{(b - a)^{8}}{3 870 720} Μ \leq \frac{(b - a)^{8}}{3 870 720} Μ \end{array}$

其中M为 $| f^{(7)} (x) |$ 在区间[a, b]上的最大值。

3 结论

通过对数值积分公式的余项的渐进估计进行分析, 进而利用渐进估计对数值积分公式进行修正和改进, 达到提高求积公式的代数精确度的目的。本文给出了Cotes求积公式余项的渐进性质以及改进的Cotes求积公式, 类似地可以讨论Simpson3/8求积公式渐进性质及其改进。

参考文献

[1]许晓阳, 陈露.两类数值积分公式的改进.科学技术与工程, 2008;8 (5) :1294—1295

[2]邱淑芳, 王泽文.数值积分公式中间点的渐进性质及其应用.数学的实践与认识, 2006;36 (5) :218—223

数值积分篇4

基于蒙特卡罗方法和区域分解法,建立低地球轨道空间环境航天器表面原子氧通量密度和积分通量的数学模型.模型考虑了航天器表面几何构型、原子氧数密度和分析热运动、地球自转对航天器速度的影响以及轨道运行参数.通量密度分布的求解是通过其微分方程的对于独立变量分子运动速度和与表面速度矢量合成的积分得到,积分通量是通过沿轨道时间积分来实现.与此同时,得到了沿入射攻角变化原子氧分布的`最大值和最小值.计算结果表明:通量分布伴随入射攻角增大而急剧下降,在迎风面达到最大值,背风面最小值.入射攻角是影响分布计算结果的重要因素.计算误差与NASA-LDEF飞行试验实验结果吻合较好.

作者：刘阳姜利祥李涛 LiuYang JiangLixiang LiTao 作者单位：刘阳,LiuYang(大连海事大学,大连116026;清华大学,北京100084)

姜利祥,李涛,JiangLixiang,LiTao(北京卫星环境工程研究所,北京,100094)

数值积分篇5

关键词：椭圆积分,转绳线,数值求解,垂直轴风力发电机

椭圆积分是从事机械制造(例如柔轮或非圆齿轮等零件的设计)、断裂力学、电磁学、回旋加速器、精密计量以及天体力学等方面工作的科技和教学人员常用的数学工具[1,2],其数值往往需要查表才能得到[2]。

尽管大部分的数学手册或数学用表都会给出完全椭圆积分的数值,但由于受到位数和步长的限制大都不可能达到很高的精确度[2]。针对椭圆积分的近似求解,文献[2]做了相应研究。与文献[2]不同的是,文中应用数值方法对第一类完全椭圆积分和不完全椭圆积分进行近似的求解,结合应用实例,对特定值k值的近似结果的误差进行计算分析,所采用的方法简单、可靠,避免了查表的麻烦以及先进的计算工具与落后的查表方法之间的不协调,同时也可以大大地减少数据的误差,方便在工程设计人员中进行推广。

1 椭圆方程的数值求解及误差分析

考虑到合理的求解速度,针对椭圆积分的数值近似求解选择了收敛速度较快的梯形法。通用的梯形公式表达为:

由文献[3]得到梯形公式的最大误差计算公式为:

由式(2)可以看出,减小误差的最为直接的方法就是取n→∞,而实际计算过程中,如果n值大小不做限制,虽然得到的结果趋于精确但计算耗时太长,因此选择合适的步长和迭代次数,得到合理精度的计算结果对工程设计人员来说是具有实用意义的。

由梯形公式法,分别取值n=10i(其中i=1,2,3,…7),进行计算,并使用式(2)对计算结果进行误差计算。随着n值的增加,计算机对椭圆积分的计算时间亦是成几何级数增加,从对应的结果可以看出,如果最大计算误差<0.01%为可接受的范围,则对应10-5的步长n-1可以最快得到近似的结果;对于步长≤10-6的近似计算来说,过小的步长对近似解的精度影响有限。图1～3为关键计算步长时的完全椭圆积分数值计算结果的误差。

2 数值解法的应用实例

垂直轴风力发电机(Vertical Axis Wind Turbine,VAWT)的风轮从叶片型式上可以分为圆柱形达里厄风轮(H型风轮)、抛物线达里厄风轮(准型风轮)和圆台型达里厄风轮,见图4。

通常φ型风轮采用Troposkien曲线作为叶片型线,如图5示。总长L0,单位长度质量m0的柔性绳索2端固定,以角速度ω绕固定轴旋转,受离心力作用自然弯曲形成Troposkien曲线,取PP0段为研究对象,P点处绳索的张力为T,PP0

段总重量为G,离心力为Fc,将PP0段的载荷进行坐标轴正交分解,得到:

式中:δ-张力T与垂直方向的夹角,(°),如图标,取值范围δ∈[0,π/2]。

式(3)中的Fc和G分别为:

式中:m0-单位长度叶片的质量,kg/m;g=9.83m/s2(重力加速度);ω-风轮旋转角速度,rad/s;T0-材料理论极限张力,N。

由图5所示的几何关系及综合式(3)～(5),得到:

式(6)为“微分-积分方程”。当ω较大时,绳索的自重相对于张力、离心力是小量,对其忽略后得到:

利用H0对r和h无因次化,得到:

取值:

对应段,取值β=R/H0,整理式(9)得到:

整理多项式后,引入变量式(11)和式(12),得到式(13)

取值x=sinψ,ψ∈,则dx=cosψdψ,替代式(13)中变数x,得到:

式(14)整理后得到:

式(15)为第一类椭圆积分的不完全椭圆积分与椭圆积分的差值。对应式(15)中的完全椭圆积分与不完全椭圆积分的差值,应用梯度公式法进行数值计算,其中取步长10-5,取k=sinψ,设定H0=27m,应用梯形公式法对式(15)在ψ1=30°和ψ2=60°时进行近似计算,结果表示为图6。

取式(16)对应ψ的值为完全椭圆积分函数的max(Mj),对应ψ1和ψ2应用梯形公式法的最大计算误差为:

3 结论

文中详细列出了不同k值和步长的完全椭圆积分的数值计算结果,并与数学用表比对,计算了与数学用表的误差和梯形公式法的理论误差。分析表明:(1)运用数值解法求解积分,步长选择对结果精度和计算耗时均有影响;(2)梯形公式法适用于第一类椭圆积分的数值计算,步长取值10-5时,与数学用表的误差小于0.01%。

椭圆积分的近似求解基于通用平台Matlab实现,可根据需要的精度调节计算步长,在计算耗时和精度间选择合适的方案。文中提出的方法简单、可靠,适用于类似积分的数值计算,方便在工程设计人员中推广。

参考文献

[1]Hancock H.Elliptic Integrals[M].New York:Dover Publications,Inc.,1958.

[2]何明高.在Matlab计算机语言程序中调用勒让德第一类完全椭圆积分[J].广东机械学院学报,1996(12):42-47.

[3]四川矿业学院数学教研组.数学手册[M].北京:科学出版社,1977.

数值积分篇6

在工程技术和科学计算中,经常会遇到求解数值积分∫ $_{a}^{b}$ f(x)dx,除了几种有限类型外,被积函数f(x)往往很难甚至无法求出它的原函数,即使有的原函数能求出来,但是也相当麻烦不适宜于计算,这时常常需要借助近似方法求解近似数值解。数值积分虽然计算方法很多,如梯形求积方法[1]、抛物线求积方法[1]、Newton-Cotes方法、Romberg 方法、Gauss 方法等[2],但这些传统方法都有一定的局限性。其中梯形求积方法和抛物线求积方法适用于光滑性较差的被积函数,但是对于其它被积函数精度较低;Newton-Cotes方法是一种利用插值多项式来构造数值积分的常用方法,但高阶Newton-Cotes方法的收敛性没有保证;Romberg方法虽然收敛速度快,计算精度较高,但是计算量较大;Gauss方法积分精度高、数值稳定、收敛速度较快,但节点与系数的计算较困难。针对这些问题,本文基于对基本人工鱼群算法进行改进提出一种基于人工鱼群算法的优化分割数值积分算法。通过仿真实验与传统的数值积分方法比较表明新算法十分有效,可以看作是传统方法的补充和拓展,在工程中有较大的应用价值。

1 人工鱼群算法

1.1 基本人工鱼群算法

在一片水域中,鱼往往能自行或尾随其它鱼找到营养物质多的地方,因而鱼生存数目最多的地方一般就是本水域中营养物质最多的地方。人工鱼群算法AFSA(Artificial Fish-school Algorithm)[3]就是根据鱼群上述特点提出的一种新型仿生优化算法,通过构造人工鱼来模仿鱼群的觅食、聚群、追尾及随机行为,从而实现寻优,是群体智能思想[4]的一个具体应用。它的主要特点是不需要了解问题的特殊信息,只需要对问题进行优劣的比较,通过各人工鱼个体的局部寻优行为,最终使全局最优值突现出来。作为一种新型的全局寻优策略,它已应用于时变系统的在线辩识、鲁棒PID的参数整定和优化前向神经网络中,取得了较好的效果。

然而随着优化问题的复杂性,基本AFSA在应用中也存在不足:当寻优的空间维数高或区域较大或处于变化平坦的区域时收敛于全局最优解的速度减慢,搜索性能劣化,甚至会陷入局部极小;算法一般在优化初期收敛快,后期却往往较慢;得到的解是满意解,精度不高。因此在具体应用中常常要对它进行改进。

1.2 柯西变异人工鱼群算法

Gunter[5]通过理论分析表明柯西变异算子具有较强的局部逃逸能力。为了增强人工鱼群算法在高维复杂环境下的全局搜索能力和提高搜索效率,本文在基本人工鱼群算法中引入柯西变异算子来对人工鱼的状态Xi=(xi1,xi2,...,xin)进行变异,定义如下:

X′i=Xi·[1+k·Cauchy(0,1)] (1)

其中Cauchy(0,1)的一维概率密度函数为 $f (x) = \frac{1}{π} \cdot \frac{1}{1 + x^{2}}, - \infty < x < + \infty$ ,k是随迭代由1线性减为0的变量,迭代前期变异幅度大,后期变异幅度小。

设ξ是[0, 1]上的服从均匀分布的随机变量,由随机变量生成函数定理[6],Cauchy(0,1)的柯西分布的随机变量生成函数为 η=tan[(ξ-0.5)π]。

为判断随迭代次数增加搜索结果是否有改进,种群迭代后将种群中最优鱼的状态的函数值与公告板[7]进行比较,如果优于公告板,则取代公告板状态。当在连续两次迭代过程中公告板没有改变或变化极小时,此时我们看作鱼群迭代停滞,则启用变异操作。此算法CMAFSA(Artificial Fish-school Algorithm Based on Cauchy Mutation)流程如下:

Step1 初始化群体:在控制变量可行域内随机生成Np条人工鱼,形成初始鱼群。

Step2 公告板赋初值:计算初始鱼群各人工鱼当前状态的函数值y,取y为最小值者进入公告板,并将此鱼位置状态也赋值给公告板。

Step3 行为选择:各人工鱼分别模拟追尾行为和聚群行为,评价行动后的值,选择y值较小的行为实际执行,缺省行为方式为觅食行为。

Step4 公告板更新:种群迭代后,如果种群最优鱼状态的y优于公告板的y,则更新公告板。

Step5 变异条件判断:若公告板在连续两次迭代过程没有改变或变化极小(<η),则执行Step6;否则执行Step7。

Step6 变异操作:用历史最优鱼替换当前群体的最差鱼形成中间种群,对中间种群的人工鱼按式(1)变异,计算变异后各状态的函数值与公告板比较,若优,更新公告板。

Step7 终止条件判断:判断是否已达到预置的最大迭代次数MaxIteration或判断最优值是否达到了满意的误差界内(<ε),若不满足,执行Step3,进行下一代鱼群优化过程;否则执行Step8。

Step8 算法终止:输出最优解(即公告板中人工鱼状态和函数值)。

1.3 CMAFSA性能测试

本文用一个典型的高维(30维)、多峰函数Ackley来测试CMAFSA的性能。算法参数设置如下:Np=100,Try_number=30,δ=0.618, η=10-4,最大迭代次数200次, Visual=5,Step=0.5。独立测试5次(随机生成5个初始群体分别用AFSA和CMAFSA运行), 运行结果如表1和图1(为了便于比较, 图中纵坐标是5次的每代最优值的平均值的常用对数)所示。

从表1和图1可以看出AFSA在优化高维、多峰问题时效果很差,而CMAFSA的精度、收敛速度、稳定性有了非常明显的提高,实际上5次实验中CMAFSA最差的一次是在第114(从图1也可看出)代达到了理论最优值。

2基于人工鱼群算法的优化分割数值积分算法

在数值积分 $\int_{a}^{b} f (x) d x \approx \sum_{i = 1}^{n} f (ξ_{i}) Δ x_{i}$ 的近似计算中,分割的优劣和f(ξi)的替代对求解结果的精度有很大的影响。我们的基本思想是在积分区间[a,b]上用柯西变异人工鱼群算法搜到一个最优的分割,考虑到函数值的变化关系和曲线的弯曲情况,在得到的最优分割对应的小区间上用一个合适的值代替f(ξi)来达到近似计算,算法步骤如下:

Step1 确定人工鱼个体的表达式:Xi=(x1,x2,...,xn),其中xi为积分区间[a,b]内的随机n个节点。

Step2 鱼群初始化:在[a,b]内随机生成Np条人工鱼, 形成初始鱼群。实际上就是在[a,b]上产生Np个随机分割。

Step3 公告板赋初值:计算初始鱼群各人工鱼当前状态的函数值y,取y为最小值者进入公告板,并将此鱼位置状态也赋值给公告板。其中函数值是这样计算的:将人工鱼的状态值置于积分区间的左端点和右端点之间,各自按照升序排好顺序。这样,积分区间[a,b]总共有n+2个节点和n+1个小区间,然后分别计算这n+2个节点相邻节点之间的长度Δxi,i=1, 2, …,n+1。再计算出这n+2个节点对应的函数值和每个小区间中点的函数值。找出每个小区间左端点、中点和右端点的函数值中的最小值wi和最大值Wi,i=1,2,…,n+1,则该个体所对应的函数值 $y = \sum_{i = 1}^{n + 1} | W_{i} - w_{i} | Δ x_{i}$ 。

Step4 行为选择:各人工鱼分别模拟追尾行为和聚群行为,评价行动后的值,选择y值较小的行为实际执行,缺省行为方式为觅食行为。

Step5 公告板更新:种群迭代后,如果种群最优鱼状态的y优于公告板的y,则更新公告板。

Step6 变异条件判断:若公告板在连续两次迭代过程没有改变或变化极小(<η),则执行Step7;否则执行Step8。

Step7 变异操作:用历史最优鱼替换当前群体的最差鱼形成中间种群,对中间种群的人工鱼按式(1)变异,计算变异后各状态的函数值与公告板比较,若优,更新公告板。

Step8 鱼群算法终止条件判断:判断是否已达到预置的最大迭代次数MaxIteration或判断y的最优值是否达到了满意的误差界内(<ε),若不满足,执行Step4,进行下一代鱼群优化过程;否则执行Step9。

Step9 得到最优人工鱼个体状态(xbest1,xbest2,...,xbestn),这n个数同a,b一起排序,从而得到最优分割。

Step10 输出积分 $\int_{a}^{b} f (x) d x \approx \sum_{i = 1}^{n + 1} \bar{f}_{i} \cdot Δ x_{b e s t_{i}}$ 。其中 $\bar{f}_{i}$ 是最优分割所得到的n+1个小区间每个区间左、右端点和三点对应的函数值的加权平均 $\bar{f}_{i} = \frac{f_{i 左} + 4 f_{中} + f_{i 右}}{6}$ ,Δxbesti是最优分割所对应小区间的长度。

3 积分仿真实验与分析

为了检验本文提出的数值积分算法的有效性和正确性,本文给出的一些实例,以下仿真均在Matlab 7.0下编程运行。参数设置如下:Np=50,Try_number=20,δ=0.618, η=10-4,ε=10-4,最大迭代次数200次,Visual=0.1,Step=0.01, n=100。每个函数独立运行5次,并与传统的梯形法、Simpson方法、Composite Simpson方法、Romberg方法和神经网络等比较。

例1 分别计算被积函数 $x^{2}, x^{4}, \sqrt{1 + x^{2}}, \frac{1}{1 + x}, \sin (x), e^{x}$ 六个函数在[0,2]上的积分。与文献[8,9]的结果统计比较如表2所示。

由表2看出传统方法中只有Simpson法在求解函数x2时的结果优外,其他的在这6个函数上效果都不好,而本文的算法对每个函数均独立运行5次皆是精确值,可见精度高且很稳定。

例2 计算积分:

∫ $_{0}^{1}$ e-x2dx

被积函数的原函数不是初等函数,本文算法(参数设置同前)5次运算与文献[9]分别用矩形法、梯形法、Simpson法给出结果统计如表3所示。

由表3看出对于被积函数的原函数不是初等函数的函数,本文的算法也是相当有效的。

例3 计算积分:

$\int_{0}^{48} \sqrt{1 + (c o s x)^{2}} d x$

用传统的方法Romberg比较麻烦,不易处理,文献[9]用神经网络算法得到的结果是58.5205,文献[10]用Composite Simpson rule计算的结果是58.47082,而此积分的精确值是58.4704691。由于此被积函数是一个以π为周期的函数,48=15π+(48-15π),此积分可化为:

$\int_{0}^{48} \sqrt{1 + (c o s x)^{2}} d x = 15 \int_{0}^{π} \sqrt{1 + (c o s x)^{2}} d x + \int_{0}^{48 - 15 π} \sqrt{1 + (c o s x)^{2}} d x$

用本文的算法(参数设置同前)运行5次求得的积分值的最差值是58.47046674929646,最佳值是58.47046909676099,平均值是58.47046871265326,由此看出对于Romberg方法不易处理的函数积分,本文的算法也是很优的。

例4 计算奇异函数:

$f (x) = {\begin{cases} e^{- x} 0 \leq x < 1 \\ e^{- \frac{x}{2}} 1 \leq x < 2 \\ e^{- \frac{x}{3}} 2 \leq x \leq 3 \end{cases}$

在[0,3]上的积分该函数的积分精确值是1.546036,文献[9]用神经网络算法算得的结果是1.5467。用本文的算法(参数设置同前)运行5次求得的积分值的最差值是1.54554025335573,最佳值是1.54606116448507,平均值是1.54600370898978,可见本文算法也可求解奇异函数的积分。

4 结论

本文在分析传统数值积分和人工鱼群算法的基础上提出了一种求解定积分数值解的新算法,此算法对被积函数没有太多的限制,通过四组积分算例实验表明此方法十分有效,此方法可以看作是传统方法的补充和拓展。在优化区间的每一段上如何再处理以及用此算法计算重积分是我们下一段要做的工作。

摘要：在分析传统求数值积分和基本人工鱼群算法不足的基础上,提出了一种基于人工鱼群算法的优化分割数值积分算法,该算法不仅能求解通常意义下函数的数值积分,还能计算奇异函数的数值积分。通过算例与传统数值积分方法比较,实验结果表明该算法是可行的和有效的。

关键词：人工鱼群算法,柯西变异,优化分割,数值积分

参考文献

[1]华东师范大学数学系.数学分析[M].北京:高等教育出版社,2001.

[2]薛密.数值数学和计算[M].上海:复旦大学出版社,1991.

[3]李晓磊,邵之江,钱积新.一种基于动物自治体的寻优模式:鱼群算法[J].系统工程理论与实践,2002,22(11):32-38.

[4]Bonabeau E,Theraulaz G.SwarmSmarts[J].Scientific American,2000,282(3):72-79.

[5]Gunter R.Local convergence rates of simple evolutionary algorithms withCauchy mutation[J].IEEE Transactions on Evolutionary computation,1997,4(1):249-258.

[6]王梓坤.概率论基础及其应用[M].北京:科技出版社,1979.

[7]李晓磊,路飞,田国会,等.组合优化问题的人工鱼群算法应用[J].山东大学学报:工学版,2004,34(5):64-67.

[8]Burden R L,Faires J D.Numerical Analysis[M].7th Edition.Brooks/Cole,Thomson Learning,Inc,2001:190.

[9]王小华,何怡刚,曾喆昭.三角基函数神经网络算法在数值积分中的应用[J].电子与信息学报,2004,26(3).

数值积分篇7

然而,由于飞机外形结构复杂,采用笛卡尔网格对飞机表面进行网格划分时,所出现的数值积分区域就不是标准的积分区域了。例如二维翼型网格划分后会出现带曲边边界的三角形和四边形网格,三维翼型网格划分后会出现带曲面的四面体网格。

目前,一般的数值积分技术中,对于带曲边边界的网格是直接线段化处理,然而这样做精度有限,而如果将网格加密处理,精度提高有限却会增加计算负担和降低计算效率,所以在本文将提出一种保留曲边边界进行处理并且提高了数值积分精度的方法。

1 复杂网格介绍

网格划分属于流体力学计算中的预处理部分,是将计算区域划分为较小的、不重叠的子域或单元(网格),从而得到有限单元及节点。而单元的数目决定了CFD的计算精度,网格的细密程度决定了数值解的精度和计算时间。

在采用笛卡尔网格进行划分时,由于飞机外形复杂,例如对飞机机翼表面进行网格划分时,会出现一边为曲边的四边形网格(或三角形网格)和两边为曲边的四边形网格(或三角形网格),而如今随着不断研究,飞机外形也是日趋复杂化,相应地对飞机表面进行网格划分所得到的网格也是日渐复杂。

当然,若是网格划分足够密,每个网格足够小,那么曲边在一定程度上也就相当于直边了,所以在目前的数值积分技术中,都是采用直接将曲边看成直线进行处理,但如此做法必然会影响到最后的数值积分精度,所以本文提出一种基于复杂网格处理的数值积分技术,来提高数值积分的精度。

2 复杂网格处理

现假设在二维平面 ℝ2上有一个带曲边边界的四边形网格,要求解函数f (x,y) 在其积分区域 Ω 上的数值积分,一般情况下,直接将曲边看成直线边进行处理,利用上面介绍的数值积分方法(Simpson或Gauss方法等)求数值积分得到所要的结果,下面来介绍一下我们的处理方法。

首先建立一般化的直角坐标,坐标平面为x Oy平面,四边形网格的四个顶点坐标所对应的向量依次为 γ1,γ2,γ3,γ4,其中向量 γ2,γ4所对应的顶点连接的正是网格的曲边。现有另一坐标系O - ξη ,其上有一个单位为1 的正方形区域 Ωˉ ,接下来找到一个一一映射:

将区域一一映射到区域Ω上。

很自然地,有如下映射关系:

那么就可以构造一个一一映射:

则有:

根据数学分析中积分变换的内容有如下关系:

其中 Ω 表示原来曲边边界网格的区域,表示单位为1的标准正方形区域,而J(φ) 为映射 φ 的Jacobi矩阵,

我们还可进行如下处理:

根据一一映射关系,可以假设 ℝ2(x,y) 中曲边的曲线函数为 Γ(ξ)(0 ≤ ξ ≤ 1) ,那么在曲边的两个端点处,就有如下关系:

则在式(8)中就有

其中

通过对 ξ 的取值讨论,容易得到上式右边等于,于是式(8)就可以写成:

对带曲边的网格区域进行如上处理后再利用已有的数值积分技术进行数值积分即可得到精度提高的数值结果了。

3 数值实验和结果

为了便于进行实验,我们取了一个相对简单的算例,并且能得到精确结果的例子,最后能够进行多方对比来验证我们的方法。

现假设有一个带曲边的四边形网格,并且曲边函数为f (x)= ex,另外三边分别为x = 0,y = 0,x = 1,区域表示为 Ω 。

现在要求解的数值积分结果,其中F(x,y)=x2+y2,我们都知道该结果完全可以精确求解出来:

数值结果为I ≈ 2.8388970421。

而在现有的数值积分技术中,都是忽略了曲边,直接连接了曲边的两点,而后利用现有的数值积分方法进行数值计算。我们选取了比较常用的Simpson方法和Gauss方法进行数值积分,得到结果如下:

下面我们先用本文第三部分介绍的方法对曲边网格进行处理后有:

则Jacobi矩阵为

则所求的数值积分变成了

其中表示标准正面形区域。

之后我们依然采用现有的Simpson方法和Gauss方法,并且取同样的节点进行数值积分计算,得到的结果为:

IS′= 2.8384170008,IG′= 2.8388960361.

我们可以对以上数据进行简单的处理和误差分析,得到结果如表1所示。

其中 α 表示对于网格曲边直接连接曲边两点做线段化处理,而 β 表示对网格曲边进行映射变换处理,处理方法见本文第三部分。

从上表可以很容易看出,对曲边边界进行映射变换后再进行数值积分的话,精确度有显著提高,尽管这只是一个简单算例,但也说明了我们对于带曲边的网格进行适当的处理再数值积分的话,的确可以提高精度,相信这种处理方法和思想在高精度数值积分技术中会发挥其应有的作用。

参考文献

[1]徐利治,周蕴时.高维数值积分[M].科学出版社,1980.

[2]徐利治,周蕴时,何天晓.高维数值积分选讲[M].安徽教育出版社,1985.

[3]Philip J.Davis,Philip Rabinowitz著,冯振兴,伍富良译.数值积分法[M].高等教育出版社,1986.

数值积分篇8

本文对Matlab中的各种数值积分函数进行了说明和分析, 并且通过对不同函数的比较总结了各个数值积分函数的优缺点。

1. Matlab数值积分函数的分析

1.1 基于梯形算法的积分函数

在Matlab中trapz (x, y, dim) 表示梯形法数值积分, 通过已知参数x, y按dim维使用梯形公式进行积分。向量x和y有相同的长度, (xi, yi) 代表曲线上的一点。曲线上点的距离不一定相等, x值也不一定有序。

另外, Matlab中函数cumtrapz (x, y, dim) 表示累积的梯形积分, 即不断的从第一个累积到当前的结果, 它的返回值是与x维数相同的向量。

1.2 基于Simpson算法的积分函数

quad (fun, a, b, tol, trace) 用自适应Simpson公式[3]进行数值积分的, 适用于精度要求低、平滑性较差的被积函数。其参数中, 被积函数fun必须是函数句柄;积分限[a, b]必须是有限的, 因此不能是inf;参数tol表示绝对误差, 默认值为10-6;当trace是非零值时则在递归过程中显示[fcnt a b-a Q]的值。值得注意的是, 该算法当积分递归计算超过10000次或当函数有不可积的奇点时, 则返回inf。

另外, Matlab中函数quadv (fun, a, b, tol, trace) 是用矢量化自适应simpson法进行数值积分的, 该函数是将quad函数矢量化, 然后一次可以计算多个积分, 所有的要求完全与quad相同。

1.3 基于Lobatto算法的积分函数

函数quadl (fun, a, b, tol, trace) 是用自适应Lobatto法进行数值积分的, 适用于精度要求高, 被积函数曲线比较光滑的数值积分。其中参数说明及注意事项与quad函数相同。

1.4 基于Gauss-Kronrod算法的积分函数

函数quadgk (fun, a, b, param1, val1, param2, val2, ...) 是用自适应Gauss-Kronrod法进行数值积分, 适用于高精度和震荡数值积分, 支持无穷区间, 并且能够处理端点包含奇点的情况, 同时还支持沿着不连续函数积分, 复数域线性路径的围道积分法。其参数中, 被积函数fun必须是函数句柄, 积分限[a, b]可以是[-inf, inf], 但必须快速衰减, param, val为函数的其它控制参数。

2. 几种数值积分函数的比较

从上述六个函数中选取具有代表性的trapz、quad、quadl以及quadgk四个函数进行分析。

2.1 基于精度的分析

使用这四个函数分别计算积分, 积分区间取[0, 1], 在Matlab中运行, 得出Matlab命令与计算结果, 见表 (2-1) 。

由表可看出, 函数trapz是定步长计算积分的, 它可以通过减小步长来提高精度, 但无法由该方法直接知道当前步长下的精度。而后面三个函数是由精度来控制运算的, quad及quadl是利用递归算法, 当达到指定精度时停止计算。由Matlab中这三种算法本身的实现过程可知, quad及quadl算法的默认精度都是10-6, 而quadgk的默认精度为10-10。

另外, 由于函数quad是用simpson算法实现的, 即用抛物线来代替被积函数, 只有三次代数精度[4], 所以对于光滑性[5]较差的函数具有较好的运算结果;而函数quadl是用Lobatto算法实现的, 故更适用于光滑性较好的函数。

2.2 基于运算速度的分析

由上面的研究可知, 用quad、quadl及quadgk进行数值计算时, 对运算结果的精度都能有较好的控制。为了进一步研究这三个函数的优缺点, 下面我们研究它们的运算速度随精度的变化。改变这三个函数的参数设置, 让它们的误差从10-4变化到10-16, 记录耗时随精度的变化情况, 得出图 (2-1) 。

分析图2-1发现, 当对精度的要求较低时, 利用函数quad计算积分耗时最短, 而随着精度增加, 它的耗时增加很快。故当精度要求较高时, quad函数的效率较低, 不建议继续使用。

观察函数quadl及quadgk (2009) 的耗时随精度的变化曲线发现的, 这两个函数的运算速度基本稳定, 但quadl当精度高于10-13时, 耗时在逐步增加。故当精度要求在10-6与10-13之间时, 建议读者优先考虑使用quadl函数计算积分。

2.3 基于适用范围的分析

由于函数quadgk在Matlab2009版才推出, 所以普及程度不够高, 但无论从精度、速度或是适用范围来看, 它是目前Matlab中最强大的数值积分函数。

3. 结论

本文对Matlab中的几种数值积分函数进行了简单阐述及说明, 对函数trapz、quad、quadl及quadgk进行具体分析及比较, 得出如下结论:

总的来说, 在Matlab的几种数值积分函数里, 函数quadgk的功能最为强大。

摘要：Matlab软件被广泛应用于科学计算, 而数值积分又是很重要的一个分支。本文首先对Matlab中的几种数值积分函数进行阐述和分析, 然后通过举例对函数trapz、quad、quadl及quadgk进行具体比较, 总结了各个数值积分函数的优缺点及适用性。

关键词：Matlab,数值积分函数,适用范围

参考文献

[1]张志涌, 等编著.精通MATLAB R2011a[M].北京:北京航空航天大学出版社, 2011, 11.

[2]华东师范大学数学系编著.数学分析上册[M].北京:高等教育出版社, 2011.