2.3 数值积分2.3.1 一元函数的数值积分函数1 quad 、quadl 、quad8功能 数值定积分，自适应Simpleson 积分法。

格式 q = quad(fun,a,b) %近似地从a 到b 计算函数fun 的数值积分，误差为10-6。

若给fun 输入向量x ，应返回向量y ，即fun 是一单值函数。

q = quad(fun,a,b,tol) %用指定的绝对误差tol 代替缺省误差。

tol 越大，函数计算的次数越少，速度越快，但结果精度变小。

q = quad(fun,a,b,tol,trace,p1,p2,…) %将可选参数p1,p2,…等传递给函数fun(x,p1,p2,…)，再作数值积分。

若tol=[]或trace=[]，则用缺省值进行计算。

[q,n] = quad(fun,a,b,…) %同时返回函数计算的次数n… = quadl(fun,a,b,…) %用高精度进行计算，效率可能比quad 更好。

… = quad8(fun,a,b,…) %该命令是将废弃的命令，用quadl 代替。

例2-40>>fun = inline(‘3*x.^2./(x.^3-2*x.^2+3)’); equivalent to: function y=funn(x)y=3*x.^2./(x.^3-2*x.^2+3);>>Q1 = quad(fun,0,2) >>Q2 = quadl(fun,0,2)计算结果为：Q1 =3.7224 Q2 =3.7224补充：复化simpson 积分法程序程序名称 Simpson.m调用格式 I=Simpson('f_name',a,b,n)程序功能 用复化Simpson 公式求定积分值输入变量 f_name 为用户自己编写给定函数()y f x 的M 函数而命名的程序文件名 a 为积分下限b 为积分上限n 为积分区间[,]a b 划分成小区间的等份数 输出变量 I 为定积分值 程序function I=simpson(f_name,a,b,n) h=(b-a)/n; x=a+(0:n)*h; f=feval(f_name,x); N=length(f)-1;if N==1fprintf('Data has only one interval') return; end if N==2I=h/3*(f(1)+4*f(2)+f(3)); return; end if N==3I=3/8*h*(f(1)+3*f(2)+3*f(3)+f(4)); return; end I=0;if 2*floor(N/2)==NI=h/3*(2*f(N-2)+2*f(N-1)+4*f(N)+f(N+1)); m=N-3; else m=N; endI=I+(h/3)*(f(1)+4*sum(f(2:2:m))+2*f(m+1)); if m>2I=I+(h/3)*2*sum(f(3:2:m)); end例题 求0sin I xdx π=⎰。

解 先编制sin y x =的M 函数。

程序文件命名为sin_x.m 。

function y=sin_x(x) y=sin(x)将区间4等份，调用格式为I=Simpson (’sin _x’,0,pi,4)计算结果为y =0 0.7071 1.0000 0.7071 0.0000I =2.0046将区间20等份，调用格式为I=Simpson (’sin _x’,0,pi,20)计算结果为y =0 0.1564 0.3090 0.4540 0.5878 0.7071 0.80900.8910 0.9511 0.9877 1.0000 0.9877 0.9511 0.8910 0.8090 0.7071 0.5878 0.4540 0.3090 0.1564 0.0000I =2.0000重做上例2—40：simpson('funn',0,2,100)函数2 trapz功能 梯形法数值积分格式 T = trapz(Y) %用等距梯形法近似计算Y 的积分。

若Y 是一向量，则trapz(Y)为Y 的积分；若Y 是一矩阵，则trapz(Y)为Y 的每一列的积分；若Y 是一多维阵列，则trapz(Y)沿着Y 的第一个非单元集的方向进行计算。

T = trapz(X,Y) %用梯形法计算Y 在X 点上的积分。

若X 为一列向量，Y 为矩阵，且size(Y,1) = length(X)，则trapz(X,Y)通过Y 的第一个非单元集方向进行计算。

T = trapz(…,dim) %沿着dim 指定的方向对Y 进行积分。

若参量中包含X ，则应有length(X)=size(Y ,dim)。

