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

第二篇 数学试验第 2 章 数学试验2.2 实验 2一元微积分的编程实现1 编程实现数值积分的几种编程实现数值积分的几种-- --方法方法 c c 语言语言 数值计算 2010-11-05 09:52:43 阅读 385 评论 1字号：大中小 订阅 复化梯形公式 在区间不大时 , 用梯形公式、 辛卜生公式计算定积分是简单实用的 , 但当区间较大时 , 用梯形公 式、辛卜生公式计算定积分达不到精确度要求 . 为了提高计算的精确度，我们将 [a,b] 区间 n 等分，在每个小区间上应用梯形 公式、辛卜生公式计算定积分，然后将其结果相加，这样就得到了复化梯形公式和复化辛卜生公式。 1. 复化梯形公式 将积分区间等分 , 设, 则节点为 对每个小区间上应用梯形公式 , 然后将其结果相加，则得 (3.14) 称 (3.14) 式为复化梯形公式 . 当在 [a,b] 上有连续的二阶导数时，则复化梯形公式 (3.14) 的余项推导如下： 第二篇 数学试验第 2 章 数学试验2.2 实验 2一元微积分的编程实现2 因为 所以在区间 [a,b] 上公式 (3.14) 的误差为 又因为在区间 [a,b] 上连续，由连续函数的性质知，在区间 [a,b] 上存在一点， 于是 （ 3.15 ） 复化梯形公式，复化抛物线公式和 Romberg 求积法的算法程序： 以下程序均定义误差限为 1*10^-5; 1)复化梯形公式： #include #include 第二篇 数学试验第 2 章 数学试验2.2 实验 2一元微积分的编程实现3 #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;ie);//自定义误差限 e printf(“%.8lf“,t);//输出积分的近似值 return 0; } 2)复化抛物线公式： #include #include #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)); 第二篇 数学试验第 2 章 数学试验2.2 实验 2一元微积分的编程实现5 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;ie);//自定义误差限 printf(“%.8lf“,s); return 0; } 第二篇 数学试验第 2 章 数学试验2.2 实验 2一元微积分的编程实现6 3)Romberg 求积法： #include #include #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;i1 的情况下的 Mathematica 运算结果 3.3.数值积分数值积分 第二篇 数学试验第 2 章 数学试验2.2 实验 2一元微积分的编程实现25 数值积分是解决求定积分的另一种有效的方法， 它可以给出一个近似解。 特别是对于 用 Integrate 命令无法求出的定积分，数值积分更是可以发挥巨大作用。 它的主要命令格 式为表 2.2.3 中的命令。 表 2.2.3 数值积分一些主要命令 Nintegrate[f,{x,a,b}]在[a,b]上求 f 数值积分 Nintegrate[f,{x,a,x1,x2,…,b}]以 x1,x2….为分割求[a,b]上的数值积分 Nintegrate[f,{x,a,b},MaxRecursionn]求数值积分时指定迭代次数n. 下面我们求 sinsinx 在[0,π]上的积分值，由于这个函数的不定积分求不出，因此使用 Integrate 命令无法得到具体结果，但可以用数值积分可以求得近似值。如例2.2.18。 例 2.2.18 求积分  0 sin(sinx)dx。 这个被积函数f(x)=sin(sinx)找不到初等的原函数， 用命令 Integrate[f,{x,min,max}] 直接输入，得不到结果。但用命令Nintegrate[f,{x,a,b}]输入，即可得到近似值。具体情 况参见图 2.2.20，其中 Pi 表示π。 第二篇 数学试验第 2 章 数学试验2.2 实验 2一元微积分的编程实现26  图 2.2.20 积分 0 sin(sinx)dx 在 Mathematica 程序下的数值解 如果积分的被函数存在不连续点，或存在奇异点我们可对积分进行分段求解。例如例 2.2.19。 例 2.2.19求积分1 1 x1 dx。 被积函数 1 x 在[-1，1]上，显然有x=0 点是一个无穷间断点。因此若要求其数值积分， 必须在其中插入点 0。 于是， 利用数值积分命令 Nintegrate[f,{x,xmin， x1， x2， ……， xmax}] 输入，其中 x1，x2，……是奇异点。这样就可得到积分的数值解。具体操作参见图2.2.21 第二篇 数学试验第 2 章 数学试验2.2 实验 2一元微积分的编程实现27 图 2.2.21 积分 11 x1 dx 在 Mathematica 程序下的数值解 对无穷积分，也可求得积分数值解，例如例2.2.20。 例 2.2.20 求积分 0 exdx。 2 对于求无穷积分的数值解，也可以直接用命令Nintegrate[f,{x,a,b}]，其中b 可以 用无穷大 Infinity 替换。具体操作参见图2.2.22。 第二篇 数学试验第 2 章 数学试验2.2 实验 2一元微积分的编程实现28 图 2.2.22积分  0 exdx在 Mathematica 程序下的数值解 2 思考题： 可以将第一篇中的习题一、习题二中的题用Mathematica 程序做一下。