Matlab积分.doc

一.数值积分的实现方法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。function f=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);(2) 调用数值积分函数quad求定积分。S,n=quad(fesin,0,3*pi)S = 0.9008n = 772牛顿柯特斯法基于牛顿柯特斯法，MATLAB给出了quad8函数来求定积分。该函数的调用格式为：I,n=quad8(fname,a,b,tol,trace)其中参数的含义和quad函数相似，只是tol的缺省值取10-6。该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。(1) 被积函数文件fx.m。function f=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x);(2) 调用函数quad8求定积分。I=quad8(fx,0,pi)I = 2.4674分别用quad函数和quad8函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。调用函数quad求定积分：format long;fx=inline(exp(-x);I,n=quad(fx,1,2.5,1e-10)I = 0.28579444254766n = 65调用函数quad8求定积分：format long;fx=inline(exp(-x);I,n=quad8(fx,1,2.5,1e-10)I = 0.28579444254754n = 333被积函数由一个表格定义在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X,Y定义函数关系Y=f(X)。用trapz函数计算定积分。命令如下：x=1:0.01:2.5; Y=exp(-X); %生成函数关系数据向量trapz(X,Y)ans = 0.285796824163938.1.3 二重定积分的数值求解使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace)该函数求f(x,y)在a,bc,d区域上的二重定积分。参数tol，trace的用法与函数quad完全相同。 计算二重定积分(1) 建立一个函数文件fxy.m：function f=fxy(x,y)global ki;ki=ki+1; %ki用于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);(2) 调用dblquad函数求解。global ki;ki=0;I=dblquad(fxy,-2,2,-1,1)kiI = 1.57449318974494ki = 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求导数dpdpx=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数值积分的一些做法Filed under: 未分类 franz 9:31 pm 今天是元宵节，突然来讲Matlab的确很奇怪，连我自己也这样感觉 由于Matlab对于计算过程的细节要求定义详细，往往让人觉得使用不方便与同样很流行的Mathematic相比，Matlab更象是程序，而不是象Mathematic那么直观的写作业 Matlab同样提供符号积分（不定积分）和数值积分（定积分）两大功能符号积分使用int命令，配合之前的函数定义使用比如： a=; b=; function f = fun(x) f=a*x2+b; 将此文件寸为一个单独m文件，再在主程序中调用即可： F = int(fun,c,d); c和d为积分上下限，通常为符号，可使用syms c;语句定义。在完成符号积分之后可以对其附值，则就完成了数值积分的任务 有一点需要注意的是，通过这样过程求的数值，其格式为符号格式sym，不能做图，不能和数值进行运算处理方法是： vpa (F);求得位符号近似解，再用double ( vpa (F) );将其转换成Matlab默认的双精度位数就可以注意，直接做 double (F) 行不通，Matlab会提示你Conversion from sym to double is not possible，不知道not possbile和impossible有什么区别，呵呵 符号积分中一个最大的问题在于存在不可积的函数，比如十分常用的Boltzman分布，或者叫Arrhenius公式：exp(- q*V/ (k *T ) 插句话，我一直以为Arrhenius是德国人，知道前几天在和老师讨论中讲起，我老师很自豪的对我说，有个著名的瑞典化学家，得过Noble Prize，叫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=2y = 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和quad8quad是变步长辛普生法，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)计算二重积分x2+y2,x从0到1，y从0到1I=dblquad(inline(x.2+y.2),0,1,0,1)symsxy;I2=int(int(x2+y2,y,0,1),0,1)(2)计算二重积分I=(1-x2-y2)dxdy,其中D=(x,y)|x2+y2=1clear;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 fun=(x) exp(-x.2./2)./(sqrt(2*pi)T=rctrap(fun,0,pi/2,14), syms tfi=int(exp(-t2)/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 =Columns 1 through 7 1.4168 1.3578 1.3313 1.3195 1.3142 1.3119 1.3109Columns 8 through 14 1.3105 1.3103 1.3102 1.3102 1.3101 1.3101 1.3101Fs = 0.4419wT =Columns 1 through 7 0.9749 0.9159 0.8894 0.8776 0.8723 0.8700 0.8690Columns 8 through 14 0.8686 0.8684 0.8683 0.8683 0.8683 0.8682 0.86822 复合辛普森（Simpson）数值积分的MATLAB主程序function y=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;for k=2:2:2*mx2=x+k*h; z2= z2+2*feval_r(fun,x2);endfor k=3:2:2*mx3=x+k*h; z3= z3+4*feval_r(fun,x3);endy=(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(wuguwucha)&(km)|(k F=inline(1./(1+x); RT,R,wugu,h=romberg(F,0,1.5,1.e-8,13)syms xfi=int(1/(1+x),x,0,1.5); Fs=double(fi),wR=double(abs(fi-R), wR1= wR - wuguRT = 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.9163R = 0.9163wugu =2.9436e-011h = 0.0469Fs = 0.9163wR =9.8007e-011wR1 =6.8571e-0114复合梯形法function I,step = CombineTraprl(f,a,b,eps)%f 被积函数%a,b 积分上下限%eps 精度%I 积分结果%step 积分的子区间数if(nargin =3) eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h;while abs(I2-I1)eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=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); endendI=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.0001endI=0;switch type case 1, 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; case 2, 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; case 3, 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 abs(I2-I1)eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/6)*(subs(sym(f),findsym(sym(f),x)+. 4*subs(sym(f),findsym(sym(f),(x+x1)/2)+. subs(sym(f),findsym(sym(f),x1); end end I=I2; step=n;end6 牛顿科茨法function I = NewtonCotes(f,a,b,type)%type = 1 科茨公式%type = 2 牛顿科茨六点公式%type = 3 牛顿科茨七点公式I=0;switch type case 1, I=(b-a)/90)*(7*subs(sym(f),findsym(sym(f),a)+. 32*subs(sym(f),findsym(sym(f),(3*a+b)/4)+. 12*subs(sym(f),findsym(sym(f),(a+b)/2)+. 32*subs(sym(f),findsym(sym(f),(a+3*b)/4)+7*subs(sym(f),findsym(sym(f),b); case 2, I=(b-a)/288)*(19*subs(sym(f),findsym(sym(f),a)+. 75*subs(sym(f),findsym(sym(f),(4*a+b)/5)+. 50*subs(sym(f),findsym(sym(f),(3*a+2*b)/5)+. 50*subs(sym(f),findsym(sym(f),(2*a+3*b)/5)+. 75*subs(sym(f),findsym(sym(f),(a+4*b)/5)+19*subs(sym(f),findsym(sym(f),b); case 3, I=(b-a)/840)*(41*subs(sym(f),findsym(sym(f),a)+. 216*subs(sym(f),findsy