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

1、第二篇 数学试验第2章数学试佥2.2实验2 一元微积分的编程实现25编程实现数值积分的几种-方法c语言数值计算 2010-11-05 09:52:43 阅读385评论1 字号：大中小订阅复化梯形公式在区间奇，8不大时，用梯形公式、辛卜生公式计算定积分是简单实用的 ，但当区间口，较大时，用梯形公 式、辛卜生公式计算定积分达不到精确度要求.为了提高计算的精确度，我们将 a,b区间n等分，在每个小区间上应用梯形公式、辛卜生公式计算定积分，然后将其结果相加，这样就得到了复化梯形公式和复化辛卜生公式。1.复化梯形公式将积分区间融等分，设翼,则节点为 =对每个小区间上应用梯形公式 ，然后将其结果相加，则得

2、，二f曲闻:M+1/（勺）十？'5。十+1/ 1J+（工J1二 4，（。）+ 2$_/（/）+/年）di-ll（3.14）称（3.14）式为复化梯形公式.当 J'lG 在a,b上有连续的二阶导数时,则复化梯形公式（3.14）的余项推导如下：因为*""暇/90与-1 M旅七4所以在区间a,b上公式（3.14）的误差为&=£&=-持£?电） JU11/ Jt-1又因为 J 1工）在区间a,b上连续，由连续函数的性质知，在区间 a,b上存在一点 肛,工28/7口 ="叮）a<r<h网U1于是网力=厅） 口

3、三町 0上12123（ 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);dot0=t;g=0;for (i

4、=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);自定义误差限eprintf("%.8lf",t);/输出积分的近似值return 0;2)复化抛物线公式: #include <stdio.h>#include <math.h>/积分下限a/积分上限b#define e 1e-5#define a 0#define b 1#define f(x) (4/(1+(x*x)/被积函数f(x)int main()int i,n;/赋

5、初值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 求积法：#incl

