数值积分与微分方程.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-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.7224Q2 =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.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.8910 0.8090 0.7071 0.5878 0.4540 0.3090 0.1564 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的第一个非单元集方向进行计算。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函数而命名的程序文件名a为积分下限b为积分上限n为积分区间划分成小区间的等份数输出变量 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补充例题 求。解 先编制的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) %返回与用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 = 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)必须返回一用于积分的向量。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-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)，或为包含一混合矩阵的方程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，用初始条件y0求解显式微分方程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（缺省值为每一元素为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）运用数学中的变量替换：yn=y(n-1),yn-1=y(n-2),y2=y1=y，把高阶（大于2阶）的方程（组）写成一阶微分方程组：，（3）根据（1）与（2）的结果，编写能计算导数的M-函数文件odefile。（4）将文件odefile与初始条件传递给求解器Solver中的一个，运行后就可得到ODE的、在指定时间区间上的解列向量y（其中包含y及不同阶的导数）。2求解器Solver与方程组的关系表见表2-3。表2-3函数指令含 义函 数含 义求解器Solverode23普通2-3阶法解ODEodefile包含ODE的文件ode23s低阶法解刚性ODE选项odeset创建、更改Solver选项ode23t解适度刚性ODEodeget读取Solver的设置值ode23tb低阶法解刚性ODE输出odeplotODE的时间序列图ode45普通4-5阶法解ODEodephas2ODE的二维相平面图ode15s变阶法解刚性ODEodephas3ODE的三维相平面图ode113普通变阶法解ODEodeprint在命令窗口输出结果3因为没有一种算法可以有效地解决所有的ODE问题，为此，MATLAB提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver。表2-4 不同求解器Solver的特点求解器SolverODE类型特点说明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算法；低精度当精度较低时，计算时间比ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比ode15s短4在计算过程中，用户可以对求解指令solver中的具体执行参数进行设置（如绝对误差、相对误差、步长等）。表2-5 Solver中options的属性属性名取值含义AbsTol有效值：正实数或向量缺省值：1e-6绝对误差对应于解向量中的所有元素；向量则分别对应于解向量中的每一分量RelTol有效值：正实数缺省值：1e-3相对误差对应于解向量中的所有元素。在每步(第k步)计算过程中，误差估计为：e(k)=max(RelTol*abs(y(k),AbsTol(k)NormControl有效值：on、off缺省值：off为on时，控制解向量范数的相对误差，使每步计算中，满足：norm(e)1缺省值：k = 1若k1，则增加每个积分步中的数据点记录，使解曲线更加的光滑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）用于解n个未知函数的方程组时，M函数文件中的待解方程组应以x的向量形式写成。例题A.7 解微分方程，其中首先，将导数表达式的右端编写成一个liA7.m函数程序：function yy=liA7(x,y)yy=sin(x);然后直接调用：x,y=ode23(liA7,0,pi,-1)plot(x,y)例4.用微分方程的数值解法求解微分方程，设自变量t的初始值为0,终值为3pi,初始条件y(0)=0,y(0)=0解：将高阶微分方程化为一阶微分方程组，即用变量代换：=这样，将导数表达式的右端编写成一个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，求二阶微分方程的数值解解：首先变量代换：这样，将导数表达式的右端编写成一个jie3.m函数程序function f=jie3(x,z)f=0 1;1/(2*x2)-1 -1/x*z;然后，在主程序中调用已有的数值积分函数进行积分：x,z=ode23(jie3,pi/2,pi,2;-2/pi)plot(x,z(:,1)例2-45 求解描述振荡器的经典的Ver der Pol微分方程y(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 = 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函数文件名x0,xend为起点和终点y0为已知初始值h为步长输出变量 X为离散的自变量Y为离散的函数值程序 functionx,y=Eulerpro(fxy,x0,y0, xend, h) n=fix(xend-x0)/h); y(1)=y0; x(1)=x0; for k=2:n x(k)=0; y(k)=0; end for 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; end plot(x,y)例题A.7 解微分方程，其中。解 先编制的M函数。程序文件命名为fxy.m。 function Z=fxy(x,y) Z=sin(x); 取步长0.1，调用格式为X,Y=Eulerpro(fxy,0,-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.8777 -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函数文件名x0,xend为起点