编程实现数值积分的几种--方法 c语言

编程实现数值积分的几种--方法c语言
数值计算2010-11-05 09:52:43 阅读385 评论1 字号：大中小订阅
复化梯形公式
在区间不大时, 用梯形公式、辛卜生公式计算定积分是简单实用的, 但当区间较大时, 用梯形公式、辛卜生公式计算定积分达不到精确度要求 . 为了提高计算的精确度，我们将[a,b] 区间n等分，在每个小区间上应用梯形公式、辛卜生公式计算定积分，然后将其结果相加，这样就得到了复化梯形公式和复化辛卜生公式。

1. 复化梯形公式
将积分区间等分, 设, 则节点为
对每个小区间上应用梯形公式, 然后将其结果相加，则得
(3.14)
称(3.14) 式为复化梯形公式 .
当在[a,b] 上有连续的二阶导数时，则复化梯形公式(3.14) 的余项推导如下：
因为
所以在区间[a,b] 上公式(3.14) 的误差为
又因为在区间[a,b] 上连续，由连续函数的性质知，在区间[a,b] 上存在一点，
于是
（ 3.15 ）
复化梯形公式，复化抛物线公式和Romberg求积法的算法程序：以下程序均定义误差限为1*10^-5;
1)复化梯形公式：
#include <stdio.h>
#include <math.h>
#define e 1e-5
#define a 0 //积分下限a
#define b 1 //积分上限b
#define f(x) (4/(1+(x*x))) //被积函数f(x)
int main()
{
int i,n;
double h,t0,t,g;
n=1; //赋初值
h=(double)(b-a)/2;
t=h*(f(a)+f(b));
do
{
t0=t;
g=0;
for (i=1;i<=n;i++)
g+=f((a+(2*i-1)*h));
t=(t0/2)+(h*g); //复化梯形公式
n*=2;
h/=2;
}
while (fabs(t-t0)>e); //自定义误差限e printf("%.8lf",t); //输出积分的近似值
return 0;
}
2)复化抛物线公式：
#include <stdio.h>
#include <math.h>
#define e 1e-5
#define a 0 //积分下限a
#define b 1 //积分上限b
#define f(x) (4/(1+(x*x))) //被积函数f(x)
int main()
{
int i,n;
double f1,f2,f3,h,s0,s;
f1=f(a)+f(b); //赋初值
f2=f(((double)(b+a)/2));
f3=0;
s=((double)(b-a)/6)*(f1+4*f2);
n=2;
h=(double)(b-a)/4;
do //复化抛物线算法
{
f2+=f3;
s0=s;
f3=0;
for (i=1;i<=n;i++)
f3+=f((a+(2*i-1)*h));
s=(h/3)*(f1+2*f2+4*f3);
n*=2;
h/=2;
}
while (fabs(s-s0)>e); //自定义误差限
printf("%.8lf",s);
return 0;
}
3)Romberg求积法：
#include <stdio.h>
#include <math.h>
#define e 1e-5
#define a 0 //积分下限a
#define b 1 //积分上限b
#define f(x) (4/(1+(x*x))) //被积函数f(x)
double t[100][100];
int main()
{
int n,k,i,m;
double h,g,p;
h=(double)(b-a)/2;
t[0][0]=h*(f(a)+f(b));
k=1;
n=1;
do //Romberg算法
{
g=0;
for (i=1;i<=n;i++)
g+=f((a+((2*i-1)*h)));
t[k][0]=(t[k-1][0]/2)+(h*g);
for (m=1;m<=k;m++)
{
p=pow(4,(double)(m));
t[k-m][m]=(p*t[k-m+1][m-1]-t[k-m][m-1])/(p-1);
}
m-=1;
h/=2;
n*=2;
k+=1;
}
while (fabs(t[0][m]-t[0][m-1])>e); //自定义误差限e
printf("%.8lf",t[0][m]);
return 0;
}
给定精度，定义误差限为1*10^-5，分别求出步长的先验估计值：用复化梯形公式计算，要求h<0. 007746。

用复化抛物线公式计算，要求h<0.131607。

