高斯-勒让德数值积分Matlab代码.doc

高斯-勒让德数值积分Matlab代码codefunction ql,Ak,xk=guasslegendre(fun,a,b,n,tol)% 高斯-勒让德数值积分% 参数说明% fun：积分表达式，可以是函数句柄、inline函数、匿名函数、字符串表达式，但是必须可以接受矢量输入% a,b：积分上下限，注意积分区间太大会降低精度，此时建议使用复化求积公式，默认-1 1% n：积分阶数，可以任意正整数，但是不建议设置过大，大不一定能得到更好的精度，默认7% tol：积分精度，默认1e-6% ql：积分结果% Ak：求积系数% xk：求积节点，满足ql=sum(Ak.*fun(xk)% 举例说明% fun=(x)exp(x).*cos(x); % 必须可以接受矢量输入% quadl(fun,0,pi) % 调用MATLAB内部积分函数检验% ql,Ak,xk=guasslegendre(fun,0,pi)% 高斯积分说明% 1.高斯积分是精度最高的插值型数值积分，具有2n+1阶精度，并且高斯积分总是稳定。一般插值型数值积分精度为至少n阶，且具有高阶不稳定性% 2.高斯求积节点就是对应n阶正交多项式的根，构建首项系数为1的正交多项式参见/thread-4798-1-1.html% 3.高斯求积系数，可以由Lagrange多项式插值系数进行积分得到% 4.由高斯求积节点为根构成的多项式，与任何不超过n阶的多项式正交% 勒让德正交多项式说明% 1.区间-1,1上关于权rho(x)=1的正交多项式系，满足% |- 2/(2j+1)(i=j)% (Pi(x)*Pj(x),-1,1)= | % |- 0 (ij)% 2.勒让德正交多项式的通式为：P0=1, Pn=1/(2n*n!) * diff(x2-1)n,n)(n=1,2,.)% 3.关于高斯-勒让德的数值积分的求积系数和节点，由于表达式不便于书写，感兴趣的网友可以参考相关书籍% by dynamic of Matlab技术论坛% see also % contact me % 2010-01-15 23:05:33% 输入参数有效性检验if nargin=1 a=-1;b=1;n=7;tol=1e-8;elseif nargin=3 n=7;tol=1e-8;elseif nargin=4 tol=1e-8;elseif nargin=2|nargin5 error(The Number of Input Arguments Is Wrong!);end% 计算求积节点syms xp=sym2poly(diff(x2-1)(n+1),n+1)/(2n*factorial(n);tk=roots(p); % 求积节点% 计算求积系数Ak=zeros(n+1,1);for i=1:n+1 xkt=tk; xkt(i)=; pn=poly(xkt); fp=(x)polyval(pn,x)/polyval(pn,tk(i); Ak(i)=quadl(fp,-1,1,tol); % 求积系数end% 积分变量代换，将a,b变换到-1,1xk=(b-a)/2*tk+(b+a)/2;% 检验积分函数fun有效性fun=fcnchk(fun,vectorize);% 计算变量代换之后积分函数的值fx=fun(xk)*(b-a)/2;% 计算积分值ql=sum(Ak.*fx);/codeMatlab只能直接计算内外积分限都是已知实数的二重积分比如 INT_(a)(b)INT_cdf(x,y)dxdy，其中a,b,c,d都是已知实数，则可以用Matlab指令:dblquad(f(x,y),c,d,a,d);但如果内积分限是外积分的积分变量函数时，Matlab不能直接运算。比如INT_(a)(b)INT_g1(y)g2(y)f(x,y)dxdy, 其中外积分限a,b都是已知实数, 而内积分限g1(y), g2(y)都是外积分变量y的函数解法： 把内积分变量x转换成y和一个新定义变量v的函数： x=(1-v)*g1(y)+v*g2(y),这样，原来的积分公式可以变为：INT_(a)(b)INT_01)f(1-v)*g1(y)+v*g2(y), y)*g2(y)-g1(y)dvdy这时内积分变量是v，内积分上下限是常数0和1. 积分函数变成f(1-v)*g1(y)+v*g2(y),y)*g2(y)-g1(y) ，注意这里f(x,y)中的x已经转换成(1-v)*g1(y)+v*g2(y)。这样，这个二重积分可以用Matalb指令算出:dblquad(f(1-v)*g1(y)+v*g2(y),y)*(g2(y)-g1(y) ,0,1,a,b);用matlab求定积分的三个实例代码一、符号积分符号积分由函数int来实现。该函数的一般调用格式为：int(s)：没有指定积分变量和积分阶数时，系统按findsym函数指示的默认变量对被积函数或符号表达式s求不定积分；int(s,v)：以v为自变量，对被积函数或符号表达式s求不定积分；int(s,v,a,b)：求定积分运算。a,b分别表示定积分的下限和上限。该函数求被积函数在区间a,b上的定积分。a和b可以是两个具体的数，也可以是一个符号表达式，还可以是无穷(inf)。当函数f关于变量x在闭区间a,b上可积时，函数返回一个定积分结果。当a,b中有一个是inf时，函数返回一个广义积分。当a,b中有一个符号表达式时，函数返回一个符号函数。例：求函数x2+y2+z2的三重积分。内积分上下限都是函数，对z积分下限是sqrt(x*y)，积分上限是x2*y；对y积分下限是sqrt(x)，积分上限是x2；对x的积分下限1，上限是2，求解如下：syms x y z %定义符号变量F2=int(int(int(x2+y2+z2,z,sqrt(x*y),x2*y),y,sqrt(x),x2),x,1,2) %注意定积分的书写格式F2 =1610027357/6563700-6072064/348075*2(1/2)+14912/4641*2(1/4)+64/225*2(3/4) %给出有理数解VF2=vpa(F2) %给出默认精度的数值解VF2 =224.92153573331143159790710032805二、数值积分1.数值积分基本原理求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。2.数值积分的实现方法基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为：I,n=quad(fname,a,b,tol,trace)基于变步长、牛顿柯特斯(Newton-Cotes)法，MATLAB给出了quadl函数来求定积分。该函数的调用格式为：I,n=quadl(fname,a,b,tol,trace)其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。例：求函数exp(-x*x)的定积分，积分下限为0，积分上限为1。fun=inline(exp(-x.*x),x); %用内联函数定义被积函数fnameIsim=quad(fun,0,1) %辛普生法Isim =0.746824180726425IL=quadl(fun,0,1) %牛顿柯特斯法IL =0.746824133988447三、梯形法求向量积分trapz(x,y)梯形法沿列方向求函数Y关于自变量X的积分(向量形式，数值方法)。d=0.001;x=0:d:1;S=d*trapz(exp(-x.2)S= 0.7468或：format long gx=0:0.001:1; %x向量，也可以是不等间距y=exp(-x.2); %y向量，也可以不是由已知函数生成的向量S=trapz(x,y); %求向量积分S =0.746824071499185int的积分可以是定积分，也可以是不定积分（即有没有积分上下限都可以积）可以得到解析的解，比如你对x2积分，得到的结果是1/3*x3，这是通过解析的方法来解的。如果int(x2,x,1,2)得到的结果是7/3quad是数值积分，它只能是定积分（就是有积分上下限的积分），它是通过simpson数值积分来求得的（并不是通过解析的方法得到解析解，再将上下限代入，而是用小梯形的面积求和得到的）。如果f=inline(x.2);quad(f,1,2)得到的结果是2.333333，这个数并不是7/3int是符号解，无任何误差，唯一问题是计算速度；quad是数值解，有计算精度限制，优点是总是能有一定的速度，即总能在一定时间内给出一个一定精度的解。 对于y=exp(-(x.2+x+1)/(1+x)，被积函数之原函数无封闭解析表达式,符号计算无法解题，这是符号计算有限性，结果如下： syms xy=exp(-(x.2+x+1)/(1+x)s=int(y,x,0,inf) y =exp(-x2-x-1)/(1+x)Warning: Explicit integral could not be found. In at 58 s =int(exp(-x2-x-1)/(1+x),x = 0 . Inf)只有通过数值计算解法 dx=0.05; %采样间隔 x=0:dx:1000; %数值计算适合于有限区间上,取有限个采样点，只要终值足够大，精度不受影响y=exp(-(x.2+x+1)./(1+x);S=dx*cumtrapz(y); %计算区间内曲线下图形面积,为小矩形面积累加得S(end) ans =0.5641 %所求定积分值或进行编程，积分上限人工输入，程序如下：%表达式保存为函数文件function y=fxy(x)y=exp(-(x.2+x+1)./(1+x); % save fxy.m% main -主程序clear,clch=.001;p=0;a=0;R=input(请输入积分上限，R=)while a0 e1=2; end end if n2=0 e2=1; else if n20 e2=2; end end if n3=0 e3=1; else if n30 e3=2; end end F=(x,y,z)sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos(n3*pi*z/9); S=triplequad(F,-6,6,-5.5,5.5,-4.5,4.5) %求三重数值积分将以上代码保存为Fn.m程序文件，即m文件，然后运行： Fn(1,1,1)S =866.9655 M-file:functionf=intg(a,b,y0,z0,h0)symsxyzhu=x+y+z+h;v=int(u,x,a,b);f=subs(subs(subs(v,y,y0),z,z0),h,h0)intg(0,3