数值积分与微分方程教程.doc

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为用户自己编写给定函数的M函数而命名的程序文件名 a为积分下限 b为积分上限 n为积分区间划分成小区间的等份数 输出变量 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==1 fprintf(Data has only one interval) return; end if N==2 I=h/3*(f(1)+4*f(2)+f(3)); return; end if N==3 I=3/8*h*(f(1)+3*f(2)+3*f(3)+f(4)); return; end I=0; if 2*floor(N/2)==N I=h/3*(2*f(N-2)+2*f(N-1)+4*f(N)+f(N+1)); m=N-3; else m=N; end I=I+(h/3)*(f(1)+4*sum(f(2:2:m))+2*f(m+1)); if m2 I=I+(h/3)*2*sum(f(3:2:m)); end 例题 求。 解 先编制的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.0000 I = 2.0046 将区间20等份，调用格式为 I=Simpson(’sin_x’,0,pi,20) 计算结果为 y = 0 0.1564 0.3090 0.4540 0.5878 0.7071 0.8090 0.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.0000 I = 2.0000 重做上例2—40：simpson(funn,0,2,100) 函数2 trapz 功能 梯形法数值积分 格式 T = trapz(Y) %用等距梯形法近似计算Y的积分。若Y是一向量，则t