而实际上，在运用后验估计的程序中，以相同的精度，得到的复化梯形步长为h= 0.003906，得到的复化抛物线步长为h=0.0625，它们大致上分别为先验估
计值的一半，符合要求。

在实践中比较上体三种算法的计算量，当取相同精度1*10^-5时，复化梯形调用函数f257次，孵化抛物线公式调用函数f17次，Romberg算法调用函数f17
次，从计算量上看，后两者较小。

2.2实验2 一元微积分的编程实现
【实验目的与要求】
实验目的：
熟悉用Mathematic进行一元微积分计算的编程方法。

先修内容：
第一篇计算机数学第1章极限与连续和第2章微分与积分。

实验要求：
掌握数学表达式的正确书写格式；熟悉Mathematic有关一元微积分的常用命令、常用数学函数。

【实验原理】
Mathematic的基本语法、数学表达式的正确书写格式；有关一元微积分的常用命令、常用数学函数。

【实验步骤】
2.2.1实验内容1 极限
Mathematica 计算极限的命令是Limit 它的使用方法主要有表2.2.1种的一些命令。

表2.2.1极限的主要命令及说明
命 令 说 明
Limit[expr,x →x0] 当x 趋向于x0时求expr 的极限 Limit[expr,x →x0,Direction →1] 当x 趋向于x0时求expr 的左极限 Limit[expr,x →x0,Direction →-1] 当x 趋向于x0时求expr 的右极限 Infinity 无穷大
趋向的点可以是常数，也可以是+∞，-∞ 。

注意Mathemica 没有区分∞和+∞，求x →∞时的极限要小心。

下面就具体操作几个运行极限的Mathemica 程序。

1．求6
32lim 2-+∞→x x x 利用Limit[expr,x → Infinity]命令，计算∞→x lim expr ；再将表达式expr 转化成 6
322-+x x ，其中Sqrt[x^2+2]是指22+x 。

具体运行程序参见图2.2.1。

图2.2.1运行632lim 2-+∞→x x x 的Mathemica 程序
2．求220sin lim x
x x → 利用Limit[expr,x →0]命令，计算0lim →x expr ；再将表达式expr 转化成22sin x x ，其中Sin[x]^2是指sin 2x 。

具体运行程序参见图2.2.2。

图2.2.2运行220sin lim x
x x →的Mathemica 程序 3．求x x ln lim 0
+→ 利用Limit[expr,x →0,Direction →-1]命令，计算+→0
lim x expr ；再将表达式expr 转化成lnx ，其中Log[x]是指lnx 。

具体运行程序参见图2.2.3。

图2.2.3运行x x ln lim 0+
→的Mathemica 程序
2.2.2实验内容2 函数的微分
在Mathematica 中，计算函数的微分或是非常方便的，命令为D[f,x],表示对x 求函数f 的导数或偏导数。

该函数的常用格式有以下几种如表2.2.2。

表2.2.2关于函数微分的几个主要命令及说明
下面就具体操作几个运行微分的Mathemica程序。

例2.2.1求函数e x sinx的导数。

利用D[f,x]命令，计算f的导数；再将表达式f转化成e x sinx，其中Exp[x]是指e x。

这样就可得到函数e x sinx导数的Mathemica程序。

具体运行程序参见图2.2.4。

图2.2.4运行e x sinx的导数的Mathemica程序
例2.2.2求函数e x sinx的二阶导数
利用D[f,{x,n}]命令，计算f的n导数；再将表达式f换成e x sinx，其中Exp[x]是指e x；n换成2。

就可得到函数e x sinx的二阶导数的Mathemica程序。

具体运行程序参见图2.2.5。

图实验2.5函数e x sinx的二阶导数的Mathemica程序
例2.2.3假设a是常数，对sinax求导。

这题仍可以利用D[f,x]命令，计算f的导数；再将表达式f转化成sinax，其中a 不做任何处理，就可视作普通的字符（或称为常数）。

具体运行程序参见图2.2.5。

图2.2.5函数sinax的导数的Mathemica程序
Mathematica也可以求函数式未知的函数微分，通常结果使用数学上的表示法。

如下面例2.2.4。

例2.2.4求xg(x)对x的导数和4阶导数。

