1、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,p
2、2,)，再作数值积分。若tol=或 trace=，则用缺省值进行计算。q,n = quad(fun,a,b,) %同时返回函数计算的次数 n = quadl(fun,a,b,) %用高精度进行计算，效率可能比 quad 更好。 = quad8(fun,a,b,) %该命令是将废弃的命令，用 quadl 代替。例 2-40fun = 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
3、.7224Q2 =3.7224补充：复化 simpson 积分法程序程序名称 Simpson.m调用格式 I=Simpson(f_name,a,b,n)程序功能 用复化 Simpson 公式求定积分值输入变量 f_name 为用户自己编写给定函数 的 M 函数而命名的程序文件名()yfxa 为积分下限b 为积分上限n 为积分区间 划分成小区间的等份数,ab输出变量 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 ha
4、s only one interval)return;endif N=2I=h/3*(f(1)+4*f(2)+f(3);return;endif N=3I=3/8*h*(f(1)+3*f(2)+3*f(3)+f(4);return;endI=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;elsem=N;endI=I+(h/3)*(f(1)+4*sum(f(2:2:m)+2*f(m+1);if m2I=I+(h/3)*2*sum(f(3:2:m);end例题 求 。0sinIxd解 先编制 的 M 函数。程序文件命
5、名为 sin_x.m。yfunction 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.8090 0.8910 0.9511 0.9877 1.0000 0.9877 0.9511 0.89100.8090 0.7071 0.5878 0.4540 0.3090 0.15
6、64 0.0000I =2.0000重做上例 240：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 的第一个非单元集
7、方向进行计算。T = trapz(,dim) %沿着 dim 指定的方向对 Y 进行积分。若参量中包含 X，则应有 length(X)=size(Y,dim)。例 2-41X = -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 为用户自己编写给定函数 的 M 函数而命名的程序文件名()yfxa 为积分下限b 为积分上限n 为积分区间 划分成小区间的等份数,ab输出变量
8、 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补
9、充例题 求 。0sinIxd解 先编制 的 M 函数。程序文件命名为 sin_x.m。yfunction 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,ra
10、ts功能 有理分式近似。虽然所有的浮点数值都是有理数，有时用简单的有理数字（分子与分母都是较小的整数）近似地表示它们是有必要的。函数 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 = ra
11、ts(X,strlen) %返回一包含简单形式的、X 中每一元素的有理近似字符串 S，若对于分配的空间中不能显示某一元素，则用星号表示。该元素与 X 中其他元素进行比较而言较小，但并非是可以忽略。参量 strlen 为函数 rats 中返回的字符串元素的长度。缺省值为 strlen=13，这允许在 78 个空格中有 6 个元素。S = rats(X) %返回与用 MATLAB 命令 format rat 显示 X 相同的结果给S。例 2-42s = 1-1/2+1/3-1/4+1/5-1/6+1/7format ratS1 = rats(s)S2 = rat(s)n,d = rat(s)PI1
12、 = rats(pi)PI2 = rat(pi)计算结果为：s =0.7595S1 =319/420 S2 =1 + 1/(-4 + 1/(-6 + 1/(-3 + 1/(-5)n =319 d =420 PI1 =355/113 PI2 =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
13、)必须返回一用于积分的向量。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,
14、)。若 tol=，method=，则使用缺省精度和算法 quad。例 2-43fun = 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 之一。Odefun 为显式常微分方程 y=f(t,y)，或为包含一混合矩
15、阵的方程 M(t,y)*y=f(t,y)。命令 ode23 只能求解常数混合矩阵的问题；命令 ode23t 与 ode15s 可以求解奇异矩阵的问题。Tspan 积分区间（即求解区间）的向量 tspan=t0,tf。要获得问题在其他指定时间点t0,t1,t2,上的解，则令 tspan=t0,t1,t2,tf（要求是单调的） 。Y0 包含初始条件的向量。Options 用命令 odeset 设置的可选积分参数。P1,p2, 传递给函数 odefun 的可选参数。格式 T,Y = solver(odefun,tspan,y0) %在区间 tspan=t0,tf上，从 t0 到 tf，用初始条件 y
16、0 求解显式微分方程 y=f(t,y)。对于标量 t 与列向量 y，函数f=odefun(t,y)必须返回一 f(t,y)的列向量 f。解矩阵 Y 中的每一行对应于返回的时间列向量 T 中的一个时间点。要获得问题在其他指定时间点t0,t1,t2,上的解，则令 tspan=t0,t1,t2,tf（要求是单调的） 。T,Y = solver(odefun,tspan,y0,options) %用参数 options（用命令 odeset 生成）设置的属性（代替了缺省的积分参数） ，再进行操作。常用的属性包括相对误差值 RelTol（缺省值为 1e-3）与绝对误差向量 AbsTol（缺省值为每一元素
17、为 1e-6） 。T,Y =solver(odefun,tspan,y0,options,p1,p2) 将参数 p1,p2,p3,等传递给函数odefun，再进行计算。若没有参数设置，则令 options=。1求解具体 ODE 的基本过程：（1）根据问题所属学科中的规律、定律、公式，用微分方程与初始条件进行描述。F(y,y,y,y(n),t) = 0y(0)=y0,y(0)=y1,y(n-1)(0)=yn-1而 y=y;y(1);y(2);,y(m-1)，n 与 m 可以不等（2）运用数学中的变量替换：y n=y(n-1),yn-1=y(n-2),y2=y1=y，把高阶（大于 2 阶）的方程（
18、组）写成一阶微分方程组： ，)y,t(f,tyn12 n10n210y)(（3）根据（1）与（2）的结果，编写能计算导数的 M-函数文件 odefile。（4）将文件 odefile 与初始条件传递给求解器 Solver 中的一个，运行后就可得到 ODE的、在指定时间区间上的解列向量 y（其中包含 y 及不同阶的导数） 。2求解器 Solver 与方程组的关系表见表 2-3。表 2-3函数指令 含 义 函 数 含 义ode23 普通 2-3 阶法解 ODE odefile 包含 ODE 的文件ode23s 低阶法解刚性 ODE odeset 创建、更改 Solver 选项ode23t 解适度刚
19、性 ODE 选项 odeget 读取 Solver 的设置值ode23tb 低阶法解刚性 ODE odeplot ODE 的时间序列图ode45 普通 4-5 阶法解 ODE odephas2 ODE 的二维相平面图ode15s 变阶法解刚性 ODE odephas3 ODE 的三维相平面图求解器Solverode113 普通变阶法解 ODE输出odeprint 在命令窗口输出结果3因为没有一种算法可以有效地解决所有的 ODE 问题，为此， MATLAB 提供了多种求解器 Solver，对于不同的 ODE 问题，采用不同的 Solver。表 2-4 不同求解器 Solver 的特点求解器 So
20、lver ODE 类型 特点 说明ode45 非刚性 一步算法；4，5 阶 Runge-Kutta方程；累计截断误差达(x) 3 大部分场合的首选算法ode23 非刚性 一步算法；2，3 阶 Runge-Kutta方程；累计截断误差达(x) 3 使用于精度较低的情形ode113 非刚性 多步法；Adams 算法；高低精度均可到 10-310 -6 计算时间比 ode45 短ode23t 适度刚性 采用梯形算法 适度刚性情形ode15s 刚性 多步法；Gears 反向数值微分；精度中等 若 ode45 失效时，可尝试使用ode23s 刚性 一步法；2 阶 Rosebrock 算法；低精度 当精度
21、较低时，计算时间比ode15s 短ode23tb 刚性 梯形算法；低精度 当精度较低时，计算时间比ode15s 短4在计算过程中，用户可以对求解指令 solver 中的具体执行参数进行设置（如绝对误差、相对误差、步长等） 。表 2-5 Solver 中 options 的属性属性名 取值 含义AbsTol 有效值：正实数或向量缺省值：1e-6 绝对误差对应于解向量中的所有元素；向量则分别对应于解向量中的每一分量RelTol 有效值：正实数缺省值：1e-3相对误差对应于解向量中的所有元素。在每步(第 k 步) 计算过程中，误差估计为：e(k)1缺省值：k = 1 若 k1，则增加每个积分步中的数
22、据点记录，使解曲线更加的光滑Jacobian 有效值：on、off缺省值：off 若为on时，返回相应的 ode 函数的 Jacobi 矩阵Jpattern 有效值：on、off缺省值：off 为on时，返回相应的 ode 函数的稀疏 Jacobi 矩阵Mass有效值：none 、M、M(t)、M(t,y)缺省值：noneM：不随时间变化的常数矩阵M(t)：随时间变化的矩阵M(t,y)：随时间、地点变化的矩阵MaxStep 有效值：正实数缺省值：tspans/10 最大积分步长注意：（1）求微分方程数值解的函数命令中，函数 odefun 必须以 dx/dt 为输出量，以 t,x 为输入量。（2
23、）用于解 n 个未知函数的方程组时，M 函数文件中的待解方程组应以 x 的向量形式写成。例题 A.7 解微分方程 ，其中sinyx0,1y首先，将导数表达式的右端编写成一个 liA7.m 函数程序：function yy=liA7(x,y)yy=sin(x);然后直接调用：x,y=ode23(liA7,0,pi,-1)plot(x,y)例 4.用微分方程的数值解法求解微分方程 ，设自变量 t 的初始值为 0,21ty终值为 3pi,初始条件 y(0)=0,y(0)=0解：将高阶微分方程化为一阶微分方程组，即用变量代换： yx21= 2121txx )21(0102tx这样，将导数表达式的右端编
24、写成一个 exf.m 函数程序function xdot=exf(t,x)u=1-(t.2)/(2*pi);xdot=0 1;-1 0*x+0 1*u;然后，在主程序中调用已有的数值积分函数进行积分：clf;t0=0;tf=3*pi;x0=0;0t,x=ode23(exf,t0,tf,x0)y=x(:,1)例 5，求二阶微分方程 的数值解2)(,)2(,0)1(22 yyxyx解：首先变量代换： z21 1221)(zxz这样，将导数表达式的右端编写成一个 jie3.m 函数程序function f=jie3(x,z)f=0 1;1/(2*x2)-1 -1/x*z;然后，在主程序中调用已有的数
25、值积分函数进行积分：x,z=ode23(jie3,pi/2,pi,2;-2/pi)plot(x,z(:,1)例 2-45 求解描述振荡器的经典的 Ver der Pol 微分方程 0)1(22ydtdtyy(0)=1，y(0)=0令 x1=y，x2=dy/dt，则dx1/dt = x2dx2/dt = (1-(x1)2)*x2-x1编写函数文件 verderpol.m：function xprime = verderpol(t,x)global MUxprime = x(2);MU*(1-x(1)2)*x(2)-x(1);再在命令窗口中执行：global MUMU = 7;Y0=1;0t,x
26、= ode45(verderpol,0,40,Y0);x1=x(:,1);x2=x(:,2);plot(t,x1,t,x2)图形结果为图 2-20。图 2-20 Ver der Pol 微分方程图A.4.1 改进的 Euler 法程序程序名称 Eulerpro.m调用格式 X,Y=Eulerpro(fxy,x0,y0, xend ,h) 程序功能 解常微分方程输入变量 fxy 为用户编写给定函数 的 M 函数文件名(,)yfxx0,xend 为起点和终点y0 为已知初始值h 为步长输出变量 X 为离散的自变量Y 为离散的函数值程序functionx,y=Eulerpro(fxy,x0,y0,
27、xend, h)n=fix(xend-x0)/h);y(1)=y0;x(1)=x0;for k=2:nx(k)=0;y(k)=0;endfor i=1:(n-1)x(i+1)=x0+i*h;y1=y(i)+h*feval(fxy,x(i),y(i);y2=y(i)+h*feval(fxy,x(i+1),y1);y(i+1)=(y1+y2)/2;endplot(x,y)例题 A.7 解微分方程 ，其中 。sinyx0,1y解 先编制 的 M 函数。程序文件命名为 fxy.m。ifunction Z=fxy(x,y)Z=sin(x);取步长 0.1，调用格式为X,Y=Eulerpro(fxy,0,
28、-1, pi ,0.1)计算结果如图 A.6 所示。图 A.6 微分方程求解结果x =0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.70000.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.50001.6000 1.7000 1.8000 1.9000 2.0000 2.1000 2.2000 2.30002.4000 2.5000 2.6000 2.7000 2.8000 2.9000 3.0000y =-1.0000 -0.9950 -0.9801 -0.9554 -0.9211 -0.877
29、7 -0.8255 -0.7650-0.6970 -0.6219 -0.5407 -0.4541 -0.3629 -0.2681 -0.1707 -0.07150.0283 0.1279 0.2262 0.3222 0.4150 0.5036 0.5872 0.66490.7359 0.7996 0.8553 0.9025 0.9406 0.9693 0.9883A.4.2 Runge-Kutta 法程序程序名称 RungKt4.m调用格式 X,YRungKt4(fxy,x0,y0,xend,M)程序功能 解常微分方程输入变量 fxy 为用户编写给定函数 的 M 函数文件名(,)yfxx0,x
30、end 为起点和终点y0 为已知初始值M 为步长数输出变量 X 为离散的自变量Y 为离散的函数值程序function X,Y=Rungkt4(fxy,x0,y0,xend,M)h=(xend-x0)/M;X=zeros(1,M+1);Y=zeros(1,M+1);X=x0:h:xend;Y(1)=y0;for i=1:Mk1=h*feval(fxy,X(i),Y(i);k2=h*feval(fxy,X(i)+h/2,Y(i)+k1/2);k3=h*feval(fxy,X(i)+h/2,Y(i)+k2/2);k4=h*feval(fxy,X(i)+h,Y(i)+k3);Y(i+1)=Y(i)+(
31、k1+2*k2+2*k3+k4)/6;endplot(X,Y)例题 A.8 解微分方程 ，其中 。sinyx0,1y解 先编制 的 M 函数。文件名取为 fxy.m。ifunction Z=fxy(x,y)Z=sin(x);取步长数为 30，调用格式为X,Y= Rungkt4(fxy,0,-1,pi,30)计算结果如图 A.7 所示。X =0 0.1047 0.2094 0.3142 0.4189 0.5236 0.6283 0.7330 0.8378 0.9425 1.0472 1.1519 1.2566 1.3614 1.4661 1.5708 1.6755 1.7802 1.8850 1
32、.9897 2.0944 2.1991 2.3038 2.4086 2.5133 2.6180 2.7227 2.8274 2.9322 3.0369 3.1416图 A.7 微分方程的求解过程Y =-1.0000 -0.9945 -0.9781 -0.9511 -0.9135 -0.8660 -0.8090 -0.7431 -0.6691 -0.5878 -0.5000 -0.4067 -0.3090 -0.2079 -0.1045 0.0000 0.1045 0.2079 0.3090 0.4067 0.5000 0.5878 0.6691 0.7431 0.8090 0.8660 0.9135 0.9511 0.9781 0.9945 1.0000
