整理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); 求得32位符号近似解,再用double(vpa(F));将其转换成Matlab默认的双精度位数就可以.注意,直接做double(F)行不通,Matlab会提示你"Conversionfromsymtodoubleisnotpossible",不知道"notpossbile"和"impossible"有什么区别,呵呵.
符号积分中一个最大的问题在于存在不可积的函数,比如十分常用的Boltzman分布,或者叫Arrhenius公式:
exp(-q*V/(k*T)).
插句话,我一直以为Arrhenius是德国人,知道前几天在和老师讨论中讲起,我老师很自豪的对我说,有个著名的瑞典化学家,得过NoblePrize,叫Arrhenius你知道不知道,才发现这个小问题,3滴汗,搞的我老师也很有挫败感...
对于不可积的函数,使用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,10是积分上下限,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