Login
升级VIP 登录 注册 安全退出
当前位置: 首页 > word文档 > 述职汇报 > [整理]Matlab积分.,matlab整理数据

[整理]Matlab积分.,matlab整理数据

收藏

本作品内容为[整理]Matlab积分.,格式为 doc ,大小 46592 KB ,页数为 20页

[整理]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.5x).sin(x+pi/6);(2)调用数值积分函数quad求定积分。[S,n]=quad(\'fesin\',0,3pi)S=0.9008n=772.牛顿-柯特斯法基于牛顿-柯特斯法,MATLAB给出了quad8函数来求定积分。该函数的调用格式为:[I,n]=quad8(\'fname\',a,b,tol,trace)其中参数的含义和quad函数相似,只是tol的缺省值取10-6。\x7f该函数可以更精确地求出定积分的值,且一般情况下函数调用的步数明显小于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.28579444254766n=65调用函数quad8求定积分:formatlong;fx=inline(\'exp(-x)\');[I,n]=quad8(fx,1,2.5,1e-10)I=0.28579444254754n=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.285796824163938.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)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+2x.^2-x+12)+(x+5).^(1/6)+5x+2\');g=inline(\'(3x.^2+4x-1)./sqrt(x.^3+2x.^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数值积分的一些做法Filedunder:未分类—franz@9:31pm今天是元宵节,突然来讲Matlab的确很奇怪,连我自己也这样感觉.由于Matlab对于计算过程的细节要求定义详细,往往让人觉得使用不方便.与同样很流行的Mathematic相比,Matlab更象是程序,而不是象Mathematic那么直观的写作业.Matlab同样提供符号积分(不定积分)和数值积分(定积分)两大功能.符号积分使用int命令,配合之前的函数定义使用.比如:a=…;b=…;functionf=fun(x)f=ax^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(-qV/(kT)).插句话,我一直以为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=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和quad8.quad是变步长辛普生法,quad8是牛顿-柯特斯法,以我今天的例子持续来看,相差很小.对被积分函数的定义同上,另外,还有一种办法,可以不用另外写一个m文件再调用,省下些小麻烦.可使用inline语句:fxy=inline(‘exp(-axy/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到1I=dblquad(inline(\'x.^2+y.^2\'),0,1,0,1)symsxy;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}symsxyz;I4=int(int(int(ysin(x)+zcos(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/√xdx,x从0到1I=quadl(inline(‘1./sqrt(x)’),eps,1)symsxI2=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:mh=h/2;n=2n;s=0;fork=1:n/2x=a+h(2k-1);s=s+feval_r(fun,x);endT(i+1)=T(i)/2+hs;endT=T(1:m);e.g运行后屏幕显示精确值Fs,用rctrap计算的递归值T和T与精确值Fs的绝对误差wT>>fun=@(x)exp((-x^.2./2)./(sqrt(2pi)))T=rctrap(fun,0,pi/2,14),symstfi=int(exp((-t^2)/2)/(sqrt(2pi)),t,0,pi/2);Fs=double(fi),wT=double(abs(fi-T))fun=@(x)exp((-x^.2./2)./(sqrt(2pi)))T=Columns1through71.41681.35781.33131.31951.31421.31191.3109Columns8through14--------------------------1.31051.31031.31021.31021.31011.31011.3101Fs=0.4419wT=Columns1through70.97490.91590.88940.87760.87230.87000.8690Columns8through140.86860.86840.86830.86830.86830.86820.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)/(2m);x=a;z2=0;z3=0;x2=0;x3=0;fork=2:2:2mx2=x+kh;z2=z2+2feval_r(fun,x2);endfork=3:2:2mx3=x+kh;z3=z3+4feval_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((wugu>wucha)&(k>F=inline(\'1./(1+x)\');[RT,R,wugu,h]=romberg(F,0,1.5,1.e-8,13)symsxfi=int(1/(1+x),x,0,1.5);Fs=double(fi),wR=double(abs(fi-R)),wR1=wR-wuguRT=1.0500000000.95360.921400000.92600.91680.91650000.91870.91630.91630.916300--------------------------0.91690.91630.91630.91630.916300.91640.91630.91630.91630.91630.9163R=0.9163wugu=2.9436e-011h=0.0469Fs=0.9163wR=9.8007e-011wR1=6.8571e-011>>--------------------------4复合梯形法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;whileabs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;fori=0:n-1x=a+hi;x1=x+h;I2=I2+(h/2)(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1));endend--------------------------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.0001endI=0;switchtypecase1,I=((b-a)/6)(subs(sym(f),findsym(sym(f)),a)+...4subs(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)+...3subs(sym(f),findsym(sym(f)),(2a+b)/3)+...3subs(sym(f),findsym(sym(f)),(a+2b)/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;whileabs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;fori=0:n-1x=a+hi;x1=x+h;I2=I2+(h/6)(subs(sym(f),findsym(sym(f)),x)+...4subs(sym(f),findsym(sym(f)),(x+x1)/2)+...subs(sym(f),findsym(sym(f)),x1));endendI=I2;step=n;end--------------------------6牛顿-科茨法functionI=NewtonCotes(f,a,b,type)%type=1科茨公式%type=2牛顿-科茨六点公式%type=3牛顿-科茨七点公式I=0;switchtypecase1,I=((b-a)/90)(7subs(sym(f),findsym(sym(f)),a)+...32subs(sym(f),findsym(sym(f)),(3a+b)/4)+...12subs(sym(f),findsym(sym(f)),(a+b)/2)+...32subs(sym(f),findsym(sym(f)),(a+3b)/4)+7subs(sym(f),findsym(sym(f)),b));case2,I=((b-a)/288)(19subs(sym(f),findsym(sym(f)),a)+...75subs(sym(f),findsym(sym(f)),(4a+b)/5)+...50subs(sym(f),findsym(sym(f)),(3a+2b)/5)+...50subs(sym(f),findsym(sym(f)),(2a+3b)/5)+...75subs(sym(f),findsym(sym(f)),(a+4b)/5)+19subs(sym(f),findsym(sym(f)),b));case3,I=((b-a)/840)(41subs(sym(f),findsym(sym(f)),a)+...216subs(sym(f),findsym(sym(f)),(5a+b)/6)+...27subs(sym(f),findsym(sym(f)),(2a+b)/3)+...272subs(sym(f),findsym(sym(f)),(a+b)/2)+...--------------------------27subs(sym(f),findsym(sym(f)),(a+2b)/3)+...216subs(sym(f),findsym(sym(f)),(a+5b)/6)+41subs(sym(f),findsym(sym(f)),b));end7高斯公式functionI=IntGauss(f,a,b,n,AK,XK)if(n<5&&nargin==4)AK=0;XK=0;elseXK1=((b-a)/2)XK+((a+b)/2);I=((b-a)/2)sum(AK.subs(sym(f),findsym(f),XK1));endta=(b-a)/2;tb=(a+b)/2;switchncase0,I=2tasubs(sym(f),findsym(sym(f)),tb);case1,I=ta(subs(sym(f),findsym(sym(f)),ta0.5773503+tb)+...subs(sym(f),findsym(sym(f)),-ta0.5773503+tb));--------------------------case2,I=ta(0.55555556subs(sym(f),findsym(sym(f)),ta0.7745967+tb)+...0.55555556subs(sym(f),findsym(sym(f)),-ta0.7745967+tb)+...0.88888889subs(sym(f),findsym(sym(f)),tb));case3,I=ta(0.3478548subs(sym(f),findsym(sym(f)),ta0.8611363+tb)+...0.3478548subs(sym(f),findsym(sym(f)),-ta0.8611363+tb)+...0.6521452subs(sym(f),findsym(sym(f)),ta0.3398810+tb)...+0.6521452subs(sym(f),findsym(sym(f)),-ta0.3398810+tb));case4,I=ta(0.2369269subs(sym(f),findsym(sym(f)),ta0.9061793+tb)+...0.2369269subs(sym(f),findsym(sym(f)),-ta0.9061793+tb)+...0.4786287subs(sym(f),findsym(sym(f)),ta0.5384693+tb)...+0.4786287subs(sym(f),findsym(sym(f)),-ta0.5384693+tb)+...0.5688889subs(sym(f),findsym(sym(f)),tb));end8区间逐次分半梯形法functionq=DblTraprl(f,a,A,b,B,m,n)if(m==1&&n==1)%梯形公式--------------------------q=((B-b)(A-a)/4)(subs(sym(f),findsym(sym(f)),{a,b})+...subs(sym(f),findsym(sym(f)),{a,B})+...subs(sym(f),findsym(sym(f)),{A,b})+...subs(sym(f),findsym(sym(f)),{A,B}));else%复合梯形公式C=4ones(n+1,m+1);C(1,:)=2;C(:,1)=2;C(n+1,:)=2;C(:,m+1)=2;C(1,1)=1;C(1,m+1)=1;C(n+1,1)=1;C(n+1,m+1)=1;%C矩阵endF=zeros(n+1,m+1);q=0;fori=0:nforj=0:mx=a+i(A-a)/n;y=b+j(B-b)/m;F(i+1,j+1)=subs(sym(f),findsym(sym(f)),{x,y});q=q+F(i+1,j+1)C(i+1,j+1);--------------------------endendq=((B-b)(A-a)/4/m/n)q;9区间逐次分半布尔法function[I,step]=DDBuer(f,a,b,eps)if(nargin==3)eps=1.0e-4;end;n=1;h=b-a;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h/2;tol=1;whiletol>epsn=n+1;h=h/4;I1=I2;I2=0;fori=0:(4^(n-1)-1)x=a+h4i;x1=x+h;x2=x1+h;x3=x2+h;x4=x3+h;--------------------------I2=I2+(2h/45)(7subs(sym(f),findsym(sym(f)),x)+...32subs(sym(f),findsym(sym(f)),x1)+...12subs(sym(f),findsym(sym(f)),x2)+...32subs(sym(f),findsym(sym(f)),x3)+...7subs(sym(f),findsym(sym(f)),x4));endtol=abs(I2-I1);endI=I2;step=n;10自适应求积法functionI=SmartSimpson(f,a,b,eps)if(nargin==3)eps=1.0e-4;end;e=5eps;I=SubSmartSimpson(f,a,b,e);functionq=SubSmartSimpson(f,a,b,eps)QA=IntSimpson(f,a,b,1,eps);QLeft=IntSimpson(f,a,(a+b)/2,1,eps);QRight=IntSimpson(f,(a+b)/2,b,1,eps);if(abs(QLeft+QRight-QA)<=eps)q=QA;else--------------------------q=SubSmartSimpson(f,a,(a+b)/2,eps)+SubSmartSimpson(f,(a+b)/2,b,eps);ende.g>>SmartSimpson(\'xsin(x)\',0,1)ans=0.3011>>-------------',)


  • 编号:1700785315
  • 分类:述职汇报
  • 软件: wps,office word
  • 大小:20页
  • 格式:docx
  • 风格:商务
  • PPT页数:46592 KB
  • 标签:

广告位推荐

相关述职汇报更多>