matlab自学经典课件4

第四章数值计算功能,4.1多项式计算4.2数值插值和曲线拟合polyfitintep14.3函数最大值4.4非线性方程问题求解fzero4.5常微分方程的数值求解4.6数值微分与积分4.7概率统计,4.1多项式运算,1多项式表示法2多项式运算函数,多项式求根多项式求值多项式乘法和除法多项式的微积分,1多项式表示法,(1)直接法:,多项式,多项式系数用行向量表示，按降幂排列；缺某幂次项，其幂次项系数为零。,P=a0,a1,.an-1,an,符串形式,x多项式变量,y=poly2str(P,x),(2)字符串表示法：,y=poly2sym(P,x),(3)符号多项式表示法：,完整形式,p=1,3,0,4,5y=poly2str(p,x)y=x4+3x3+4x+5symsxy=poly2sym(p,x)y=x4+3*x3+4*x+5,p=1,3,0,4,5p=13045y=poly2str(p,x)y=x4+3x3+4x+5rp=roots(p)rp=-3.23460.5594+1.1980i0.5594-1.1980i-0.8843,P=1,8,3,4,-10;pp=poly2str(P,X)pp=X4+8X3+3X2+4X-10,例,2.多项式的运算函数,多项式的值；根和微分；拟合曲线；部分分式,polyval,polyval,polyvalm,covolution,（1）多项式求根（roots）,（1）多项式的根=0n次多项式n个根，或实根，或若干对共轭复根。,r=roots(P),P：多项式系数向量；r：n个根b(1),b(2),b(n)的向量。,r=roots(P)r=-7.6998-0.5572+1.1335i-0.5572-1.1335i0.8141,例,求多项式x4+8x3+3x2+4x-10的根,P=1,8,3,4,-10;r=roots(P),多项式根的r向量：r(1),r(2),r(n),（2）根行向量创建多项式（poly）,P=poly(r),PP=poly(r)PP=1.00008.00003.00004.0000-10.0000formatrat避免浮点显示p=poly(r)p=1834-10,多项式的向量系数：,(3)多项式求值（polyval）,Y=polyval(P,x),若x数值，则求多项式在该点x的值；若X是向量或矩阵，每元素x(i,j)多项式求值。,代数多项式P求值,数组乘方运算:,Y=a0*X.n+a1*X.(n-1)an+1,p=1,3,0,2;poly2str(p,x)ans=x3+3x2+2,a=1,2;3,4a=1234,y=polyval(p,a)y=62256114,y=polyval(p,2)y=22,例,例,Y=polyvalm(P,X),矩阵多项式求值（polyvalm）,Y=a0*Xn+a1*Xn-1an+1,X必为方阵，作自变量代入多项式求值；结果阶数还是保持方阵阶数。,相当于矩阵X乘方运算：,设A为方阵，P代表多项式x3-5x2-2：pp=polyval(P,A)的含义pp=A.*A.*A-5*A.*A-2Pm=polyvalm(P,A)的含义：Pm=A*A*A-5*A*A-2,p=1,-5,0,-2;a=2,4;1,0a=2410pp=polyval(p,a)pp=-14-18-6-2pm=polyvalm(p,a)pm=-18-8-2-14,矩阵乘方运算:,数组乘方运算:,例,例,14,(4)多项式加减运算,相同次数多项式，加减其系数向量,不同次数,低次多项式中高次项用0补足。,p2=0 x3-0 x2+2x+1,p1+p2=2x3-x2+2x+4,2,-1,0,3,2,1,0,0,2,1,2,-1,2,4,p1=2x3-x2+3,例,例,(5)多项式乘法-向量卷积,w=conv(P1,P2),多项式x4+8x3-10与2x2-x+3的乘积。,p1=1,8,0,0,-10;p2=2,-1,3;p12=conv(p1,p2)p12=Columns1through5215-524-20Columns6through710-30,例,例,h=3,2,1,-2,1,0,-4,0,3;x=1,-2,3,-4,3,2,1;y=conv(h,x);n=0:14;stem(n,y);%杆图xlabel(Timeindexn);%标坐标轴ylabel(Amplitude);title(OutputObtainedbyConvolution)%图形标题grid;,卷积是两个变量在某范围内相乘后求和的结果。卷积的变量是序列x(n)和h(n)，y=conv(x,h)来实现卷级的，输出的结果个数等于x的长度与h的长度之和减去1。卷积公式：z(n)=x(n)*y(n)=x(m)y(n-m)dm.,例,(6)多项式除法-向量解卷积,Q,r=deconv(P1,P2),Q：多项式P1除以P2的商式，r：P1除以P2的余式，Q和r仍是多项式系数向量。,P1=conv(P2,Q)+r,deconv是conv的逆函数,求多项式x4+8x3-10除以2x2-x+3的结果。,q,r=deconv(p12,p2)q=1800-10r=Columns1through500000Columns6through700q,r=deconv(p2,p12)q=0r=2-13,例,(7)多项式的一阶导数（polyder）,k=polyder(P),k=polyder(P,Q),多项式乘积PQ的导函数：,多项式P的导函数：,导函数的分子存入p，分母存入q。,p,q=polyder(P,Q),多项式除法P/Q的导函数：,例：求有理分式的导数。命令如下：P=1;Q=1,0,5;p,q=polyder(P,Q),例,例,21,(8)多项式积分(polyint),I=polyint(2,-1,0,3,5);,例：p(x)=2x3-x2+3,求常数项5,k=polyint(P,c),k=polyint(P),不定积分，常数项取c,不定积分，常数项取零,例,4.2曲线拟合和数据插值,4.2.1曲线拟合polyfit4.2.2数据插值interp,拟合成多项式,推算未知点值,4.2.1.多项式曲线拟合（polyfit）,最小二乘法(误差平方和最小)；x、y已知数据的横、纵坐标；n为多项式次数;,p=polyfit(x,y,n),1.N次多项式曲线拟合,S结构体数组(struct)，估计预测误差，含R，df和normr。R：先输入x构建范德蒙矩阵V，后QR分解，得上三角矩阵。df：自由度，df=length(y)-(n+1)。df0时，超定方程组求解，拟合点数比未知数(p(1)p(n+1)多。normr：标准偏差、残差范数，normr=norm(y-V*p)，此处的p为求解之后的数值。,p,S=polyfit(x,y,n),对不同年份人口数据分别进行1次、2次和9次多项式拟合（polyfit），用poly2str表示多项式完整形式法；1次、2次和9次多项式估计2000年人口(polyval)，结合美国2000年人口普查截至2000年4月1日美国人口2.81421906亿数据绘制19002000年间的时间-人口数(polyval)曲线，要求用plot不同线型（LineSpec）绘制原始数据点、拟合的1次、2次和9次多项式曲线，标注图例，说明是否阶数越到高越好,美国从1900年1990年每隔10年进行人口普查的数据如下表所示：（百万）,例,plot(x,y,LineSpec)线型属性,例如r-.*、-.r*、*-.r等形式是等效的,属性控制符顺序不受限制，可有可无。,plot(x,y,-gh),clearn=1900:10:1990;r=7599,9197,10571,12320,13166,15069,17932,20321,22650,24963;nrf1=polyfit(n,r,1)%1次多项式拟合nrfs1=poly2str(nrf1,n)%1次多项式完整字符串表达式nrf2=polyfit(n,r,2)%2次多项式拟合nrf9=polyfit(n,r,9)%9次多项式拟合nrf10=polyfit(n,r,10)%10次多项式拟合nrf1_2000=polyval(nrf1,2000)%一次拟合得到的2000年人口数nrf2_2000=polyval(nrf2,2000)%2次拟合得到的2000年人口数nrf9_2000=polyval(nrf9,2000)%9次拟合得到的2000年人口数n20=1900:4:2000;%1900年到2000年的线性等分数组nrfv1=polyval(nrf1,n20);%从1900年到2000年间1次拟合人口数nrfv2=polyval(nrf2,n20);%从1900年到2000年间2次拟合人口数nrfv9=polyval(nrf9,n20);%从1900年到2000年间9次拟合人口数plot(n,r,or,n20,nrfv1,-p4=polyfit(x,y,11)y4=polyval(p4,xx);Warning:Polynomialisnotunique;degree=numberofdatapoints.,例,holdon;plot(x,y,*);plot(xx,y1,r-);plot(xx,y2,g-);plot(xx,y3,b-),多项式次数要适当，过低误差大，过高波动明显,已知：10个实验数据，分别采用2次、5次、9次和11次拟合，选出最佳拟合次数,3.函数拟合,clearx=0:5;y=0.2097,0.3523,0.4339,0.5236,0.7590,0.8998;ly=log(y);p=polyfit(x,ly,1);b=p(1);la=p(2);a=exp(la);Xx=linspace(0,5,30);Yy=a*exp(b*Xx);plot(x,y,ro,Xx,Yy),例,例,插值构造的函数必须通过已知数据点；拟合则不要求,只要均方差最小。,4.数据拟合和插值,相同点,都需根据已知数据构造函数；可使用得到函数计算未知点的函数值。,不同点,4.2.2数据插值,构造函数近似表达式的方法。常用多项式作插值函数,称多项式插值。多项式插值基本思想：高次多项式或分段的低次多项式为被插值函数f(x)的近似解析表达式。,1.插值法,根据已知输入/输出数据集和当前输入估计输出值,2.插值函数,多项式插值法：拉格朗日插值法牛顿插值法埃尔米特插值法分段插值法样条插值法等,3.一维插值,函数y=f(x)一维插值原理：调用格式：,yi=interp1(x,y,xi,method),已知X,Y两个等长向量，采样点和样本值；Xi是向量或标量，欲插值点；Yi是一个与Xi等长的插值结果。,分段线性插值linear：插值点样本值落在两相邻数据点的连线上.(缺省).最邻近插值法nearest：两点间插值点对应值就是离两点最近那点值。3次多项式插值cubic：已知数据求3次多项式，用多项式求插值。3次样条插值spline：每分段内构造一3次多项式，插值函数除满足插值条件外节点处光滑条件，再按照样条函数插值。,插值方法method,yy=spline(x,y,xx),clearx=0:1:10;y=sin(x);xx=0:0.2:10;yy=interp1(x,y,xx)subplot(1,4,1)plot(x,y,-*,xx,yy,dr);subplot(1,4,2);y2=interp1(x,y,xx,nearest);plot(x,y,-*,xx,y2,dr);subplot(1,4,3);y3=interp1(x,y,xx,cubic);plot(x,y,-*,xx,y3,dr)subplot(1,4,4);y4=interp1(x,y,xx,spline);plot(x,y,-*,xx,y4,dr),X1取值范围不能超出X给定范围，否则会给出“NaN”错误。,例,采用差值法估计美国的人口数量2000年人口及19001996年每隔4的人口数(interp1)；将上题原始数据点、2阶拟合曲线和插值曲线绘制在同一图，标注坐标轴为时间和人口（xlabel、ylabel）、图形标题（title）为插值和拟合和图例legend,例,clearn=1900:10:1990;r=7599,9197,10571,12320,13166,15069,17932,20321,22650,24963;n4=1900:4:1996;%19001996年每隔4的线性等分数组nrf2=polyfit(n,r,2)%2次多项式拟合nrfv2=polyval(nrf2,n4)nr4=interp1(n,r,n4,spline)%插值19001996年每隔4的人口数nr2000=interp1(n,r,2000)%线性插值2000年人口数nr2000=interp1(n,r,2000,spline)%插值2000年人口数plot(n,r,or,n4,nrfv2,-*k,n4,nr4,-dg)%绘图legend(原始数据,2阶多项式拟合曲线,插值曲线)%图例,例,38.02522138.07526738.12531838.17548738.22581538.275150238.325316238.375620738.425841438.475815638.525597538.575347338.625160138.67571038.72538638.77532538.82524738.87523238.92520738.97518239.025157,插值解法XRD1：nano-sic/Alcomposiste,例,plot(xrd1(:,1),xrd1(:,2),-*),loadxrd1.txt,xx=linspace(38,39,80);yy=interp1(xrd1(:,1),xrd1(:,2),xx,spline);plot(xx,yy,-*g),xrd1(:,1),)x=38.025038.075038.125038.175038.225038.275038.325038.375038.425038.475038.525038.575038.625038.675038.725038.775038.825038.875038.925038.975039.0250,xrd1(:,2),y=22126731848781515023162620784148156597534731601710386325247232207182157,X1取值范围不能超出X给定范围，否则会给出“NaN”错误。,xx=linspace(38,39,80);yy=interp1(xrd1(:,1),xrd1(:,2),xx);yy=interp1(xrd1(:,1),xrd1(:,2),xx,);,三次样条插值,4.二维插值,函数z=f(x,y)二维插值的原理：,向量的采样点X,Y，与采样点对应函数值Z；向量或标量的欲插值点Xi,Yi，插值结果Zi。指定插值方法method：linearnearestcubicsplineXi,Yi取值不超X,Y给定范围，否则“NaN”错误。,zi=interp2(x,y,z,xi,yi,method),函数interp2二维插值：,例,运行结果如下图所示。,向量的采样点X,Y，与采样点对应函数值Z；向量或标量的欲插值点Xi,Yi，插值结果Zi。指定插值方法method：linear(缺省)nearestcubicv4griddata算法Xi,Yi取值不超X,Y给定范围，否则“NaN”错误。,xi,yi,zi=griddata(x,y,z,xi,yi,method),函数griddata二维栅格插值：,输入参量（XI,YI）通常是规则的格点,clearp=rand(30,3)X=p(:,1)Y=p(:,2)Z=p(:,3)Xi,Yi=meshgrid(linspace(min(X),max(X),100)Xi,Yi,Zi=griddata(X,Y,Z,Xi,Yi,cubic)mesh(Xi,Yi,Zi),Xi,Yi,Zi=griddata(X,Y,Z,Xi,Yi,v4)mesh(Xi,Yi,Zi),随机生成30*3矩阵,例,1.函数最小值和非线性方程求解,调用格式一样,x,fval,exitflag,output=fminbnd(fun,x1,x2,options),最小值的解或零点根x、该点x的函数值fval、程序退出标志exitflag输出结果output、待求根函数fun、初始值x0、左右边界x1,x2,fval,exitflag,output和options可省略,x,fval,exitflag,output=fminsearch(fun,x0,options),x,fval,exitflag,output=fzero(fun,x0,options),迭代解法,一元函数最小值：,多元函数最小值：,x,fval,exitflag,output,jacobian=fsolve(fun,x0,options),非线性方程组求解,非线性方程求解,option：配置优化参数,optimset函数定义参数;,options=optimset(param1,value1.),optimset(display,off),4.3函数最小值,输出量x：最小值的x解或零点的根(自变量值)；fval=f(x)：最小值或零点时函数值f(x)；fun：待求最小值或根函数f(x)：函数名(少用)；字符串表达式，匿名对象、函数句柄、内联函数；exitflag：程序退出标志,1收敛到解x；output：选定输出结果；x0：搜索初始值；x1，x2：左右边界,output=迭代次数iterations:24函数评价次数funcCount:48算法algorithm:bisection,interpolation对分插值退出信息message:Zerofoundintheinterval-28,190.51,option：配置优化参数,optimset函数定义参数;最优化工具箱的(20多个)选项设定;optimset命令:将option显示出来;改变其中某个选项;display:选项函数调用时中间结果显示方式off不显示;iter每步都显示;final只显示最终结果。,optimset(display,off),options=optimset(param1,value1.),2.待求函数fun的形式,函数文件名(少用)；myfun函数句柄匿名对象字符串表达式内联函数,fun,+函数文件名,(变量表)函数表达式,expression,inline(expression),inline,本质和function一样，它直接内嵌在命令行里，不用单独定义function。,函数表达式,内联函数inline,inline(expression,x),f=inline(x-10.x+2,x),f(3)ans=-995,functionf=myfunn(x)f=x-10.x+2;end,函数文件名myfunn,函数文件名(少用)；函数句柄匿名对象字符串表达式内联函数,Fh=myfunn,Fn=(x)x-10x+2,Fs=x-10x+2,Fi=inline(x-10x+2),whosNameSizeBytesClassAttributesFh1x132function_handleFi1x11144inlineFn1x132function_handleFs1x816charfhh1x18double,x,fval,exitflag,output=fminbnd(fun,x1,x2,options),3.求一元函数最小值,求区间x1x2内函数f(x)最小值时x：,X：返回最小值的解;fval=F(x)，即最小值；fun：用于定义需求解函数;x1，x2：范围；option：最优化工具箱的选项设定,x=2.5148f=-0.4993e=1o=iterations:13funcCount:14algorithm:goldensectionsearch,parabolicinterpolationmessage:优化已终止:当前的x满足使用1.000000e-04的OPTIONS.Tol.,functionf=mmfun(x)f=exp(-0.1*x)*(sin(x)2)-0.5*(x+0.1)*sin(x)end,x,f,e,o=fminbnd(mmfun,-10,10,optimset(display,iter),求,在（-10，10）,最小值,例,最优化的结构体output，最优化取下列域：algorithm使用算法funcCount函数个数评估intervaliterations发现区间的迭代次数iterations发现零值点迭代次数message退出信息,x,f=fminbnd(mmfun,6,10,optimset(display,iter),x=8.0236f=-3.5680,x=-10:0.05:10;y=exp(-0.1*x).*(sin(x).2)-0.5*(x+0.1).*sin(x);plot(x,y)plot(x,y),Fminbnd在指定区域找到真解,求,在（6，10）,最小值,x,f,e,o=fminbnd(mmfun,6,10,optimset(display,iter),x,f=fminbnd(exp(-0.1*x)*(sin(x)2)-0.5*(x+0.1)*sin(x),6,10),x,f=fminbnd(inline(exp(-0.1*x)*(sin(x)2)-0.5*(x+0.1)*sin(x),6,10),最小,最小,x,fval=fminbnd(-exp(-0.1*x)*(sin(x)2)+0.5*(x+0.1)*sin(x),-8,-2)x=-4.830203748934343fval=-3.947274022275747,最大值求解：-f(x)在区间(a,b)上的最小值,x=-10:0.05:10;y=-exp(-0.1*x).*(sin(x).2)+0.5*(x+0.1).*sin(x);plot(x,y)plot(x,y),最大,例,3.求多元函数的最小值,在初始x0附近寻找局部最小值；使用options选项来指定优化器的参数；fval,exitflag,output和options可以省略。,x,fval,exitflag,output=fminsearch(fun,x0,options),著名RosenbrocksBanana测试函数,理论极小值x=1,y=1.求极小值点,x=11fval=0exitflag=1iterations:24funcCount:49algorithm:Nelder-Meadsimplexdirectsearchmessage:优化已终止:当前的x满足使用1.000000e-04的OPTIONS.TolX的终止条件，.,例,x,fval,exitflag,output=fminsearch(fxy,1;1),functionf=fxy(x)f=100*(x(2)-x(1).2).2+(1-x(1).2end,x,fval,exitflag,output=fminsearch(100*(x(2)-x(1).2).2+(1-x(1).2,1;1),方程,代数方程,微分方程,线性方程,非线性方程,常微分方程,偏微分方程,4.5常微分方程的求解,2.5线性方程组求解,4.4非线性方程求解,4.4非线性方程数值求解,4.4.1单变量非线性方程求解(fzero)4.4.2非线性方程组的求解(fsolve),迭代解法,4.4.1单变量非线性方程求解,(fzero),x,fval,exitflag,output=fzero(fun,x1,x2,options),x,fval,exitflag,output=fzero(fun,x0,options),2.调用格式：,函数是否穿越横轴决定零点数值解，求方程f(x)=0根,1.解方程原理：,fval,exitflag,output和options可省略,3.求根函数fun【f(x)】的调用,求f(x)=x-10 x+2=0在x0=0.5附近的根。,函数名：,fzero（funname，x0）,例,x,f,e,o=fzero(fun,0.5,optimset(display,iter),x=0.3758；f=0；e=1o=intervaliterations:8iterations:5funcCount:21algorithm:bisection,interpolation二分法message:在区间0.34,0.613137中发现零,建立函数文件fun.m。functionfv=fun(x)fv=x-10.x+2;,x=-4:0.05:1;f=x-10.x+2;plot(x,f)x,y=ginput(1)ans=-2.0012,x,fval=fzero(fun,-1)x=-1.989761447718557fval=0,(1)建立函数文件fun.m。functionfv=fun(x)fv=x-10.x+2;(2)调用fzero函数求根。x,fval,e,output=fzero(fun,0.5)x=0.3758fval=0,x,y=fzero(fun,-3,0.9)x=-1.9898y=0,例,多个根只求x0最近的根,x,f,e,o=fzero(fun,8,optimset(display,iter)x=NaNfval=NaNexitflag=-3,正在退出fzero:终止搜索包含符号变化的区间因为在搜索期间遇到NaN或Inf函数值。请检查函数或使用其他起始值重试。o=intervaliter:22iterations:0funcCount:44algorithm:bisection,interpolationmessage:正在退出fzero:终止搜索包含符号变化的区间因为在搜索期间遇到NaN或Inf函数值。,例,无根为Nan,4.多项式求的根（roots）,roots(1,0,-2,-5)ans=2.0946-1.0473+1.1359i-1.0473-1.1359ix,fval=fzero(x3-2*x-5,3)x=2.0946fval=-8.8818e-016,求解多项式,例,4.4.2非线性方程组求解,x向量；F(x)函数向量。,x：返回的向量解;fval=F(x)：函数值向量；Jacobian：解x处的Jacobian阵；fun：用于定义需求解函数;x0：初始估计值，列向量；option：最优化工具箱的选项设定,方程组的标准形式：,F(x)=0,x,fval,exitflag,output,jacobian=fsolve(fun,x0,options),functionF=myfun(x)F=2*x(1)-x(2)-exp(-x(1);-x(1)+2*x(2)-exp(-x(2);end,x,fval=fsolve(myfun,-5;-5),求方程组解,例,FunctionF=myfun(x)。%x为自变量所构成的数组。F=表达式1；表达式2；表达式m,在x1=-5,x2=-5解,functionF=myfun2(x)F=zeros(2,1)F(1)=2*x(1)-x(2)-exp(-x(1);F(2)=-x(1)+2*x(2)-exp(-x(2);end,x,fval=fsolve(myfun2,-5;-5),FunctionF=myfun(x)。%x为m个自变量所构成的数组。m个F表达式：F（1）=表达式1；F（2）=表达式2；F（m）=表达式m,Function函数的格式,列向量表达式,(1)建立函数文件myfun.mfunctiony=myfun(x)y(1)=x(1)-0.6*sin(x(1)-0.3*cos(x(2);y(2)=x(2)-0.6*cos(x(1)+0.3*sin(x(2);,x,fval=fsolve(myfun,0.5,0.5,optimset(Display,iter),在(0.5,0.5)附近解。,例,(2)初值x0=0.5,y0=0.5下，调用fsolve求方程的根。x=fsolve(myfun,0.5,0.5,optimset(display,iter)x=0.63540.3734,将求得的解代回原方程，检验结果是否正确，可见得到了较高精度的结果。,q=myfun(0.6354,0.3734)q=1.0e-004*-0.2744-0.6564,例,4.5常微分方程,微分方程：y=f(t,y),t0ttn初始条件：y(t0)=y0,1.一阶常微分方程：,方程的初值问题数值解：,求解y(t)在节点t0,t1,t2,t3.tn,处近似值y0,y1,y2，y3.yn；采用等距节点tn=t0+nh，n=0,1,.n,h是步长,微分方程描述：,(ODE),t：时间列向量;y:对应t的数值解(列向量);odefun：显式方程y=f(t,y),或含混合矩阵方程M(t,y),tspan：时间t范围，形式t0,tf,或t0,tf,y0：初始条件的值,options：用odeset设置可选参数,龙格库塔法,t,y=ode23(odefun,tspan,y0,options)t,y=ode45(odefun,tspan,y0,options),一阶常微分方程数值解,设有初值问题，y=(y2-t-2)/4/(t+1);0t1;y0(t=0)=2试求其数值解，并与精确解(y(t)=sqrt(t+1)+1)比较。(1)建立函数文件funt.m。functionyp=funt(t,y)yp=(y2-t-2)/4/(t+1);end(2)求解微分方程。t0=0;tf=1;y0=2;t,y=ode23(funt,t0,tf,y0);%求数值解y1=sqrt(t+1)+1;%求精确解tyy1y为数值解，y1为精确值，显然两者近似。,t=linspace(0,1,11),例,t,y=ode45(t,y)(y2-t-2)/4/(t+1),0,1,2),求解著名的Rossler微分方程组a=b=0.2,c=5.7,x(0)=y(0)=z(0)=0,functiondw=rosswc(t,w)dw=-w(2)-w(3);w(1)+0.2*w(2);0.2+(w(1)-5.7)*w(3);end,x0=0,0,0;t,y=ode45(rosswc,0,100,x0)plot(t,y)figureplot3(y(:,1),y(:,2),y(:,3),例,w=x,y,zdw=dx/dt;dy/dt;dz/dt,一阶常微分方程组数值解,t0,100,解法一,t,y=ode45(t,w)-w(2)-w(3);w(1)+0.2*w(2);0.2+(w(1)-5.7)*w(3),0,100,0,0,0),functiondw=rossf(t,w)dw=zeros(3,1)%构建一阶导数列向量dw(1)=-w(2)-w(3);dw(2)=w(1)+0.2*w(2);dw(3)=0.2+(w(1)-5.7)*w(3);end,x0=0,0,0;t,y=ode45(rossf,0,100,x0)plot(t,y)figureplot3(y(:,1),y(:,2),y(:,3),解法二,FunctionF=myfun(x)。%x为m个自变量所构成的数组。F=zeros(m,1)；m个F表达式：F（1）=表达式1；F（2）=表达式2；F（m）=表达式m,求解著名的Rossler微分方程组a=b=0.2,c=5.7,x(0)=y(0)=z(0)=0,functiondw=rossw(t,w,flag,a,b,c)dw=-w(2)-w(3);w(1)+a*w(2);b+(w(1)-c)*w(3);end,x0=0,0,0;t,y=ode45(rossw,0,100,x0,0.2,0.2,5.7)plot(t,y)figureplot3(y(:,1),y(:,2),y(:,3),解法三,w=x,y,zdw=dx/dt;dy/dt;dz/dt,一阶常微分方程组数值解,例,例,2.二阶常微分方程：,dx2/dt2=f(t,x,dx/dt),22(12)+=0,求解振荡器的经典Verderpol的微分方程,令y(1)=x,y(2)=dx/dt,y=y(1),y(2),y=y(1),y(2),一阶y的微分方程组:,初始条件：,x(t=0)=1，y(1)(t=0)=1x(t=0)=0,y(2)(t=0)=0,Y0=1;0,dx/dt=,y(2),dx2/dt2=,mu*(1-,*y(2),+y(1),y(1)2),dy(1)/dt=,dy(2)/dt=,globalMUMU=1t,y=ode23(verderpol1,0,40,1;0);plot(t,y),functionyy=verderpol1(t,y)globalMUyy=y(2);MU*(1-y(1).2).*y(2)-y(1);,时间0,40,s=dsolve(D2x-(1-x2)*Dx+x=0,x(0)=1,Dx(0)=0,t),22(12)+=0,Indsolveat197s=emptysym,functionyy=ver(t,y)yy=y(2);1000*(1-y(1).2).*y(2)-y(1);,t,y=ode23(ver,0,40,2;0);plot(t,y),已知微分方程,采用数值运算ode和符号运算dsolve求解,s=dsolve(D2x-1000*(1-x2)*Dx-x=0,x(0)=2,Dx(0)=0),Indsolveat197s=emptysym,t,y=ode45(t,y)y(2);1000*(1-y(1).2).*y(2)-y(1),0,40,2;0),3.高于2阶常微分方程的求解基本过程,（1）规律、定律、公式的描述形式：,初始条件：,微分方程：,（2）,高阶方程(组)转换一阶：,一阶微分方程组：,初始条件：,列向量,（3）根据(1)与(2)，编写导数M函数文件；（4）将M文件与初始条件传递给求解器Solver;（5）运行后得到ODE的、在指定时间区间解列向量y（其中包含y及不同阶的导数）。,3.显式常微分方程组解法对比,t,y=solver(odefun,tspan,y0),solver,solver,求解器solver,奇异矩阵,奇异矩阵,4.求解器Solver与方程组的关系,5.不同求解器Solver的特点,4.6数值积分和微分,4.6.1数值积分4.6.2数值微分,将区间a,b等分n个子区间xi,xi+1，i=1n，x1=a，xn+1=b。求定积分就求和问题。数值积分方法：简单的梯形法trapz辛普生(Simpson)法牛顿柯特斯(Newton-Cotes)法,4.6.1数值积分,1.数值定积分基本原理,向量自变量X和应变量Y,(1)梯形法数值积分trapz,用trapz函数计算定积分。X=1:0.01:2.5;Y=exp(-X);%生成函数关系数据向量trapz(X,Y)ans=0.28579682416393,2.数值积分的实现方法,z=trapz(Y)z=28.5797,Z=trapz(Y),Z=trapz(X,Y),y向量和,q,n=quad(exp(-x),1,2.5)q=0.2858n=13,例,fun：被积函数；n：被积函数调用次数；a和b：定积分下限和上限；tol：控制积分精度，缺省1e-6；trace,q,n=quad(fun,a,b,tol,trace),q=quad(fun,a,b,tol,trace),(2).变步长辛普生法（自适应Simpleson）,非0展现积分过程；0不展现，缺省时trace=0；返回参数I即定积分值,是否展现积分过程,P为压力kpa，V为体积，m3；n为摩尔数kmol；R为理想气体常数8.314kpam3/kmolK;T为温度。气缸内1mol气体，温度为300k，气体在整个过程恒温，体积由V1=1m3膨胀到V2=5m3,求解气缸活塞做功,n=input(n=);T=input(T=);R=8.314;pp=(v)n*R*T./v;W=quad(pp,1,5),例,(1)建立被积函数fesin.m。functionf=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);(2)调用数值积分函数quadS,n=quad(fesin,0,3*pi)S=0.9008n=77,求定积分,quad(-0.5*x).*sin(x+pi/6),0,3*pi),q,n=quad(sqrt(1+cos(x).2),0,pi/2,1.0e-6,1),例,s,n=quad(0.2+sin(x),0,pi/2,1e-8);%辛普生法求面积x=0:pi/10:pi/2;y=0.2+sin(x);s1=sum(y);ss=s1*pi/10;%求和求面积st=trapz(x,y);%求梯形面积formatshortdisp(quad积分,blanks(4),sum求积分,blanks(4),trapz求积分)disp(s,ss,st),quad积分sum求积分trapz求积分1.31421.35781.3059,例,参数含义和quad相似；tol缺省值10-6；调用步数明显小于quad，高效率求值；积分精度更高。,I,n=quadl(fun,a,b,tol,trace),(3)牛顿柯特斯法,求定积分,(1)被积函数文件fx.m。functionf=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x);(2)调用函数quadl求定积分。I=quadl(fx,0,pi)I=2.4674,例,分别用quad和quadl求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。,调用函数quad求定积分：formatlong;fx=inline(exp(-x);I,n=quad(fx,1,2.5,1e-10)I=0.28579444254766n=65,调用函数quadl求定积分：formatlong;fx=inline(exp(-x);I,n=quadl(fx,1,2.5,1e-10)I=0.28579444254754n=33,例,在xmin,xmaxymin,ymax区域二重定积分；参数tol，trace的用法与函数quad相同；tol或method可以忽略。,3.二重定积分的数值求解,q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method),(1)建立函数文件fxy.m：functionf=fxy(x,y)globalki;ki=ki+1;%ki用于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);(2)调用dblquad函数globalki;ki=0;I=dblquad(fxy,-2,2,-1,1)kiI=1.57449318974494ki=1038,计算二重积分,f=inline(exp(-x.2/2).*sin(x.2+y),x,y);dblquad(f,-2,2,-1,1)ans=1.5745dblquad(f,-2,2,-1,1)ans=1.5745,dblquad(exp(-x.2/2).*sin(x.2+y),-2,2,-1,1),例,4.三重积分数值求解,triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method),例计算三重积分,triplequad(x,y,z)y.*sin(x)+z.*cos(x),0,pi,0,1,-1,1)ans=1.999999994362637,triplequad(inline(y.*sin(x)+z.*cos(x),0,pi,0,1,-1,1)ans=1.999999994362637,triplequad(y.*sin(x)+z.*cos(x),0,pi,0,1,-1,1),tol或method可以忽略,例,4.7.2数值微分,向前差分f(x)=f(X+h)-f(x)向后差分f(x)=f(X)-f(x-h)中心差分f(x)=f(X+h/2)-f(x-h/2),h充分小,差分,1.数值差分与差商,差商,导数,2.数值微分函数diff计算差分,向量X向前差分：,Y=diff(X),向量X的n阶向前差分：,Y=diff(X,n),矩阵A的n阶差分：,dX(i)=X(i+1)-X(i)，i=1,2,n-1。,diff(X,2)=diff(diff(X),Y=diff(A,n,dim),dim=1时(缺省状态)，按行计算；dim=2按列计算。,差分运算。,例,a=magic(3)a=816357492ad=diff(a)%一阶行差分ad=-54114-5adl=diff(