文档视界 最新最全的文档下载
当前位置:文档视界 › 数值分析变步长求解积分

数值分析变步长求解积分

数值分析变步长求解积分
数值分析变步长求解积分

#include

#include

# define N 100

double function(double s)

{

double result;

result=sin(s)/s;

return result;

}

int explode2(int k)

{

int i,s=1;

for(i=1;i

s=s*2;

return s;

}

int explode4(int k)

{

int i,s=1;

for(i=1;i

s=s*4;

return s;

}

void main()

{

int i,j,k=1,m=1;

double sum=0,sum1=0,sum2=0,sum3=0,d[4],a=1,b=2,f[N+1],h,T[4][4],c[4]={1,1,1,1};

h=(b-a)/N;

for(i=0;i

{

f[i]=function(a+i*h);

}

for(i=0;i

{

sum=sum+h*(f[i]+f[i+1])/2;

}

T[0][0]=sum;

for(i=0;i<2*N-1;i++)

{

sum1+=function(a+(i+i/2)*h/2);

}

for(i=0;i<2*2*N-1;i++)

{

sum2+=function(a+(i+i/2)*h/4);

}

for(i=0;i<2*2*2*N-1;i++)

{

sum3+=function(a+(i+i/2)*h/8);

}

d[0]=sum;d[1]=sum1;d[2]=sum2;d[3]=sum3;

for(j=1;j<4;j++)

T[i][j]=1/2*T[i][j-1]+h/explode2(j)*d[j];

for(i=1;i<4;i++)

{

for(j=0;j<4-i;j++)

T[i][j]=explode4(i)/(explode4(i)-1)*T[i-1][j+1]-1/explode4(i)*T[i-1][j];

}

printf("Romberg=%lf",T[3][0]);

}

数值分析第四章数值积分与数值微分习题答案

第四章 数值积分与数值微分 1.确定下列求积公式中的特定参数,使其代数精度尽量高,并指明所构造出的求积公式所具有的代数精度: 101210121 12120 (1)()()(0)(); (2)()()(0)(); (3)()[(1)2()3()]/3; (4)()[(0)()]/2[(0)()]; h h h h h f x dx A f h A f A f h f x dx A f h A f A f h f x dx f f x f x f x dx h f f h ah f f h -----≈-++≈-++≈-++''≈++-?? ?? 解: 求解求积公式的代数精度时,应根据代数精度的定义,即求积公式对于次数不超过m 的多项式均能准确地成立,但对于m+1次多项式就不准确成立,进行验证性求解。 (1)若101(1) ()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =,则 1012h A A A -=++ 令()f x x =,则 110A h Ah -=-+ 令2 ()f x x =,则 3 221123 h h A h A -=+ 从而解得 011431313A h A h A h -?=?? ? =?? ?=?? 令3 ()f x x =,则 3()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++=

令4()f x x =,则 455 1012()5 2 ()(0)()3 h h h h f x dx x dx h A f h A f A f h h ---== -++=? ? 故此时, 101()()(0)()h h f x dx A f h A f A f h --≠-++? 故 101()()(0)()h h f x dx A f h A f A f h --≈-++? 具有3次代数精度。 (2)若 21012()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =,则 1014h A A A -=++ 令()f x x =,则 110A h Ah -=-+ 令2 ()f x x =,则 3 2211163 h h A h A -=+ 从而解得 1143 8383A h A h A h -?=-?? ? =?? ?=?? 令3 ()f x x =,则 22322()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++=

复化梯形法 复化矩形法 变步长梯形 变步长辛普森

陕西科技大学 机械教改班 用C++的积分 其实积分的思想就是,微分—>求和—>取极限,如果是用纯手工法那就是先对一个函数微分,再求出它的面积,在取极限,因为我们的计算速度和计算量有限,现在有了计算机这个速度很快的机器,我们可以把微分后的每个小的面积加起来,为了满足精度,我们可以加大分区,即使实现不了微分出无限小的极限情况,我们也至少可以用有限次去接近他,下面我分析了四种不同的积分方法,和一个综合通用程序。 一.积分的基本思想 1、思路:微分—>求和—>取极限。 2、Newton —Leibniz 公式 ?-=b a a F b F dx x f ) ()()( 其中,)(x F 被积函数)(x f 的原函数。 3、用计算机积分的思路 在积分区间内“微分—>求和—>控制精度”。因为计算机求和不可以取极限,也就是不可以无限次的加下去,所以要控制精度。 二.现有的理论 1、一阶求积公式---梯形公式 ?=+-=b a T b f a f a b dx x f )]()([2 )( 他只能精确计算被积函数为0、1次多项式时的积分。 2、二阶求积分公式——牛顿、科特斯公式 ?=+++-=b a S b f a b f a f a b dx x f )]()2(4)([6)( 他只能精确计算被积函数为0、1、2、3次多项式时的积分。 三.四种实现方法 1.复化矩形法 将积分区间[a,b]等分成n 个子区间: ],[],[],[],[],[112322110n n n n x x x x x x x x x x ---、、、 则h=(b-a)/n,区间端点值k x =a+kh

DEFORM步长控制

DEFORM模拟控制(三):步进控制 DEFORM的步进控制是用来控制DEFORM计算每一步的递增量的,这个递增量可以是模具的位移,也可以是时间,甚至是温度。其分别对应的意义为每步模具走了多大距离,每步表示了多长时间,每步变化了多少温度。 如果指定了每步行程,则主模(一般为上模)将在每个步进中移动指定的量。主模具的总运动是每步位移乘以步总数。如果指定了每步时间,则将使用每步时间间隔。每步模具位移将是时间步乘以模具速度。 从DEFORM 3D V6.1开始,步进增量控制的定义已得到增强,以包括与时间和行程有关的步进功能。这意味着步长(每步时间和每步行程)现在可以定义为时间或冲程的函数。此功能可以在需要的地方更好地解析保存的模型信息。 定义每步行程通常更直观。但是,对于没有模具移动(例如传热)的任何问题或使用力控制的任何问题,必须指定每步时间。 1 Solution step definition(步进模式) 这里可以选择定义步进的模式,分别是模具行程(Die displacement)、时间(Time)、温度(Temperature)。温度选项不常用。 2 Step increment control(步长控制) 这里分情况讨论,当步进模式选择的是行程控制的时候,太大的步长会导致大的结果误差、网格快速变形或者难以收敛等问题。而太小的步长会浪费时间。所以,通常情况下,行程控制的步长应该设置成不超过最小网格单元边长的1/3。对于狭窄拐角的流动,飞边成形或类似的高度局部变形,可能需要控制步长小至网格单元边长的1/10 当步进模式选择的是时间控制的时候,通常按以下步骤进行计算设置: 1 使用测量工具,测量变形体的最小单元边长(划分完网格之后测量) 2 估计变形体所有区域的最大变形速度(大多数情况下,这个值就是模具的速度,对于挤压问题来说,这个值是模具的速度乘以挤压比)。假如你的模拟已经进行了一个阶段,则可以通过查看节点数据查看。

变步长的梯形积分方法的应用

CENTRAL SOUTH UNIVERSITY 数值分析实验报告

变步长的梯形积分方法的应用 一、问题背景 实际问题当中常常需要计算积分,有些数值方法,如微分方程和积分方程的求解,也都和积分计算相关联。 依据人们所熟知的微积分基本定理,对于积分 ()dx x f I a ?=b , 只要找到被积分函数()x f 的原函数()x F ,便有下列牛顿-莱布尼茨(Newton-Leibniz )公式: ()()()a F b F dx x f b -=?a . 但实际使用这种求积分方法往往有困难,因为大量的被积函数,诸如 ()0sin ≠x x x ,2x e -等,其原函数不能用初等函数表达,故不能用上述公式计算。即使能求得原函数的积分,有时计算也十分困难。例如对于被积函数 ()6 11x x f +=,其原函数 ()C x x x x x x x x F ++-+++??? ??-+=1 313ln 3411arctan 61arctan 3122, 计算()a F ,()b F 仍然很困难,另外,当()x f 是由测量或数值计算给出的一张数据表时,牛顿-莱布尼茨公式也不能直接运用。因此有必要研究积分的数值计算问题。 二、数学模型 由于牛顿-科特斯积分公式在8≥n 时不具有稳定性,故不能通过提高阶数的方法来提高求积精度。为了提高精度通常可以把积分区间划分成若干的子区间(通常是等分),再在每个子区间上用低阶求积公式。这种方法称为复合求积法。 复合梯形法虽然方法简单,但是却不能估计积分精度,这有时候是很不方便的。要想控制积分精度,可以采用如下的方法,设积分区间已经划分为n 个子区间,这时再把区间划分更细,给出新的积分结果,如果前后两次积分的差比给定的误差容限小的话,则停止细华否则继续增加积分区间。这种方法原理很简单也 容易实现,但是实际计算中一般采用的比较少,因为这种方法比较机械效率不是太高,实际上采用比较多的通常是Romberg 方法。 三、算法及流程 给定义误差容限小量TOL ,对于()dx x f b a ?,有复合梯形公式

FLUENT控制步长时间courant数的有效的经验

" What is the difference between the time accurate solution of Navier-Stokes equations and the DNS solution? " I'm not sure I understand this question correctly, but as far as I am aware DNS (Direct Numerical Simulation) is defined as a time accurate solution of the Navier-Stokes equations. Perhaps you mean "What is the difference between laminar & turbulent DNS?" ? Believing this to be so from reading the rest of your message I wrote the following: When the DNS is of turbulence rather than a laminar flow the turbulence requires initialization in some way. The same is true of LES (LES is always turbulent, because laminar LES = DNS by definition). I have found little reported work on the proceedures used to accomplish this turbulence initialization. I can only speak for the LES code I use (and the generations of code that preceded it). The turbulence is initialized by setting an initial flow field that has a random fluctation velocity component added to the inital mean velocity. ----------------- For example: U_initial_cell = U_initial_mean + (Random_number * U_initial_mean * 0.20) ( Random_number has a value in the range -1 to +1 ) This function sets the initial cell velocity to that of the initial_mean with a tolerance of 20% (i.e. + or - 20%). So if the U_initial_mean was 1.0 then the initial velocity of the cell could be anywhere between 0.8 and 1.2 depending on the Random_number (an intrinsic computer function). ----------------- The value for initial fluctuation (20% etc) is a fairly arbitrary value just required to `kick-start' the turbulence. Once the simulation has been kick-started and run sufficiently long enough for the correct energy cascade to be observed (by monitoring k.e. of the flow) the statistics data from that point onward is o.k to be used for results. This accumulation of statistical data is one of the reasons why LES/DNS turbulent simulations require so much more time to run. Providing the same random_number is used on the same cells during initialization the computations of two DNS cases will be exactly* the same when all other conditions (boundary, geometry, etc) are equivalent. *exactly is defined as Phi(x,t)_simulation_A = Phi(x,t)_simulation_B

数值分析常微分数值解的求法C语言

本科生课程设计报告实习课程数值分析 学院名称管理科学学院 专业名称信息与计算科学 学生姓名 学生学号 指导教师 实验地点 实验成绩 二〇一六年六月二〇一六年六 摘要 常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题,实用上许多很有价值的常微分方程的解不能用初等函数来表示,常常需要求其数值解.所谓数值解,是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.

关键词:数值解法;常微分

1. 实验内容和要求 常微分方程初值问题 有精确解2()cos(2)x y x x e x -=+。 要求:分别取步长h = ,,,采用改进的Euler 方法、4阶经典龙格-库塔R -K 方法和4阶Adams 预测-校正方法计算初值问题。以表格形式列出10个等距节点上的计算值和精确值,并比较他们的计算精确度。其中多步法需要的初值由经典R-K 法提供。运算时,取足以表示计算精度的有效位。 2. 算法说明 函数及变量说明 表1 函数及变量定义 1、欧拉方法: 1()()(,())i i i i y x y x hf x y x +=+ (i=0,1,2,3,......n-1) (0)y a = (其中a 为初值)

2、改进欧拉方法: ~ 1~111()()(,()) ()()[(,())(,())] 2 (0)i i i i i i i i i i y x y x hf x y x h y x y x f x y x f x y x y a ++++=+=++=(i=0,1,2......n-1) (其中a 为初值) 3、经典K-R 方法: 1 1213 243()6 (,)(,)22(,)22 (,) i i i i i i i i i i h y y K f x y h h K f x y K h h K f x y K K f x h y hK +? =+?? =???=++?? ? =++??=++?? 4、4阶adams 预测-校正方法 预测: 校正: Adsms 内插外插公式联合使用称为Adams 预测-校正系统,利用外插公式计算预测,用内插公式进行校正。计算时需要注意以下两点: 1、外插公式为显式,内插公式为隐式。故用内插外插公式时需要进行迭代。 2、这种预测-校正法是四步法,计算Yn+1时,不但用到前一步信息,而且要用到更前三步信息1-n f ,2-n f 3-n f ,因此它不是自动开始的,实际计算时必须借助某种单步法,用Runge-Kutta 格式为它提供初始值 3 21,,y y y ,依据局部截断误差公式得:

变步长梯形法

#include #include float f(float x) { float s; s=1/(x*x); return(s); } main() { float a,b,c,h,x,T1,T2,S; float T0; printf("变步长梯形法求积分:\n"); printf("需要求解的积分式为f(x)=1/(x*x)\n"); printf("请输入a: "); scanf("%f",&a); printf("请输入b: "); scanf("%f",&b); printf("请输入c: "); scanf("%f",&c); h=b-a; T1=h*(f(a)+f(b))/2; S=0; x=a+h/2; do { S=S+f(x); x=x+h; }while(x=c) { T1=T0; S=0; x=a+h/2; do { S=S+f(x); x=x+h; }while(x

T2=T1/2+h/2*S; printf("步长为%f时\t Tn=%f \t T2n=%f \t 误差变化:%f\n",h,T1,T2,fabs(T1-T2)); h=h/2;T0=T2; } printf("************************************************************** **************\n"); printf("所求的结果为=%f\n",T2); }

光伏发电系统新型变步长MPPT控制方法

- 178 - 光伏发电系统新型变步长MPPT 控制方法 范佳佳 宋平岗 (华东交通大学电气与电子工程学院,江西 南昌 330013) 【摘 要】分析并网型光伏发电系统最大功率点跟踪原理,针对含有直流-直流(DC/DC)升压环节的两级系统,文章提出了一种新型简单实用变步长最大功率点跟踪(MPPT)控制方法。系统只检测输出并网电流,并根据电流变化率di/dt 变步长改变占空比,实现最大功率点的快速跟踪。仿真分析得,系统追踪最大功率点速度快,当外界条件发生变化时,系统能快速追踪到新的最大功率点并保持不变,稳定性好,输出功率平稳,波动小,证明了理论研究的正确性与可行性。 【关键词】并网;光伏发电;MPPT;变步长 【中图分类号】TM464 【文献标识码】B 【文章编号】1008-1151(2010)10-0178-03 (一)引言 太阳能具有清洁、无污染、取之不尽用之不竭、分布面广等优点,已成为世界各国普遍研究的热点。与此同时,光伏发电装机总量每年成倍增趋势,在世界发电总量中所占的比重也越来越大。但是,光伏发电系统的一个主要缺点是其价格昂贵而发电效率却较低,这已成为限制其发展的主要因素之一[1]。所有的光伏系统都希望电池阵列在同样的日照、温度条件下输出尽可能多的电能,提高发电效率,这也就是理 论和实践上的太阳能电池阵列最大功率点跟踪MPPT(Maximum Power Point Tracking)控制。 目前,国内外常用的最大功率点跟踪控制方法大多是直 接采样太阳能电池的输出电压和输出电流,计算输出功率并改变占空比,比较功率变化趋向。比较成熟的MPPT 方法有扰 动观察法(爬山法)、增量电导法(INC)、恒定电压法(CVT)、滞环比较法、模糊逻辑控制法等。其中工程中应用较多的控制方法为:扰动观察法和恒定电压法。此类方法存在一个缺点,系统在追踪最大功率点时只能定步长跟踪,追踪速度较慢,会产生一定的能量流失。另外,此类最大功率点跟踪控 制方法一般直接采样太阳能电池的输出电压和输出电流,计 算输出功率并改变占空比,比较功率变化趋向[2]。 针对上述缺点,本文采用一种直接电流控制法,系统只 采样逆变器输出并网电流,根据并网电流的大小判断功率趋 向并改变占空比。与传统的直接电流控制法不同的是,系统 会根据采样得到的并网电流的变化率改变占空比的变化率, 实现系统最大功率点的快速跟踪。由于只对输出并网电流进 行采样,系统在不降低性能的情况下,节省了一个传感器,成本得以降低,软硬件设计也得以简化。 (二)太阳能光伏电池模型及输出特性 光伏电池是一种能吸收太阳能并将其转换为电能的半导体装置,由多个p-n 结串联而成,它的输出具有非线性特征,受环境因素影响很大,如光照强度和环境温度。随着光照强度和温度的变化,光伏电池的输出电压也会随之产生较大变化。当光伏电池处于暗处时,其输出的伏安特性与二极管指数特性相似。当光伏电池受一定强度的太阳光照射时,p-n 结两边会产生电势差,将电势差引出外导体产生电流,相当于电流源输出。太阳能电池板开路时,电流被其自身固有的p-n 结二极管分流,此时这个二极管的特性可被看成是光伏电 池的开路电压特性[3-4] 。根据光伏电池的内部结构和输出伏安 特性得到光伏电池的等效电路如图1所示,它由一个受光强和温度影响的电流源并联一个二极管再串联一个电阻组成。 图1 光伏阵列的等效电路模型 其等效数学模型为: ()()sh s s ph R I R U q nkT I R U q I I I +??????????????+?=1exp 0 (1) 式中:I ——光伏电池输出电流(工作电流);U ——光伏电池输出电压(工作电压);I ph ——光生电流;I 0——二极管饱和电流;q ——电子的电荷量(1.6×1910?C);R s ——光伏电池的串联电阻; n ——二极管特性因子;k ——玻尔兹曼常数(1.38×2310?J/K);T ——光伏电池温度;R sh ——光伏电池 的并联电阻。 基于式(1),我们使用Matlab 搭建了光伏电池的仿真 模型,如图2所示。仿真采用无锡尚德太阳能电力有限公司生产的单晶硅光伏电池组件(STP175S-24/Ab),参数为:U oc =44.2V,U m =35.2V,I sc =5.2A,I m =4.95A,T ref =25℃,R S =0.1Ω。 图2 PV 模块 (三)新型变步长MPPT 算法 大多并网光伏发电系统,变流器输出端经过滤波环节后与电网相连,通常电网可视为无穷大电压源系统,其特性较为稳定[5],如电网电压与频率等参数,当系统并入电网时,系统输出端并网电压被强制箝位在电网电压上。由于变换器自 【收稿日期】2010-07-09 【作者简介】范佳佳(1987-),男,安徽芜湖人,华东交通大学电气与电子工程学院硕士研究生,研究方向为光伏变流器拓扑结构及控制方法;宋平岗,男,华东交通大学电气与电子工程学院硕士研究生导师,研究方向为新能源及大功率变流技术。

simulink步长设置

simulink仿真设置 一、算法设置 1.变步长(Variable—Step)求解器 可以选择的变步长求解器有:ode45,ode23,ode113,odel5s,ode23s 和discret.缺省情况下,具有状态的系统用的是ode45;没有状态的系统用的是discrete。 1)ode45基于显式Runge—Kutta(4,5)公式,Dormand—Prince对.它是—个单步求解器(solver)。也就是说它在计算y(tn)时,仅仅利用前一步的计算结果 y(tn-1).对于大多数问题.在第一次仿真时、可用ode45试一下。 2)ode23是基于显式Runge—Kutta(2,3).Bogackt和Shampine对.对于宽误差容限和存在轻微刚性的系统、它比ode45更有效一些.ode23也是单步求解器。 3)odell3是变阶Adams-Bashforth—Moulton PECE求解器.在误差容限比较严时,它比ode45更有效.odell3是一个多步求解器,即为了计算当前的结果y(tn),不仅要知道前一步结果y(tn-1),还要知道前几步的结果y(tn-2),y(tn-3),…; 4)odel5s是基于数值微分公式(NDFs)的变阶求解器.它与后向微分公式BDFs(也叫Gear方法)有联系.但比它更有效.ode15s是一个多步求解器,如果认为一个问题是刚性的,或者在用ode45s时仿真失败或不够有效时,可以试试odel5s。odel5s是基于一到五阶的NDF公式的求解器.尽管公式的阶数越高结果越精确,但稳定性会差一些.如果模型是刚性的,并且要求有比较好的稳定性,应将最大的阶数减小到2.选择odel5s求解器时,对话框中会显示这一参数.可以用ode23求解器代替。del5s,ode23是定步长、低阶求解器。 5)ode23s是基于一个2阶改进的Rosenbrock公式.因为它是一个单步求解器,所以对于宽误差容限,它比odel5s更有效.对于一些用odel5s不是很有效的刚性问题,可以用它解决。 6)ode23t是使用“自由”内插式梯形规则来实现的.如果问题是适度刚性,而且需要没有数字阻尼的结果,可采用该求解器。 7)ode23tb是使用TR—BDF2来实现的,即基于隐式Runge—Kutta公式,其第一级是梯形规则步长和第二级是二阶反向微分公式.两级计算使用相同的迭代矩阵.与ode23s相似,对于宽误差容限,它比odtl5s更有效。 8)discrete(变步长)是simulink在检测到模型中没有连续状态时所选择的一种求解器。