6、ude <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 t100100;int main()int n,k,i,m;double h,g,p;h=(double)(b-a)/2;t00=h*(f(a)+f(b);k=1;n=1;do/Romberg 算法g=0;for (i=1;i<=n;i+)g+=f(a+(2牛1)*h);tk0=(tk-1+(h*g)；for (m=1;m<

7、=k;m+)p=pow(4,(double)(m);tk-mm=(p*tk-m+1m-1-tk-mm-1)/(p-1);m-=1;h/=2;n*=2;k+=1;while (fabs(t0m-t0m-1)>e);自定义误差限 eprintf("%.8lf",t0m);return 0;给定精度，定义误差限为1*10八-5 ,分别求出步长的先验估计值：用复化梯形公式计算，要求h<0. 007746 。用复化抛物线公式计算，要求 h<0.131607 。而实际上，在运用后验估计的程序中， 以相同的精度，得到的复化梯形步长为h= 0.003906，得到的复化抛物

8、线步长为 h=0.0625 ,它们大致上分别为先验估计值的一半，符合要求。在实践中比较上体三种算法的计算量， 当取相同精度1*10A-5时，复化梯形调用函数f257次，孵化抛物线公式调用函数 f17次，Romberg算法调用函数f17次，从计算量上看，后两者较小。2.2实验2 一元微积分的编程实现【实验目的与要求】实验目的：熟悉用Mathematic进行一元微积分计算的编程方法。先修内容：第一篇 计算机数学 第1章极限与连续和第2章微分与积分。实验要求：掌握数学表达式的正确书写格式；熟悉Mathematic有关一元微积分的常用命令、常用数学函数。【实验原理】Mathematic的基本语法、数学

9、表达式的正确书写格式；有关一元微积分的常用命令、常 用数学函数。【实验步骤】2.2.1实验内容1极限Mathematica计算极限的命令是 Limit它的使用方法主要有表2.2.1种的一些命令。表2.2.1极限的主要命令及说明命令说明Limitexpr,xx0当x趋向于x0时求expr的极限Limitexpr,xx0,Direction1当x趋向于x0时求expr的左极限Limitexpr,xx0,Direction-1当x趋向于x0时求expr的右极限Infinity无穷大趋向的点可以是常数，也可以是+ oo, -oo。注意 Mathemica 没有区分8和+oo,求X 8时的极限要小心。卜

10、面就具体操作几个运行极限的Mathemica 程序。，x2 23x 62.2.1 。2x 21 .求 limx 3x 6利用Limitexpr,x Infinity命令，计算lim expr ;再将表达式expr转化成其中SqrtxA2+2是指Jx2 2 。具体运行程序参见图3x 6图 2.2.1 运彳f lim 2- 的 Mathemica 程序 x2.求2sin x lim 2 x 0 x禾U用 Limitexpr,x0命令，计算lim expr;再将表达式x 02V一.expr转化成,其中xSinxA2是指sin2x。具体运行程序参见图2.2.2 。图2.2.2运行limx 0.2sin

11、 x .2的 Mathemica 程序x3.求lim In xx 0禾U用 Limitexpr,x0,Direction-1命令，计算lim expr ;再将表达式expr转化 x 0成Inx ,其中Logx是指Inx 。具体运行程序参见图2.2.3 。图 2.2.3 运行 lim in x 的 Mathemica 程序x 02.2.2 实验内容2函数的微分在Mathematica中，计算函数的微分或是非常方便的，命令为Df,x,表示对x求函数f的导数或偏导数。该函数的常用格式有以下几种如表2.2.2。表2.2.2关于函数微分的几个主要命令及说明命令说明Df,x计算导数更 dxDf,x,nn

12、n计算n阶导数f.n dx卜面就具体操作几个运行微分的Mathemica程序。例2.2.1 求函数exsinx的导数。利用Df,x命令，计算f的导数；再将表达式f转化成exsinx ,其中Expx是指ex。这样就可得到函数 exsinx导数的Mathemica程序。具体运行程序参见图2.2.4 。图2.2.4运行exsinx的导数的Mathemica程序例2.2.2 求函数exsinx的二阶导数利用Df,x,n 命令，计算f的n导数；再将表达式f换成exsinx ,其中Expx是指ex; n换成2。就可得到函数 exsinx的二阶导数的 Mathemica程序。具体运行程序参见图 2.2.5

13、。图实验2.5函数exsinx的二阶导数的 Mathemica程序例2.2.3假设a是常数，对sinax求导。这题仍可以利用 Df,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的导数。首先利用Df,x命令，将f转化成表达式xg(x)。其运行结果是g冈+xg x。具体运行程

14、序参见图2.2.6 。再求函数xg(x)对x的4阶导数。首先利用 Df,x,n命令，将f转化成表达式xg(x), n换成4。其运行结果是 4g冈+xg冈。具体运行程序参见图2.2.6 。图2.2.6函数xg(x)对x的导数和4阶导数的Mathemica程序对复合函数的求导上面的方法仍适用。如例 2.2.5 。例2.2.5 求函数ghx对x求导。为了求函数ghx对x求导，首先利用 Df,x命令，将f转化成表达式ghx 其运行结果是g hxhx。具体运行程序参见图2.2.7 。图2.2.7函数ghx 对x求导的 Mathemica程序如果要得到函数在某一点的导数值可以把这点代入导数即可。如例2.2

15、.6中求f (2)。例 2.2.6 设 f(x)=e xsinx,求 f (2)。首先，利用Df,x/x 2命令求f (2)。再将f用表达式Expx*sinx 替换，其中Expx表示exo具体运行程序参见图2.2.8 。图 2.2.8 求 f (2) 的 Mathemica 程序2.2.3 实验内容3计算积分1 .不定积分在Mathematica中计算不定积分命令为Integeratef,x,当然也可使用工具栏直接输入不定积分式。来求函数的不定积分。当然并不是所有的不定积分都能求出来。例如求sin(sin x)dxMathematica就无能为力。其运算结果见图2.2.9 。图 2.2.9 积

16、分 sin(sin x)dx Mathematica 运算结果但对于一些手工计算相当复杂的不定积分，MatheMatica还是能轻易求得，例如下面例题 2.2.7 。例2.2. 7求u 1 u2 , 7 du O2 11u若用手工运算相当麻烦，但用MatheMatica程序运算就容易多了。先用工具栏直接输入2.2.11 。不定积分式，在运行结果。如图 2.2.10的MatheMatica程序设计结果。图 2.2.11 积分 sin(sinx)d(sinx)的 Mathematica 运算结果对于在函数中出现的除积分变量外的符号，统统当作常数处理，例如例2.2.9积分(ax2+bx+c)dx中a

17、、b、c在运算过程中就当常数处理的。例 2.2.9 计算积分(ax2+bx+c)dx。首先，输入 (ax 2+bx+c)dx 的 Mathematica 运行程序后，即 可以得到结果2 3cx+ 二+之。具体运行情况参见图2.2.12 。3图 2.2.12 积分（ax2+bx+c）dx 的 Mathematica 运算结果2.定积分定积分的求解主要命令也是用Integrate只是要在命令中加入积分限Integratef,x,min,max ,或者使用工具栏输入也可以。例 2.2.12 定积分 4 x2eaxdx。 0利用命令Integratef,x,min,max ,将f转化成表达式 x2ea

18、x,其中Expax表示eax。128eax其运行结果是 8e。具体运行程序参见图2.2.13 。图2.2.13定积分x2eaxdx的Mathematica运算结果0命令Integratef,x,min,max ,也可以求广义积分。例如下面例2.2.13 。41例2.2.13 求积分 y dx。0 (x 2)2这是一个瑕积分，照例用命令Integratef,x,min,max直接输入即可得到结果。详见图2.2.14 。图2.2.14积分 1dx的Mathematica运算结果0 (x 2)2同样，命令Integratef,x,min,max也可求无穷积分。例如例 2.2.14 。1例2.2.14

19、 求积分rdx °1 x4这是一个广义积分，照例用命令Integratef,x,min,max直接输入即可得到结果。详见图2.2.15。图2.2.15积分工dx的Mathematica运算结果1 x4如果广义积发散也能给出结果，例如例 2.2.15。例2.2.15求积分 -J-dx °1 x2这是一个发散的广义积分，照例用命令Integratef,x,min,max直接输入也即可得到结果。详见图 2.2.16图2.2.16积分4 dx的Mathematica运算结果1x2对于广义积分如果 Mathematics运算程序无法判定敛散性，就用给出一个提示， 例如例2.2.16

20、,其运算结果参见图2.2.17 。2 1例 2.2.16 1dx。 0 x照例用命令Integratef,x,min,max直接输入，因为x=0是奇异点，按一般广义积分的定义，这个积分发散，这是出现提示，返回的只是原输入式的输出形式。具体结 果参见图2.2.17 。图2.2.17积分 Idx的Mathematica运算结果0 x如果广义积分敛散性与某个符号的取值有关，它也能给出在不同情况下的积分结果例如例 2.2.17。例2.2.17 求积分 dxo ixp用命令Integratef,x,min,max直接输入，+用Infinity 代替。输出结果含有一个字母参数p,这时返回的结果是一个条件表

21、达式，当 p的实部大于1时值为 ,否则发散。具体结果参见图2.2.18 。图2.2.18 积分dx的Mathematica运算结果1 xp在 Integrate中可加两个参数 Assumptions 和 GenerateConditions 例如例 2.2.17 中,只要用Assumptions Rep>1就可以得到收敛情况的解。图2.2.19积分 dx在p>1的情况下的 Mathematica运算结果1 xp3.数值积分数值积分是解决求定积分的另一种有效的方法，它可以给出一个近似解。特别是对于用Integrate 命令无法求出的定积分，数值积分更是可以发挥巨大作用。它的主要命令格式为表2.2.3中的命令。表2.2.3数值积分一些主要命令11Nintegratef,x,a,b在a,b上求f数值积分Nintegratef,x,a,x1,x2,b以x1,x2.为分割求a,b上的数值积分Nintegratef,x,a,b,MaxRecursionn求数值积分时指定迭代次数n.下面我们求sinsinx 在0,兀上的积分值，由于这个函数的不定积分求不出，因此使用Integrate命令无法得到具体结果，但可以用数彳1积分可以求得近似值。如例 2.2.18。例 2.2.18 求积分 sin(sinx)dx 。这个被积函数f(x)=sin(sinx)