例2-41>>X = -1:.1:1;>>Y = 1./(1+25*X.^2); >>T = trapz(X,Y)计算结果为：T =0.5492补充： 复化梯形积分法程序程序名称 Trapezd.m调用格式 I=Trapezd('f_name',a,b,n) 程序功能 用复化梯形公式求定积分值输入变量 f_name 为用户自己编写给定函数()y f x 的M 函数而命名的程序文件名 a 为积分下限b 为积分上限n 为积分区间[,]a b 划分成小区间的等份数 输出变量 I 为定积分值 程序function I=Trapezd(f_name,a,b,n) h=(b-a)/n;x=a+(0:n)*h; f=feval(f_name,x);I=h*(sum(f)-(f(1)+f(length(f)))/2); hc=(b-a)/100; xc=a+(0:100)*hc; fc=feval(f_name,xc); plot(xc,fc,'r');hold on ;title('Trapezoidal Rule');xlabel('x');ylabel('y'); plot(x,f);plot(x,zeros(size(x))) ; for i=1:n;plot([x(i),x(i)],[0,f(i)]); end补充例题 求0sin I xdx π=⎰。

解 先编制sin y x =的M 函数。

程序文件命名为sin_x.m 。

function y=sin_x(x) y=sin(x);将区间4等份，调用格式为I=Trapezd(’sin _x’,0,pi,4)计算结果为I=1.8961将区间20等份，调用格式为I=Trapezd(’sin _x’,0,pi,20)计算结果为I= 1.9959图A.5表示了复化梯形求积的过程。

（1）区间4等份（2）区间20等份重做上例2-41：function y=li2_41(x)y = 1./(1+25*x.^2);I=Trapezd(’li2_41’,-1,1,100)函数3 rat,rats功能有理分式近似。

虽然所有的浮点数值都是有理数，有时用简单的有理数字（分子与分母都是较小的整数）近似地表示它们是有必要的。

函数rat将试图做到这一点。

对于有连续出现的小数的数值，将会用有理式近似表示它们。

函数rats调用函数rat，且返回字符串。

格式[N,D] = rat(X) %对于缺省的误差1.e-6*norm(X(:),1)，返回阵列N与D，使N./D近似为X。

[N,D] = rat(X,tol) %在指定的误差tol范围内，返回阵列N与D，使N./D近似为X。

rat(X)、rat(X…) %在没有输出参量时，简单地显示x的连续分数。

S = rats(X,strlen) %返回一包含简单形式的、X中每一元素的有理近似字符串S，若对于分配的空间中不能显示某一元素，则用星号表示。

该元素与X中其他元素进行比较而言较小，但并非是可以忽略。

参量strlen为函数rats中返回的字符串元素的长度。

缺省值为strlen=13，这允许在78个空格中有6个元素。

S = rats(X) %返回与用MA TLAB命令format rat显示 X相同的结果给S。

例2-42>>s = 1-1/2+1/3-1/4+1/5-1/6+1/7>>format rat>>S1 = rats(s)>>S2 = rat(s)>>[n,d] = rat(s)>>PI1 = rats(pi)>>PI2 = rat(pi)计算结果为：s =0.7595S1 =319/420S2 =1 + 1/(-4 + 1/(-6 + 1/(-3 + 1/(-5))))n =319d =420PI1 =355/113PI2 =3 + 1/(7 + 1/(16))2.3.2 二元函数重积分的数值计算函数1 dblquad功能矩形区域上的二重积分的数值计算格式q = dblquad(fun,xmin,xmax,ymin,ymax) %调用函数quad在区域[xmin,xmax, ymin,ymax]上计算二元函数z=f(x,y)的二重积分。

输入向量x，标量y，则f(x,y)必须返回一用于积分的向量。

q = dblquad(fun,xmin,xmax,ymin,ymax,tol) %用指定的精度tol代替缺省精度10-6，再进行计算。

q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method) %用指定的算法method代替缺省算法quad。

method的取值有@quadl或用户指定的、与命令quad与quadl有相同调用次序的函数句柄。

q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method,p1,p2,…) %将可选参数p1,p2,..等传递给函数fun(x,y,p1,p2,…)。

若tol=[]，method=[]，则使用缺省精度和算法quad。

例2-43>>fun = inline(’y./sin(x)+x.*exp(y)’);>>Q = dblquad(fun,1,3,5,7)计算结果为：Q =3.8319e+0032.4 常微分方程数值解函数ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb功能常微分方程（ODE）组初值问题的数值解参数说明：solver为命令ode45、ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一。