变步长梯形公式(C语言)

变步长梯形公式(C语言) 程序: // cehngxu.cpp : 定义控制台应用程序的入口点。 // #include"stdafx.h" #include"stdio.h" #include"math.h" double f(double x) //这里自定义了函数,因为出现sin(0)/0的情况,系统无法计算;{ double y; if (x == 0) y = 1; //把x=0作为一种情况,单独拿出来; else y = sin(x) / x; //正常情况下的函数; return(y); } //编写这个自定义函数便于你在作业中的计算,对于不用的题目只需改动一下函数即可计算; void main() //主函数; { double a , b ,h,k,s,t[4998]; //数组的定义中不能出现变量,故对其任意取值,一般取一个很大的数; int n = 1, m = 0; printf("please input a="); scanf_s("%lf", &a); //在VC环境下用scanf输入没有问题,但是在我这编译环境visual studio 2013下需要用scanf_s()输入; printf("please input b="); scanf_s("%lf", &b); t[1] = (b - a)*(f(a) + f(b)) / 2; do { h = (b - a) / n; s = 0; for (k = 0; k < n; k++) { s += f(a + (k + 0.5)*h); } t[2 * n] = t[n] / 2 + h / 2 * s; //梯形的递推公式 n = 2*n; m = m + 1;

数值分析变步长求解积分

