数值积分微分方程

《数值积分微分方程》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

-2.3数值积分2.3.1一元函数的数值积分函数1quad、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)’);equivalentto:functiony=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为定积分值程序functionI=simpson(f_name,a,b,n)h=(b-a)/n;x=a+(0:n)*h;f=feval(f_name,x);N=length(f)-1;.-- -ifN==1fprintf('Datahasonlyoneinterval')return;endifN==2I=h/3*(f(1)+4*f(2)+f(3));return;endifN==3I=3/8*h*(f(1)+3*f(2)+3*f(3)+f(4));return;endI=0;if2*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));ifm>2I=I+(h/3)*2*sum(f(3:2:m));end例题求。解先编制的M函数。程序文件命名为sin_x.m。functiony=sin_x(x)y=sin(x)将区间4等份，调用格式为I=Simpson(’sin_x’,0,pi,4)计算结果为y=00.70711.00000.70710.0000I=.-- -2.0046将区间20等份，调用格式为I=Simpson(’sin_x’,0,pi,20)计算结果为y=00.15640.30900.45400.58780.70710.80900.89100.95110.98771.00000.98770.95110.89100.80900.70710.58780.45400.30900.15640.0000I=2.0000重做上例2—40：simpson('funn',0,2,100)函数2trapz功能梯形法数值积分格式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为用户自己编写给定函数的M函数而命名的程序文件名a为积分下限b为积分上限n为积分区间划分成小区间的等份数输出变量I为定积分值程序functionI=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');holdon;title('TrapezoidalRule');xlabel('x');ylabel('y');plot(x,f);plot(x,zeros(size(x)));fori=1:n;plot([x(i),x(i)],[0,f(i)]);end补充例题求。解先编制的M函数。程序文件命名为sin_x.m。functiony=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：functiony=li2_41(x)y=1./(1+25*x.^2);I=Trapezd(’li2_41’,-1,1,100)函数3rat,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命令formatrat显示 X相同的结果给S。例2-42>>s=1-1/2+1/3-1/4+1/5-1/6+1/7>>formatrat>>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二元函数重积分的数值计算函数1dblquad.-- -功能矩形区域上的二重积分的数值计算格式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之一。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-3～10-6计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gear’s反向数值微分；精度中等若ode45失效时，可尝试使用ode23s刚性一步法；2阶Rosebrock算法；低精度当精度较低时，计算时间比ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比ode15s短4．在计算过程中，用户可以对求解指令solver中的具体执行参数进行设置（如绝对误差、相对误差、步长等）。表2-5Solver中options的属性属性名取值含义AbsTol有效值：正实数或向量缺省值：1e-6绝对误差对应于解向量中的所有元素；向量则分别对应于解向量中的每一分量RelTol有效值：正实数缺省值：1e-3相对误差对应于解向量中的所有元素。在每步(第k步)计算过程中，误差估计为：e(k)<=max(RelTol*abs(y(k)),AbsTol(k))NormControl有效值：on、off缺省值：off为‘on’时，控制解向量范数的相对误差，使每步计算中，满足：norm(e)<=max(RelTol*norm(y),AbsTol)Events有效值：on、off为‘on’时，返回相应的事件记录OutputFcn有效值：odeplot、odephas2、odephas3、odeprint缺省值：odeplot若无输出参量，则solver将执行下面操作之一：画出解向量中各元素随时间的变化；画出解向量中前两个分量构成的相平面图；.-- -画出解向量中前三个分量构成的三维相空间图；随计算过程，显示解向量OutputSel有效值：正整数向量缺省值：[]若不使用缺省设置，则OutputFcn所表现的是那些正整数指定的解向量中的分量的曲线或数据。若为缺省值时，则缺省地按上面情形进行操作Refine有效值：正整数k>1缺省值：k=1若k>1，则增加每个积分步中的数据点记录，使解曲线更加的光滑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函数程序：functionyy=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函数程序functionxdot=exf(t,x)u=1-(t.^2)/(2*pi);xdot=[01;-10]*x+[01]'*u;然后，在主程序中调用已有的数值积分函数进行积分：clf;t0=0;tf=3*pi;x0=[0;0][t,x]=ode23('exf',[t0,tf],x0)y=x(:,1).-- -例5，求二阶微分方程的数值解解：首先变量代换：这样，将导数表达式的右端编写成一个jie3.m函数程序functionf=jie3(x,z)f=[01;1/(2*x^2)-1-1/x]*z;然后，在主程序中调用已有的数值积分函数进行积分：[x,z]=ode23('jie3',[pi/2,pi],[2;-2/pi])plot(x,z(:,1))例2-45求解描述振荡器的经典的VerderPol微分方程y(0)=1，y’(0)=0令x1=y，x2=dy/dt，则dx1/dt=x2dx2/dt=μ(1-(x1)^2)*x2-x1编写函数文件verderpol.m：functionxprime=verderpol(t,x)globalMUxprime=[x(2);MU*(1-x(1)^2)*x(2)-x(1)];再在命令窗口中执行：>>globalMU>>MU=7;>>Y0=[1;0]>>[t,x]=ode45(‘verderpol’,0,40,Y0);>>x1=x(:,1);x2=x(:,2);>>plot(t,x1,t,x2)图形结果为图2-20。图2-20VerderPol微分方程图A.4.1改进的Euler法程序程序名称Eulerpro.m调用格式[X,Y]=Eulerpro('fxy',x0,y0,xend,h)程序功能解常微分方程.-- -输入变量fxy为用户编写给定函数的M函数文件名x0,xend为起点和终点y0为已知初始值h为步长输出变量X为离散的自变量Y为离散的函数值程序function[x,y]=Eulerpro(fxy,x0,y0,xend,h)n=fix((xend-x0)/h);y(1)=y0;x(1)=x0;fork=2:nx(k)=0;y(k)=0;endfori=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解微分方程，其中。解先编制的M函数。程序文件命名为fxy.m。functionZ=fxy(x,y)Z=sin(x);取步长0.1，调用格式为[X,Y]=Eulerpro(‘fxy’,0,-1,pi,0.1)计算结果如图A.6所示。.-- -图A.6微分方程求解结果x=00.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.20001.30001.40001.50001.60001.70001.80001.90002.00002.10002.20002.30002.40002.50002.60002.70002.80002.90003.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.02830.12790.22620.32220.41500.50360.58720.66490.73590.79960.85530.90250.94060.96930.9883A.4.2Runge-Kutta法程序程序名称RungKt4.m调用格式[X,Y]＝RungKt4('fxy',x0,y0,xend,M)程序功能解常微分方程输入变量fxy为用户编写给定函数的M函数文件名x0,xend为起点和终点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;fori=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)+(k1+2*k2+2*k3+k4)/6;.-- -endplot(X,Y)例题A.8解微分方程，其中。解先编制的M函数。文件名取为fxy.m。functionZ=fxy(x,y)Z=sin(x);取步长数为30，调用格式为[X,Y]=Rungkt4('fxy',0,-1,pi,30)计算结果如图A.7所示。X=00.10470.20940.31420.41890.52360.62830.73300.83780.94251.04721.15191.25661.36141.46611.57081.67551.78021.88501.98972.09442.19912.30382.40862.51332.61802.72272.82742.93223.03693.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.10450.00000.10450.20790.30900.40670.50000.58780.66910.74310.80900.86600.91350.95110.97810.99451.0000.--