首先，求函数xg(x)对x的导数。

首先利用D[f,x]命令，将f转化成表达式xg(x)。

其运行结果是g[x]+xg'[x]。

具体运行程序参见图2.2.6。

再求函数xg(x)对x的4阶导数。

首先利用D[f,{x,n}]命令，将f转化成表达式xg(x)，n换成4。

其运行结果是4g（3）[x]+xg（4）[x]。

具体运行程序参见图2.2.6。

图2.2.6函数xg(x)对x的导数和4阶导数的Mathemica程序
对复合函数的求导上面的方法仍适用。

如例2.2.5。

例2.2.5求函数g[h[x]]对x求导。

为了求函数g[h[x]]对x求导，首先利用D[f,x]命令，将f转化成表达式g[h[x]]。

其运行结果是g' [h[x]]h ' [x]。

具体运行程序参见图2.2.7。

图2.2.7函数g[h[x]]对x求导的Mathemica程序
如果要得到函数在某一点的导数值可以把这点代入导数即可。

如例2.2.6中求f'(2)。

例2.2.6设f(x)=e x sinx,求f'(2)。

首先，利用D[f,x]/x→2命令求f'(2)。

再将f用表达式 Exp[x]*sinx替换，其中Exp[x]表示e x。

具体运行程序参见图2.2.8。

图2.2.8求f '(2) 的Mathemica 程序
2.2.3实验内容3计算积分
1.不定积分
在Mathematica 中计算不定积分命令为Integerate[f,x],当然也可使用工具栏直接输入不定积分式。

来求函数的不定积分。

当然并不是所有的不定积分都能求出来。

