整理Matlab积分.docx
《整理Matlab积分.docx》由会员分享，可在线阅读，更多相关《整理Matlab积分.docx（18页珍藏版）》请在冰豆网上搜索。
整理Matlab积分
一.数值积分的实现方法
1．变步长辛普生法
基于变步长辛普生法，MATLAB给出了quad函数来求定积分。
该函数的调用格式为：
[I,n]=quad（'fname',a,b,tol,trace）
其中fname是被积函数名。
a和b分别是定积分的下限和上限。
tol用来控制积分精度，缺省时取tol=0.001。
trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。
返回参数I即定积分值，n为被积函数的调用次数。
例8-1求定积分。
（1）建立被积函数文件fesin.m。
functionf=fesin（x）
f=exp（-0.5*x）.*sin（x+pi/6）;
（2）调用数值积分函数quad求定积分。
[S,n]=quad（'fesin',0,3*pi）
S=0.9008
n=77
2．牛顿－柯特斯法
基于牛顿－柯特斯法，MATLAB给出了quad8函数来求定积分。
该函数的调用格式为：
[I,n]=quad8（'fname',a,b,tol,trace）
其中参数的含义和quad函数相似，只是tol的缺省值取10-6。
该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。
（1）被积函数文件fx.m。
functionf=fx（x）
f=x.*sin（x）./（1+cos（x）.*cos（x））;
（2）调用函数quad8求定积分。
I=quad8（'fx',0,pi）
I=2.4674
分别用quad函数和quad8函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。
调用函数quad求定积分：
formatlong;
fx=inline（'exp（-x）'）;
[I,n]=quad（fx,1,2.5,1e-10）
I=0.28579444254766
n=65
调用函数quad8求定积分：
formatlong;
fx=inline（'exp（-x）'）;
[I,n]=quad8（fx,1,2.5,1e-10）
I=0.28579444254754
n=33
3．被积函数由一个表格定义
在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz（X,Y）函数。
其中向量X,Y定义函数关系Y=f（X）。
用trapz函数计算定积分。
命令如下：
x=1:
0.01:
2.5;
Y=exp（-X）;%生成函数关系数据向量
trapz（X,Y）
ans=0.28579682416393
8.1.3二重定积分的数值求解
使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解。
该函数的调用格式为：
I=dblquad（f,a,b,c,d,tol,trace）
该函数求f（x,y）在[a,b]×[c,d]区域上的二重定积分。
参数tol，trace的用法与函数quad完全相同。
计算二重定积分
（1）建立一个函数文件fxy.m：
functionf=fxy（x,y）
globalki;
ki=ki+1;%ki用于统计被积函数的调用次数
f=exp（-x.^2/2）.*sin（x.^2+y）;
（2）调用dblquad函数求解。
globalki;ki=0;
I=dblquad（'fxy',-2,2,-1,1）
ki
I=1.57449318974494
ki=1038
二.数值微分
数值微分的实现
在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为：
DX=diff（X）：
计算向量X的向前差分，DX（i）=X（i+1）-X（i），i=1,2,…,n-1。
DX=diff（X,n）：
计算X的n阶向前差分。
例如，diff（X,2）=diff（diff（X））。
DX=diff（A,n,dim）：
计算矩阵A的n阶差分，dim=1时（缺省状态），按列计算差分；dim=2，按行计算差分。
生成以向量V=[1,2,3,4,5,6]为基础的范得蒙矩阵，按列进行差分运算。
命令如下：
V=vander（1:
6）
DV=diff（V）%计算V的一阶差分
例8-7用不同的方法求函数f（x）的数值导数，并在同一个坐标系中做出f'（x）的图像。
程序如下：
f=inline（'sqrt（x.^3+2*x.^2-x+12）+（x+5）.^（1/6）+5*x+2'）;
g=inline（'（3*x.^2+4*x-1）./sqrt（x.^3+2*x.^2-x+12）/2+1/6./（x+5）.^（5/6）+5'）;
x=-3:
0.01:
3;
p=polyfit（x,f（x）,5）;%用5次多项式p拟合f（x）
dp=polyder（p）;%对拟合多项式p求导数dp
dpx=polyval（dp,x）;%求dp在假设点的函数值
dx=diff（f（[x,3.01]））/0.01;%直接对f（x）求数值导数
gx=g（x）;%求函数f的导函数g在假设点的导数
plot（x,dpx,x,dx,'.',x,gx,'-'）;%作图
Matlab数值积分的一些做法
Filedunder:
未分类—franz@9:
31pm
今天是元宵节，突然来讲Matlab的确很奇怪，连我自己也这样感觉．
由于Matlab对于计算过程的细节要求定义详细，往往让人觉得使用不方便．与同样很流行的Mathematic相比，Matlab更象是程序，而不是象Mathematic那么直观的写作业．
Matlab同样提供符号积分（不定积分）和数值积分（定积分）两大功能．符号积分使用int命令，配合之前的函数定义使用．比如：
a=…;
b=…;
functionf=fun（x）
f=a*x^2+b;
将此文件寸为一个单独m文件，再在主程序中调用即可：
F=int（fun,c,d）;
c和d为积分上下限，通常为符号，可使用symsc;语句定义。
在完成符号积分之后可以对其附值，则就完成了数值积分的任务．
有一点需要注意的是，通过这样过程求的数值，其格式为符号格式’sym’，不能做图，不能和数值进行运算．处理方法是：
vpa（F）; 求得３２位符号近似解，再用double（vpa（F））;将其转换成Matlab默认的双精度位数就可以．注意，直接做double（F）行不通，Matlab会提示你＂Conversionfromsymtodoubleisnotpossible＂，不知道＂notpossbile"和"impossible"有什么区别，呵呵．
符号积分中一个最大的问题在于存在不可积的函数，比如十分常用的Boltzman分布，或者叫Arrhenius公式：
exp（-q*V/（k*T））．
插句话，我一直以为Arrhenius是德国人，知道前几天在和老师讨论中讲起，我老师很自豪的对我说，有个著名的瑞典化学家，得过NoblePrize，叫Arrhenius你知道不知道，才发现这个小问题，３滴汗，搞的我老师也很有挫败感．．．
对于不可积的函数，使用int命令之后，得到的表达式中会含有 Ei 项，其本身也是一个函数，以你给定的符号积分上下限为变量．在含有 Ei 项后，对符号上下限附值亦无法得到数值积分解，因为Ei 项也需要解．解Ei 项的方法如下：
result = str2num （maple（‘evalf（Ei（a,b））’））;
％此语句可以直接使用，其中的a和b就是你所要解数值解的积分上下限，并且只能是具体数值，不能为符号．
若是积分上限或者下限是一个变化的值 （比如今天我做的一个很简单的计算就是这样的情况），则可以使用以下的方法：
>>result=maple（‘evalf’,'（Ei（1,y））’）
result=Ei（1,y）
>>y=2
y=2
>>result=subs（result）
result=Ei（1,2）
>>result=vpa（result）
result=.48900510708061119567239835228050e-1
>>result=str2num（maple（‘evalf’,'（Ei（1,2））’））
result=0.0489
比如其中的y就可以是一个变化的值，写成一个循环，从而计算相应的积分值．
Matlab也直接提供两种主要的数值积分函数：
quad 和 quad8．quad是变步长辛普生法，quad8是牛顿－柯特斯法，以我今天的例子持续来看，相差很小．
对被积分函数的定义同上，另外，还有一种办法，可以不用另外写一个m文件再调用，省下些小麻烦．可使用inline语句：
fxy=inline（‘exp（-a*x*y/b）’,‘x’,'y’）;
Matlab将字符串中的表达式录入内存，准备之后使用．这个语句有一个很大的缺点是，表达式中，除了变量之外，其他数值（如上式中的a和b）不能使用字母等符号表示，而必须为数值，即使你已经在之前定义过，也不行，因为 inline语句将引号中的表达式录为符号格式，其中任何的字母或者字母组合，都会被认为是变量，即使你在语句之后只写了真正的变量，程序还是会全部误读．这就在做积分时候带来错误，明明是常量的a和b，Matlab还是要求你给他们一个积分空间进行积分，从而出错．
定义完函数之后，直接使用quad 或者 quad8，进行数值积分：
F=quad（fx,1,10,1e-10）;%1,1０是积分上下限，1e-10是积分精度．
元宵节说了一堆电脑计算的事情，真的很奇怪．．．
祝大家元宵节快乐！
四数值积分
trapz（）用复化梯形公式求解定积分
命令格式:
I=trapz（x,y）
X是积分区间内的离散数据点向量，y是x的各分量函数值构成的向量，输出项I为积分的近似值
例：
计算积分exp（x）在0到1上
clear;
clc;
x=0:
0.2:
1;
y=exp（x）;
I=trapz（x,y）
quad（）采用自适应步长的辛普森公式求定积分
命令格式:
I=quad（fun,a,b,tol）
Fun是被积函数，a，b是积分区间的左右端点，tol为积分的精度要求，缺省值是1e-6，输出项为积分的近似值。
例：
计算积分exp（x）在0到1上
fun=inline（‘exp（x）’）;I=quad（fun,0,1）;
quadl（）采用自适应步长的Lobatto公式求定积分
命令格式：
I=quadl（fun,a,b,tol）
dblquad（）在矩形区域上求二重积分
命令格式:
I=dblquad（fun,a,b,c,d,tol,method）
Fun是二元被积函数f（x,y），a,b是变量x的上下限，cd是变量y的上下限，tol微积分的精度要求，缺省值是1e-6，method是积分方法，一种是@quad,另一种是@quadl，缺省值为@quad.
例子：
（1）计算二重积分x^2+y^2,x从0到1，y从0到1
I=dblquad（inline（'x.^2+y.^2'）,0,1,0,1）
syms x y;
I2=int（int（x^2+y^2,'y',0,1）,0,1）
（2）计算二重积分I=∫∫√（1-x^2-y^2）dxdy,其中D={（x,y）|x^2+y^2<=1}
clear;
clc;
fun1=inline（'sqrt（max（1-（x.^2+y.^2）,0））','x','y'）;
fun2=inline（'sqrt（1-（x.^2+y.^2））.*（x.^2+y.^2<=1）','x','y'）;
I1=dblquad（fun1,-1,1,-1,1）
I2=dblquad（fun2,-1,1,-1,1）
triplequad（）在立方体区域上求三重积分
命令格式:
I=triplequad（fun,a,b,c,d,e,f,tol,method）
Fun是三元被积分函数;a，b是变量x的上下限；cd是y的上下限；ef是变量z的上下限；tol为积分进度要求，缺省值为1e-6；method是积分方法，一种是@quad，另一种为@quadl
例子：
计算三重积分:
I=∫∫∫[ysin（x）+zcos（x）]dv,式中区域{（x,y,z）|0<=x<=pi,0<=y<=1,-1<=z<=1}
syms x yz;I4=int（int（int（y*sin（x）+z*cos（x））,'z',-1,1）,'y',0,1）,'x',0,pi）
I3=triplequad（inline（'y.*sin（x）+z.*cos（x）'）,0,pi,0,1,-1,1）
广义积分的数值方法：
（1）无界函数的广义积分。
对被积函数在积分区间某点附近无界的（奇点），然后用数值方法计算。
例子：
计算积分∫1/√x dx,x从0到1
I=quadl（inline（‘1./sqrt（x）’）,eps,1）
syms x
I2=int（1/sqrt（x）,0,1）
1梯形数值积分的MATLAB主程序
functionT=rctrap（fun,a,b,m）
%fun函数，a积分上限b积分下限m递归次数
n=1;h=b-a;T=zeros（1,m+1）;x=a;
T
（1）=h*（feval_r（fun,a）+feval_r（fun,b））/2;
fori=1:
m
h=h/2;n=2*n;s=0;
fork=1:
n/2
x=a+h*（2*k-1）;s=s+feval_r（fun,x）;
end
T（i+1）=T（i）/2+h*s;
end
T=T（1:
m）;
e.g
运行后屏幕显示精确值Fs，用rctrap计算的递归值T和T与精确值Fs的绝对误差wT
>>fun=@（x）exp（（-x^.2./2）./（sqrt（2*pi）））
T=rctrap（fun,0,pi/2,14）,symst
fi=int（exp（（-t^2）/2）/（sqrt（2*pi））,t,0,pi/2）;
Fs=double（fi）,wT=double（abs（fi-T））
fun=
@（x）exp（（-x^.2./2）./（sqrt（2*pi）））
T=
Columns1through7
1.4168 1.3578 1.3313 1.3195 1.3142 1.3119 1.3109
Columns8through14
1.3105 1.3103 1.3102 1.3102 1.3101 1.3101 1.3101
Fs=
0.4419
wT=
Columns1through7
0.9749 0.9159 0.8894 0.8776 0.8723 0.8700 0.8690
Columns8through14
0.8686 0.8684 0.8683 0.8683 0.8683 0.8682 0.8682
>>
2复合辛普森（Simpson）数值积分的MATLAB主程序
functiony=comsimpson（fun,a,b,n）
%fun函数a积分上限b积分下限n分割小区间数
z1=feval_r（fun,a）+feval_r（fun,b）;m=n/2;
h=（b-a）/（2*m）;x=a;
z2=0;z3=0;x2=0;x3=0;
fork=2:
2:
2*m
x2=x+k*h;z2=z2+2*feval_r（fun,x2）;
end
fork=3:
2:
2*m
x3=x+k*h;z3=z3+4*feval_r（fun,x3）;
end
y=（z1+z2+z3）*h/3;
由于Matlab自带了quad就是这个算法所以比较少自己编
3龙贝格数值积分的MATLAB主程序
function[RT,R,wugu,h]=romberg（fun,a,b,wucha,m）
%fun被积函数a，b积分上下限wucha两次相邻迭代绝对差值m龙贝格积分表最大行数
%RT龙贝格积分表R数值积分结果wucha误差估计h最小步长
n=1;h=b-a;wugu=1;x=a;k=0;RT=zeros（4,4）;
RT（1,1）=h*（feval_r（fun,a）+feval_r（fun,b））/2;
while（（wugu>wucha）&（k k=k+1; h=h/2;s=0;
forj=1:
n
x=a+h*（2*j-1）;s=s+feval_r（fun,x）;
end
RT（k+1,1）=RT（k,1）/2+h*s;n=2*n;
fori=1:
k
RT（k+1,i+1）=（（4^i）*RT（k+1,i）-RT（k,i））/（4^i-1）;
end
wugu=abs（RT（k+1,k）-RT（k+1,k+1））;
end
R=RT（k+1,k+1）;
e.g.
>>F=inline（'1./（1+x）'）;[RT,R,wugu,h]=romberg（F,0,1.5,1.e-8,13）
symsx
fi=int（1/（1+x）,x,0,1.5）;Fs=double（fi）,
wR=double（abs（fi-R））,wR1=wR-wugu
RT=
1.0500 0 0 0 0 0
0.9536 0.9214 0 0 0 0
0.9260 0.9168 0.9165 0 0 0
0.9187 0.9163 0.9163 0.9163 0 0
0.9169 0.9163 0.9163 0.9163 0.9163 0
0.9164 0.9163 0.9163 0.9163 0.9163 0.9163
R=
0.9163
wugu=
2.9436e-011
h=
0.0469
Fs=
环境影响的经济损益分析，也称环境影响的经济评价，即估算某一项目、规划或政策所引起的环境影响的经济价值，并将环境影响的经济价值纳入项目、规划或政策的经济费用效益分析中去，以判断这些环境影响对该项目：
规划或政策的可行性会产生多大的影响。
对负面的环境影响估算出的是环境费用，对正面的环境影响估算出的是环境效益。
0.9163
四、环境影响的经济损益分析
wR=
9.8007e-011
1.依法评价原则；
wR1=
6.8571e-011
（3）是否符合区域、流域规划和城市总体规划。
>>
第五章 环境影响评价与安全预评价4 复合梯形法
function[I,step]=CombineTraprl（f,a,b,eps）
%f 被积函数
%a,b积分上下限
%eps精度
%I 积分结果
%step积分的子区间数
if（nargin==3）
eps=1.0e-4;
end
n=1;
h=（b-a）/2;
I1=0;
I2=（subs（sym（f）,findsym（sym（f））,a）+subs（sym（f）,findsym（sym（f））,b））/h;
whileabs（I2-I1）>eps
n=n+1;
h=（b-a）/n;
I1=I2;
I2=0;
fori=0:
n-1
x=a+h*i;
x1=x+h;
I2=I2+（h/2）*（subs（sym（f）,findsym（sym（f））,x）+subs（sym（f）,findsym（sym（f））,x1））;
end
end
I=I2;
step=n;
5 辛普森法
function[I,step]=IntSimpson（f,a,b,type,eps）
%type=1辛普森公式
%type=2辛普森3/8公式
%type=3复合辛普森公式
四、安全预评价if（type==3&&nargin==4）
eps=1.0e-4; %缺省精度为0.0001
end
（1）资质等级。
评价机构的环评资质分为甲、乙两个等级。
环评证书在全国范围内使用，有效期为4年。
I=0;
switchtype
case1,
I=（（b-a）/6）*（subs（sym（f）,findsym（sym（f））,a）+...
4*subs（sym（f）,findsym（sym（f））,（a+b）/2）+...
subs（sym（f）,findsym（sym（f））,b））;
step=1;
case2,
I=（（b-a）/8）*（subs（sym（f）,findsym（sym（f））,a）+...
3*subs（sym（f）,findsym（sym（f））,（2*a+b）/3）+...
3*subs（sym（f）,findsym（sym（f））,（a+2*b）/3）+subs（sym（f）,findsym（sym（f））,b））;
step=1;
case3,
n=2;
h=（b-a）/2;
I1=0;
I2=（subs（sym（f）,findsym（sym（f））,a）+subs（sym（f）,findsym（sym（f））,b））/h;
while
升级会员 