#include #include # define N 100 double function(double s) { double result; result=sin(s)/s; return result; } int explode2(int k) { int i,s=1; for(i=1;i

} for(i=0;i<2*2*2*N-1;i++) { sum3+=function(a+(i+i/2)*h/8); } d[0]=sum;d[1]=sum1;d[2]=sum2;d[3]=sum3; for(j=1;j<4;j++) T[i][j]=1/2*T[i][j-1]+h/explode2(j)*d[j]; for(i=1;i<4;i++) { for(j=0;j<4-i;j++) T[i][j]=explode4(i)/(explode4(i)-1)*T[i-1][j+1]-1/explode4(i)*T[i-1][j]; } printf("Romberg=%lf",T[3][0]); }

变步长,辛普森,梯形公式

#include #include #include #include #include #define ESC 27 #define EPS 0.5e-7 #define f(x1) (-2.0/(x1*x1-1)) #define g(x2) (-2.0/(x2*x2-1)) #define s(x0) (-2.0/(x0*x0-1)) char t; void tixing() { double a1=2,b1=3; double T,h1,x1; int n,i; printf("please input n:"); scanf("%d",&n); h1=(b1-a1)/n; x1=a1;T=0; for(i=1;i

数值分析第4章数值积分与数值微分

第4章 数值积分与数值微分 1 数值积分的基本概念 实际问题当中常常需要计算定积分。在微积分中,我们熟知,牛顿—莱布尼兹公式是计算定积分的一种有效工具,在理论和实际计算上有很大作用。对定积分()b a I f x dx =?,若()f x 在区间 [,]a b 上连续,且()f x 的原函数为()F x ,则可计算定积分 ()()()b a f x dx F b F a =-? 似乎问题已经解决,其实不然。如 1)()f x 是由测量或数值计算给出数据表时,Newton-Leibnitz 公式无法应用。 2)许多形式上很简单的函数,例如 2 22sin 1(),sin ,cos ,,ln x x f x x x e x x -= 等等,它们的原函数不能用初等函数的有限形式表示。 3)即使有些被积函数的原函数能通过初等函数的有限形式表示,但应用牛顿—莱布尼兹公式计算,仍涉及大量的数值计算,还不如应用数值积分的方法来得方便,既节省工作量,又满足精度的要求。例如下列积分 2 41arc 1)arc 1)1dx tg tg C x ? ?=+++-+??+? 对于上述这些情况,都要求建立定积分的近似计算方法——数值积分法。 1.1 数值求积分的基本思想 根据以上所述,数值求积公式应该避免用原函数表示,而由被积函数的值决定。由积分中值定理:对()[,]f x C a b ∈,存在[,]a b ξ∈,有 ()()()b a f x dx b a f ξ=-? 表明,定积分所表示的曲边梯形的面积等于底为b a -而高为()f ξ的矩形面积(图4-1)。问题在于点ξ的具体位置一般是不知道的,因而难以准确算出()f ξ。我们将()f ξ称为区间[,]a b 上的平均高度。这样,只要对平均高度()f ξ提供一种算法,相应地便获得一种数值求积分方法。 如果我们用两端的算术平均作为平均高度()f ξ的近似值,这样导出的求积公式 [()()]2 b a T f a f b -= + (4-1) 便是我们所熟悉的梯形公式(图4-2)。而如果改用区间中点2 a b c +=的“高度”()f c 近似地取代 平均高度()f ξ,则可导出所谓中矩形公式(简称矩形公式) ()2a b R b a f +?? =- ??? (4-2)

相关文档