例如求⎰dx x )sin(sin Mathematica 就无能为力。

其运算结果见图2.2.9。

图2.2.9积分⎰dx x )sin(sin Mathematica 运算结果
但对于一些手工计算相当复杂的不定积分，MatheMatica 还是能轻易求得，例如下面例题2.2.7。

例2.2.7求⎰++du u u u 221121。

若用手工运算相当麻烦，但用MatheMatica 程序运算就容易多了。

先用工具栏直接输入不定积分式，在运行结果。

如图2.2.10的MatheMatica 程序设计结果。

图2.2.10积分⎰++du u u u 221121的Mathematica 运算结果
积分变量的形式也可以是一函数，例如例2.2.8中计算
⎰sin(sinx)d(sinx)也可以
利用MatheMatica 程序设计。

例2.2.8计算⎰sin(sinx)d(sinx)。

先用工具栏直接输入不定积分式，即可求得正确结果。

具体结果参见图2.2.11。

图2.2.11积分⎰sin(sinx)d(sinx)的Mathematica运算结果
对于在函数中出现的除积分变量外的符号，统统当作常数处理,例如例 2.2.9积分⎰(ax2+bx+c)dx中a、b、c在运算过程中就当常数处理的。

例2.2.9计算积分⎰(ax2+bx+c)dx。

首先，输入⎰(ax2+bx+c)dx的Mathematica运行程序后，即可以得到结果
cx+
22
bx
+
33
ax。

具体运行情况参见图2.2.12。

图2.2.12积分⎰(ax 2+bx+c)dx 的Mathematica 运算结果
2.定积分
定积分的求解主要命令也是用Integrate 只是要在命令中加入积分限
Integrate[f,{x,min,max}]，或者使用工具栏输入也可以。

例2.2.12 定积分⎰4
0x 2e ax
dx 。

利用命令Integrate[f,{x,min,max}]，将f 转化成表达式x 2e ax ，其中Exp[ax]表示e ax。

其运行结果是3128ax
e 。

具体运行程序参见图2.2.13。

图2.2.13定积分⎰40x 2e ax dx 的Mathematica 运算结果
命令Integrate[f,{x,min,max}]，也可以求广义积分。

例如下面例2.2.13。

例2.2.13 求积分⎰-402)2(1
dx x 。

这是一个瑕积分，照例用命令Integrate[f,{x,min,max}]直接输入即可得到结果。

详见图2.2.14。

图2.2.14积分⎰-402)2(1
dx x 的Mathematica 运算结果
同样，命令Integrate[f,{x,min,max}]也可求无穷积分。

例如例2.2.14。

例2.2.14 求积分⎰+∞141
dx x 。

这是一个广义积分，照例用命令Integrate[f,{x,min,max}]直接输入即可得到结果。

详见图2.2.15。

图2.2.15积分⎰+∞141
dx x 的Mathematica 运算结果
如果广义积发散也能给出结果，例如例2.2.15。

例2.2.15 求积分⎰-1121
dx x 。

这是一个发散的广义积分，照例用命令Integrate[f,{x,min,max}]直接输入也即可得到结果。

详见图2.2.16
图2.2.16积分⎰-1121
dx x 的Mathematica 运算结果
对于广义积分如果Mathematica 运算程序无法判定敛散性，就用给出一个提示，
例如例2.2.16，其运算结果参见图2.2.17。

例2.2.16 dx x
⎰2
01。

照例用命令Integrate[f,{x,min,max}]直接输入，因为x=0是奇异点，按一般广
义积分的定义，这个积分发散，这是出现提示，返回的只是原输入式的输出形式。

具体结果参见图2.2.17。

图2.2.17积分dx x
⎰201的Mathematica 运算结果 如果广义积分敛散性与某个符号的取值有关，它也能给出在不同情况下的积分结果例如例2.2.17。

例2.2.17 求积分⎰+∞11
dx x p 。

用命令Integrate[f,{x,min,max}]直接输入，+∞用Infinity 代替。

输出结果含有一个字母参数p ，这时返回的结果是一个条件表达式，当p 的实部大于1时值为
p
+-11，否则发散。

具体结果参见图2.2.18。

图2.2.18积分⎰+∞11dx
的Mathematica运算结果
x p
在Integrate中可加两个参数Assumptions 和 GenerateConditions例如例2.2.17中，只要用Assumptions→{Re[p]>1}就可以得到收敛情况的解。

图2.2.19积分⎰+∞11dx
在p>1的情况下的Mathematica运算结果
x p
3.数值积分
数值积分是解决求定积分的另一种有效的方法，它可以给出一个近似解。

特别是对于用Integrate命令无法求出的定积分，数值积分更是可以发挥巨大作用。

它的主要命令格式为表2.2.3中的命令。

表2.2.3 数值积分一些主要命令
下面我们求sinsinx在[0,π]上的积分值，由于这个函数的不定积分求不出，因此使用Integrate命令无法得到具体结果，但可以用数值积分可以求得近似值。

如例2.2.18。

例2.2.18求积分⎰π0sin(sinx)dx。

这个被积函数f(x)=sin(sinx)找不到初等的原函数，用命令Integrate[f,{x,min,max}]直接输入，得不到结果。

但用命令Nintegrate[f,{x,a,b}]输入，即可得到近似值。

具体情况参见图2.2.20，其中Pi表示π。

图2.2.20积分⎰π0sin(sinx)dx 在Mathematica 程序下的数值解
如果积分的被函数存在不连续点，或存在奇异点我们可对积分进行分段求解。

例如例
2.2.19。

例2.2.19 求积分⎰-11x 1
dx 。

被积函数x 1
在[-1，1]上，显然有x=0点是一个无穷间断点。

因此若要求其数值积分，
必须在其中插入点0。

于是，利用数值积分命令Nintegrate[f,{x,xmin ，x 1，x 2，……，xmax}]输入，其中x 1，x 2，……是奇异点。

这样就可得到积分的数值解。

具体操作参见图2.2.21
图2.2.21积分⎰-11x 1
dx 在Mathematica 程序下的数值解
对无穷积分，也可求得积分数值解，例如例2.2.20。

例2.2.20 求积分dx e x ⎰∞
-02。

对于求无穷积分的数值解，也可以直接用命令Nintegrate[f,{x,a,b}]，其中b 可以用无穷大Infinity 替换。

具体操作参见图2.2.22。

图2.2.22积分dx e x ⎰∞
-02
在Mathematica 程序下的数值解 思考题：
可以将第一篇中的习题一、习题二中的题用Mathematica 程序做一下。