matlab数值运算课件

1、第四章第四章 数值计算功能数值计算功能4.1 多项式计算多项式计算4.2 数值插值和曲线拟合数值插值和曲线拟合4.3 函数极值函数极值4.4 非线性方程问题求解非线性方程问题求解4.5 常微分方程的数值求解常微分方程的数值求解4.6 数值微分与积分数值微分与积分4.7 概率统计概率统计4.1 多项式运算多项式运算 1多项式表示法多项式表示法 2多项式运算函数多项式运算函数 1多项式表示法多项式表示法 (1) 直接法直接法:多项式多项式p 多项式系数用行向量行向量表示，按降幂排列降幂排列；p 缺某幂次项，其幂次项系数为零幂次项系数为零P=a0,a1,.an-1,an符串形式符串形式,x多项式变量

2、y=poly2str(P, x)(2) 字符串表示法：字符串表示法：y= poly2sym(P, x) (3) 符号多项式表示法：符号多项式表示法：完整形式完整形式p=1,3,0,4,5y=poly2str(p,x)y = x4 + 3 x3 + 4 x + 5syms xy=poly2sym(p, x)y=x4 + 3*x3 + 4*x + 5 p=1,3,0,4,5p = 1 3 0 4 5 y=poly2str(p,x)y = x4 + 3 x3 + 4 x + 5 rp=roots(p)rp = -3.2346 0.5594 + 1.1980i 0.5594 - 1.1980i -0.

3、8843例例2.多项式的运算函数多项式的运算函数多项式的值；根和微分；拟合曲线；部分分式多项式的值；根和微分；拟合曲线；部分分式polyvalpolyvalpolyvalmcovolution n次多项式次多项式n个根，或个根，或实根实根，或若干对，或若干对共轭复根共轭复根。p P：多项式系数向量；：多项式系数向量；p r：n个根个根的的。 r=roots(P)r = -7.6998 -0.5572 + 1.1335i -0.5572 - 1.1335i 0.8141 例例求多项式求多项式x4+8x3+3x2+4x-10的根的根P=1,8,3,4,-10;r=roots(P)多项式根的r向量：

4、 创建多项式创建多项式P=poly(r)PP = 1.0000 8.0000 3.0000 4.0000 -10.0000避免浮点显示 p=poly(r)p = 1 8 3 4 -10多项式的向量系数：，则求多项式在该点，则求多项式在该点 的值；的值；每元素每元素x(i i,j)多项式求值。多项式求值。Y=a0*X.n+a1*X.(n-1)an+1 p=1,3,0,2; poly2str(p,x)ans = x3 + 3 x2 + 2 a=1,2;3,4a =1 2 3 4 y=polyval(p,a)y =6 22 56 114 y=polyval(p,2)y = 22例例例例 Y = po

5、lyvalm(P,X)矩阵多项式求值矩阵多项式求值Y=a0*Xn+a1*Xn-1an+1p X必为方阵，作必为方阵，作自变量自变量代入多项式求值；代入多项式求值；p 结果阶数还是保持方阵阶数。结果阶数还是保持方阵阶数。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 = 2 4 1 0 pp=polyval(p,a)pp = -14 -18 -6 -2 pm=polyvalm(p,a)pm =-18 -

6、8 -2 -14 例例例例14p 相同次数多项式，加减其系数向量相同次数多项式，加减其系数向量, , p 不同次数不同次数, ,低次多项式中低次多项式中。 p2 = 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例例例例多项式x4+8x3-10与2x2-x+3的乘积。 p1=1,8,0,0,-10; p2=2,-1,3; p12=conv(p1,p2)p12 = Columns 1 through 5 2 15 -5 24 -20 Col

7、umns 6 through 7 10 -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(Time index n); %标坐标轴 ylabel(Amplitude); title(Output Obtained by Convolution)%图形标题grid;卷积是两个变量在某范围内相乘后求和的结果。卷积的变量是序列x(n)和h(n)，y=conv(x,h) 来实现卷级的，输出的来实现卷级的，输出的结果个数等于结果个数等于x x的长度与的

8、长度与h h的长度之和减去的长度之和减去1 1。卷积公式：卷积公式：z(n)=x(n)z(n)=x(n)* *y(n)= x(m)y(n-m)dm.y(n)= x(m)y(n-m)dm.例例p Q：多项式P1除以P2的商式，p r：P1除以P2的余式，p Q和r仍是多项式系数向量。积这两类问题都称作解卷探、石油勘探等问题。地震信号处理、地质勘求多项式求多项式x4+8x3-10除以除以2x2-x+3的结果。的结果。 q,r=deconv(p12,p2)q = 1 8 0 0 -10 r = Columns 1 through 5 0 0 0 0 0 Columns 6 through 7 0 0

9、 q,r=deconv(p2,p12)q = 0 r = 2 -1 3 例例 多项式乘积PQ的导函数：多项式P的导函数：导函数的分子存入p，分母存入q。多项式除法P/Q的导函数：例：求有理分式的导数。命令如下：P=1;Q=1,0,5;p,q=polyder(P,Q)例例例例21I=polyint(2,-1,0,3,5); 例：例： p(x) = 2x3 - x2 + 3求求 常数项常数项 5( ) dp xx 不定积分，常数项取不定积分，常数项取 c不定积分，常数项不定积分，常数项取取 零零例例曲线拟合和数据插值曲线拟合和数据插值4.2.1 曲线拟合曲线拟合 4.2.2 数据插值数据插值 p

10、最小二乘法(误差平方和最小) ；p x、y已知数据的横、纵坐标；p n为多项式次数;N次次p S结构体数组(struct)，估计预测误差，含R，df和normr。p R：先输入x构建范德蒙矩阵V，后QR分解，得上三角矩阵。p df：自由度，df=length(y)-(n+1)。df0时，超定方程组求解，拟合点数比未知数(p(1)p(n+1)多。p normr：标准偏差、残差范数，normr=norm(y-V*p)， 此处的p为求解之后的数值。 年份年份1900191019201930194019501960197019801990人口人口75.9991.97105.71123.20131.66

11、150.69179.32203.21226.50249.63 对不同年份人口数据分别进行1次、2次和9次多项式拟合（polyfit），用poly2str表示多项式完整形式法； 1次、2次和9次多项式估计2000年人口(polyval)，结合美国2000年人口普查截至2000年4月1日美国人口2.81421906亿数据 绘制1900 2000年间的时间-人口数(polyval)曲线，要求用plot不同线型（LineSpec）绘制原始数据点、拟合的1次、2次和9次多项式曲线，标注图例 ，说明是否阶数越到高越好 美国从1900年1990年每隔10年进行人口普查的数据如下表所示：（百万）例例线型线型说

12、明说明数据点数据点标记符标记符说明说明颜色颜色说明说明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,

13、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,n

14、20); % 从1900年到2000年间9次拟合人口数 plot(n,r,or, n20,nrfv1,- p4=polyfit(x,y,11) y4=polyval(p4,xx);例例hold on;plot(x,y, *); plot(xx,y1,r-);plot(xx,y2,g-);plot(xx,y3,b-)多项式次数要适当，多项式次数要适当，过低误差大，过高波过低误差大，过高波动明显动明显已知：10个实验数据，分别采用2次、5次、9次和11次拟合，选出最佳拟合次数3. 3. 函数拟合函数拟合 clear x=0:5; y=0.2097,0.3523,0.4339,0.5236,0.75

15、90,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)例例*eb xyaln( )ln( )*yab xly 221/()1/yaxbxcdyyaxbxc例例 插值构造的函数必须通过已知数据点； 拟合则不要求,只要均方差最小。 都需根据已知数据构造函数； 可使用得到函数计算未知点的函数值。4.2.2 数据插值数据插值p构造函数近似表达式的方法。p常用多项式作插值函数,称多项式插值。p多项式插值基本思想：高次多项式或

16、分段的低次多项式为被插值函数f(x)的近似解析表达式。1. 插值法插值法拉格朗日插值法 牛顿插值法 埃尔米特插值法分段插值法样条插值法等 函数函数y=f(x) 一维插值原理：一维插值原理： 调用格式调用格式： yi = interp1(x,y,xi,method)p 已知X,Y两个等长向量，采样点和样本值；p Xi是向量或标量，欲插值点；p Yi是一个与Xi等长的插值结果。采用差值法估计美国的人口数量(1)2000年人口及 1900 1996年每隔4的人口数(interp1)；(2)将上题原始数据点、2阶拟合曲线和插值曲线绘制在同一图，标注坐标轴为时间和人口（xlabel、ylabel）、图形

17、标题（title）为插值和拟合和图例legend年份年份1900191019201930194019501960197019801990人口人口75.9991.97105.71123.20131.66150.69179.32203.21226.50249.63例例 clear n=1900:10:1990; r=7599,9197,10571,12320,13166,15069,17932,20321,22650,24963;n4=1900:4:1996; % 1900 1996年每隔4的线性等分数组nrf2=polyfit(n,r,2) %2次多项式拟合nrfv2=polyval(nrf2,

18、n4)nr4=interp1(n,r,n4,spline) % 插值1900 1996年每隔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阶多项式拟合曲线,插值曲线) %图例 例例 插值点样本值落在两相邻数据点的连线上. (缺省). 两点间插值点对应值就是离两点最近那点值。 已知数据求3次多项式，用多项式求插值。 每分段内构造一3次多项式，插值函数除满足插值条

19、件和节点处光滑条件，再按照样条函数插值。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

20、(x,y,-*,xx,y4,dr)X X1 1取值范围不能超出取值范围不能超出X X给给定范围，否则会给出定范围，否则会给出“NaN”“NaN”错误。错误。例例38.02522138.07526738.12531838.17548738.22581538.275150238.325316238.375620738.425841438.475815638.525597538.575347338.625160138.67571038.72538638.77532538.82524738.87523238.92520738.97518239.025157例例3838.238.438.638.8393

21、9.239.40100020003000400050006000700080009000plot(xrd1(:,1),xrd1(:,2) xx=linspace(38,39,80);yy=interp1(xrd1(:,1),xrd1(:,2),xx,spline); plot(xx,yy,-*g) x=xy(:,1)x = 38.0250 38.0750 38.1250 38.1750 38.2250 38.2750 38.3250 38.3750 38.4250 38.4750 38.5250 38.5750 38.6250 38.6750 38.7250 38.7750 38.8250 3

22、8.8750 38.9250 38.9750 39.0250 y=xy(:,2)y = 221 267 318 487 815 1502 3162 6207 8414 8156 5975 3473 1601 710 386 325 247 232 207 182 1573838.138.238.338.438.538.638.738.838.9390100020003000400050006000700080009000 yys=spline(xrd1(:,1),xrd1(:,2),xx) 函数函数z=f(x,y) 二维插值的原理：二维插值的原理：p 向量的采样点向量的采样点X,Y，与采样点对

23、应函数值，与采样点对应函数值Z；p 向量或标量的欲插值点向量或标量的欲插值点Xi,Yi，插值结果，插值结果Zi。p Xi,Yi取值不超取值不超X,Y给定范围，否则给定范围，否则“NaN”错误。错误。 zi = interp2(x,y,z,xi,yi,method)例例 运行结果如下图所示。p 向量的采样点向量的采样点X,Y，与采样点对应函数值，与采样点对应函数值Z；p 向量或标量的欲插值点向量或标量的欲插值点Xi,Yi，插值结果，插值结果Zi。griddata 算法p Xi,Yi取值不超取值不超X,Y给定范围，否则给定范围，否则“NaN”错误。错误。 xi,yi,zi = griddata(x

24、,y,z,xi,yi,method)输入参量（XI,YI）通常是规则的格点 clear p=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,)mesh(Xi,Yi,Zi)Xi,Yi,Zi=griddata(X,Y,Z,Xi,Yi,v4)mesh(Xi,Yi,Zi) 随机生成30*3矩阵例例4.3 函数极值函数极值调用格式一样调用格式一样最小值的解最小值的解或零点零点根根x、该点的函数值、该点的函数值fval、程序退出标志、程序退出

25、标志exitflag输出结果输出结果output、待求根函数fun、初始值x0、左右边界x1,x2x,fval,exitflag,output=)迭代解法迭代解法：最小值的x解或零点的根(自变量值)；：最小值或零点时函数值f(x)；待求最小值或根函数f(x)：函数名(少用)；字符串表达式，匿名对象、函数句柄、内联函数；程序退出标志, 1收敛到解x；选定输出结果；搜索初始值；左右边界 output = 迭代次数 iterations: 24函数评价次数 funcCount: 48算法 algorithm: bisection, interpolation对分插值退出信息 message: Zer

26、o found in the interval -28, 190.51l 配置优化参数, optimset函数定义参数;l 最优化工具箱的(20多个)选项设定;l 将option显示出来;l 改变其中某个选项;选项函数调用时中间结果显示方式l off不显示;l iter每步都显示;l final只显示最终结果。optimset(display,off)options = optimset(param1,value1.)求区间求区间x1 x2内函数内函数f(x)最小值时最小值时p X：返回最小值的解;p fval=F(x)，即最小值；p fun：用于定义需求解函数;p x1，x2：范围；p op

27、tion：的选项设定x = 2.5148f = -0.4993e = 1o = iterations: 13 funcCount: 14 algorithm: golden section search, parabolic interpolation message: 优化已终止: 当前的 x 满足使用 1.000000e-04 的 OPTIONS.Tol.function f=mmfun(x)f=exp(-0.1*x)*(sin(x)2)-0.5*(x+0.1)*sin(x)endx,f,e,o=fminbnd(mmfun,-10,10,optimset(display,iter)0.12

28、sin( )0.5 (0.1) sin( )xyexxx求在（-10，10） 最小值例例x,f,e,o=fminbnd(mmfun,6,10,optimset(display,iter)x = 8.0236f = -3.5680e = 1o = iterations: 9 funcCount: 10 algorithm: golden section search, parabolic interpolation message: 优化已终止: 当前的 x 满足使用 1.000000e-04 的 OPTIONS.Tol.黄金分割搜索，抛物线插值 x=-10:0.05:10; y=exp(-0.

29、1*x).*(sin(x).2)-0.5*(x+0.1).*sin(x);plot(x,y) plot(x,y)最小最小最小最小0.12sin( )0.5 (0.1) sin( )xyexxx求在（6，10）最小值-10-8-6-4-20246810-4-3-2-101234 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)x,fval=fminbnd(exp(-0.1*x)*(sin(x)2)-0.5*(x+0.1)*sin(x),-10,10)x = 2.514797840754

30、235fval = -0.499312445280039最小最小x,fval=fminbnd(exp(-0.1*x)*(sin(x)2)-0.5*(x+0.1)*sin(x),6,10)x = 8.0236；fval =-3.5680例例x,fval=fminbnd(-exp(-0.1*x)*(sin(x)2)+0.5*(x+0.1)*sin(x),-8,-2)x = -4.830203748934343fval = -3.947274022275747-10-8-6-4-20246810-4-3-2-101234最大值求解：最大值求解：-f(x)在区间(a,b)上的最小值 x=-10:0.0

31、5:10; y=-exp(-0.1*x).*(sin(x).2)+0.5*(x+0.1).*sin(x);plot(x,y) plot(x,y)最大最大例例p在初始x0附近寻找局部最小值；p使用options选项来指定优化器的参数；著名著名Rosenbrocks Banana 测试函数测试函数,理论极小值理论极小值x=1,y=1.求求极小值点极小值点x = 1 1fval = 0exitflag =1iterations: 24 funcCount: 49 algorithm: Nelder-Mead simplex direct search message: 优化已终止: 当前的 x 满足

32、使用 1.000000e-04 的 OPTIONS.TolX 的终止条件，.例例x,fval,exitflag,output=fminsearch(fxy,1;1)function f=fxy(x)f=100*(x(2)-x(1).2).2+(1-x(1).2endx,fval,exitflag,output=fminsearch(100*(x(2)-x(1).2).2+(1-x(1).2,1;1) 代数方程代数方程微分方程微分方程线性方程线性方程非线性方程非线性方程常微分方程常微分方程偏微分方程偏微分方程4.4 非线性方程数值求解非线性方程数值求解4.4.1 单变量非线性方程求解单变量非线性

33、方程求解 ()4.4.2 非线性方程组的求解非线性方程组的求解()迭代解法迭代解法(fzero)x,fval,exitflag,output= fzero(fun, x1,x2, options)x,fval,exitflag,output=fzero(fun, x0, options)2.调用格式：调用格式：函数是否穿越横轴决定零点数值解，求方程函数是否穿越横轴决定零点数值解，求方程f(x)=0根根 1.解方程原理：解方程原理：3. 求根函数求根函数fun【f(x)】的的调用调用求f(x)=x-10 x+2=0在x0=0.5附近的根。例例最优化的结构体最优化的结构体output，最优化取下列

34、域：，最优化取下列域： algorithm 使用算法使用算法 funcCount 函数个数评估函数个数评估 interval iterations 发现区间的迭代次数发现区间的迭代次数 iterations 发现零值点迭代次数发现零值点迭代次数 message 退出信息退出信息 x,f,e,o=fzero(fun,0.5,optimset(display,iter)x =0.3758；f =0；e = 1o = interval iterations: 8 iterations: 5 funcCount: 21 algorithm: bisection, interpolation二分法 me

35、ssage: 在区间 0.34, 0.613137 中发现零建立函数文件fun.m。function fv=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.0012x,fval=fzero(fun,-1)x = -1.989761447718557fval = 0 (1) 建立函数文件建立函数文件fun.m。function fv=fun(x)fv=x-10.x+2;(2) 调用调用fzero函数求根。函数求根。 x,fval,e,output=fzero(fun,0.5)x = 0.3758 f

36、val=0 x,y=fzero(fun,-3,0.9)x = -1.9898y = 0例例x,f,e,o=fzero(fun,8,optimset(display,iter) x = NaNfval = NaNexitflag =-3正在退出 fzero: 终止搜索包含符号变化的区间 因为在搜索期间遇到 NaN 或 Inf 函数值。请检查函数或使用其他起始值重试。o = intervaliter: 22 iterations: 0 funcCount: 44 algorithm: bisection, interpolation message: 正在退出 fzero: 终止搜索包含符号变化的

37、区间 因为在搜索期间遇到 NaN 或 Inf 函数值。例例无根为Nan,（Function handle) x,y=fzero(fun,0.9)x = 0.3758y = 2.2204e-016 x,fval=fzero(x)x-10 x+2,-1,6)x =0.3758fval = 2.2204e-016class(f)function_handle 待求最小值或根函数f(x)：函数名(少用)；字符串表达式，匿名对象、函数句柄、内联函数； x,y,e,o=fzero(x-10.x+2,-1,6,optimset(Display,iter)x =-1.9898y =0e =1o = inter

38、val iterations: 12 iterations: 6 funcCount: 31 algorithm: bisection, interpolation message: Zero found in the interval 0.28, -2.28例例例例inline本质和本质和function一样，它直接内嵌在命令一样，它直接内嵌在命令行里，不用单独定义行里，不用单独定义function。函数表达式函数表达式inlineinline(expression(expression) ) f=inline(x-10.x+2); f(3) ans = -995f=inline(x-10.

39、x+2,x)Z=fzero(f,0.5)或或z=fzero(inline(x-10.x+2),0.5)hd = funfa =(x)x-10 x+2fs =x-10 x+2fi=inline(x-10 x+2) Name Size Bytes Class fa 1x1 16 function_handle fi 1x1 832 inline fs 1x8 16 char hd 1x1 16 function_handle class(f)ans =inline例例 roots(1,0,-2,-5)ans = 2.0946 -1.0473 + 1.1359i -1.0473 - 1.1359i

40、x,fval=fzero(x3-2*x-5,3)x = 2.0946fval = -8.8818e-016例例4.4.2 非线性方程组求解非线性方程组求解x向量；向量；F(x) 函数向量函数向量。p x：返回的向量向量解;p fval=F(x)：函数值向量；函数值向量；p Jacobian：解x处的Jacobian阵；p fun：用于定义需求解函数;p x0：初始估计值，列向量列向量；p option：的选项设定方程组的标准形式：方程组的标准形式：F(x) = 0 x,fval,exitflag,output,jacobian=fsolve(fun,x0,options) function F

41、 = myfun(x)F = 2*x(1)-x(2)-exp(-x(1); -x(1)+2*x(2)-exp(-x(2);endx,fval = fsolve(myfun, -5; -5) 求方程组解求方程组解例例Function F = myfun (x)。 % x为自变量所构成的数组。为自变量所构成的数组。F=表达式表达式1；表达式；表达式2；表达式表达式m 122 12012 20 xxxxexxe在在x1=-5,x2=-5解解 (1) 建立函数文件建立函数文件myfun.m function y=myfun(x)y(1)=x(1)-0.6*sin(x(1)-0.3*cos(x(2);y

42、(2)=x(2)-0.6*cos(x(1)+0.3*sin(x(2);1 0.6 sin( 1) 0.3 cos( 2)02 0.6 cos( 1) 0.3 sin( 2)0 xxxxxxx,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.6354 0.3734将求得的解代回原方程，检验结果是否正确，将求得的解代

43、回原方程，检验结果是否正确， 可见得到了较高精度的结果。可见得到了较高精度的结果。q=myfun(0.6354,0.3734)q = 1.0e-004 * -0.2744 -0.6564例例4.5 常微分方程常微分方程微分方程：y=f(t,y), t0ttn初始条件：y( t0 )= y01. 一阶常微分方程一阶常微分方程 ： 方程的初值问题数值解：p 求解y(t)在节点t0,t1 ,t2,t3.tn, 处近似值y0, y1 , y2，y3. yn；p 采用等距节点tn=t0+nh，n=0,1,.n,h是步长是步长 微分方程描述：时间列向量;对应t的数值解(列向量);显式方程y=f(t,y),

44、或含混合矩阵方程M(t,y), 时间t范围，形式t0,tf, 或 初始条件的值,用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。function yp=funt(t,y)yp=(y2-t-2)/4/(t+1);(2) 求解微分方程。t0=0;tf=

45、1;y0=2;t,y=ode23(funt,t0,tf,y0); %求数值解y1=sqrt(t+1)+1; %求精确解tyy1 y为数值解，y1为精确值，显然两者近似。t=linspace(0,1,11)例例求解著名的求解著名的Rossler微分方程组微分方程组a=b=0.2,c=5.7,x(0)=y(0)=z(0)=0function ddx=rossler(t,x,flag,a,b,c ) ddx=-x(2)-x(3);x(1)+a*x(2);b+(x(1)-c)*x(3);end a=0.2;b=0.2;c=5.7; x0=0,0,0;t,y=ode45(rossler,0,100,x0

46、,a,b,c) plot(t,y) figure plot3(y(:,1),y(:,2),y(:,3)0102030405060708090100-15-10-50510152025-10-5051015-20-100100510152025( )( )( )( )( )( )( ) ( ) ( )x ty tz ty tx tay tz tbx tc z t 例例例例2. 二阶常微分方程二阶常微分方程 ：x=f(t,x,x)求解振荡器的经典求解振荡器的经典Ver der pol的微分方程的微分方程 令y(1)=x, y(2)=d x/dt2222(1) y(2)(1x )x= (1y(1)

47、)y(2)y(1)(2)dydxdtdtdxdyd xdtdtdt一阶微分方程组一阶微分方程组:初始条件初始条件：x(t=0)=1，y(1)(t=0)=1x(t=0)=0, y(2)(t=0)=00510152025303540-3-2-10123global MUMU=1t,y=ode23(verderpol1,0,40,1;0);plot(t,y)function yy=verderpol1(t,y)global MUyy=y(2); MU*(1-y(1).2).*y(2)-y(1);0510152025303540-6-4-20246MU=32. 高于高于2阶常微分方程阶常微分方程的的求

48、解求解基本过程基本过程 （1）规律、定律、公式）规律、定律、公式的的描述描述形式：形式：( )( ,., )0nF y y yyt(1)000101( ),( ),.( )nny tyy tyyty初始条件初始条件：微分方程微分方程：（2）(1)12,.,nnyy yyyy高阶方程高阶方程(组组)转换一阶：转换一阶：),(),(),(2121ytfytfytfyyyynn110002010)()()()(nnyyytytytyty一阶微分方程组： 初始条件： 列向量列向量 (2)(1)(1)( )( )(2)=z(1)(3)z(2)-nz(1)(2)kknnkyyyzzzzzzzzzzyyyy

49、yfzzz(k1)=(k-2)(k)=(k-1)(n)( -1)（，. (n),t）=(n).(k+1)=(k) (0)y=y=z(1)f(z(1)(2)zzz，. (n),t）= (n)( )(1)( )(1)=kkkkyyyzzyzz(k+1)(k)(k+1) (k)n( -1)f,.,nyy y yyt （）（3）根据(1)与(2)，编写导数M函数文件；（4）将M文件与初始条件传递给求解器Solver;（5）运行后得到ODE的、在指定时间区间解列向量y（其中包含y及不同阶的导数）。t,y=solver(odefun,tspan,y0)solversolver求解器求解器solver奇异矩

50、阵奇异矩阵 奇异矩阵奇异矩阵 函数指令含 义函 数含 义求解器S o lverode23普通2-3阶法解ODEodefile包含ODE的文件ode23s低阶法解刚性ODE选项odeset创 建 、 更 改Solver选项ode23t解适度刚性ODEodeget读取Solver的设置值ode23tb低阶法解刚性ODE输出odeplotODE的时间序列图ode45普通4-5阶法解ODEodephas2ODE的二维相平面图ode15s变阶法解刚性ODEodephas3ODE的三维相平面图ode113普通变阶法解ODEodeprint在命令窗口输出结果求 解 器SolverODE类型特点说明ode45

51、非刚性一步算法；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短6. 常

52、微分方程组解法器参数常微分方程组解法器参数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)法法 1.数值定积分基本原理数值定积分基本原理 向量自变量X和应变量Y(1) 梯形法数值积分梯形法数值积分用 trapz函数计算定积分。 X=1:0.01:2.5;Y=exp

53、(-X); %生成函数关系数据向量trapz(X,Y)ans = 0.285796824163932. 数值积分的实现方法数值积分的实现方法z=trapz(Y)z =28.5797Z = trapz(Y) Z=trapz(X,Y) q,n=quad(exp(-x),1,2.5)q = 0.2858n = 13例例 ：被积函数； ：被积函数调用次数；：定积分下限和上限；：控制积分精度，缺省1e-6；q,n=quad(fun,a,b,tol,trace) q=quad(fun,a,b,tol,trace) 非非0展现积分过程；展现积分过程；0不展现，缺省时不展现，缺省时trace=0；返回参数返回

54、参数I即定积分值即定积分值是否展现积分过程是否展现积分过程P为压力kpa，V为体积，m3；n为摩尔数kmol；R为理想气体常数8.314kpam3/kmolK;T为温度。气缸内1mol气体，温度为300k，气体在整个过程恒温，体积由V1=1m3膨胀到V2=5m3(/)ln( 2/ 1)WPdVnRT V dVnRTVVPV nRT求解气缸活塞做功n=input(n=);T=input(T=);R=8.314;pp=(v)n*R*T./v;W=quad(pp,1,5)例例 (1) 建立被积函数建立被积函数fesin.m。 function f=fesin(x) f=exp(-0.5*x).*si

55、n(x+pi/6); (2) 调用数值积分函数调用数值积分函数quad S,n=quad(fesin,0,3*pi) S = 0.9008 n =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=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)format short disp(quad积分积分,blanks(4),sum求

56、积分求积分,blanks(4),trapz求积分求积分)disp(s,ss,st)hold onplot(x,y,linewidth,2) stairs(x,y,-r*) stem(x,y,-gh)quad积分 sum求积分 trapz求积分 1.3142 1.3578 1.3059例例 参数含义和参数含义和quad相似；相似； tol缺省值缺省值10-6； 调用调用步数明显小于步数明显小于quad，高效率高效率求值；求值； 积分精度更高。积分精度更高。 求定积分求定积分(1) 被积函数文件被积函数文件fx.m。function f=fx(x)f=x.*sin(x)./(1+cos(x).*c

57、os(x);(2) 调用函数调用函数quad8求定积分。求定积分。I=quad8(fx,0,pi)I = 2.4674例例分别用分别用quad和和quadl求定积分的近似值，并在相同的求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。积分精度下，比较函数的调用次数。调用函数quad求定积分：format long;fx=inline(exp(-x);I,n=quad(fx,1,2.5,1e-10)I = 0.28579444254766n = 65调用函数quad8求定积分：format long;fx=inline(exp(-x);I,n=quad8(fx,1,2.5,1e-10)

58、I = 0.28579444254754n = 33例例 在在 xmin,xmax ymin,ymax 区域二重定积分；区域二重定积分； 参数参数toltol，tracetrace的用法与函数的用法与函数quadquad相同；相同； toltol或或methodmethod可以忽略。可以忽略。q = dblquad(fun,xmin,xmax,ymin,ymax,q = dblquad(fun,xmin,xmax,ymin,ymax,tol,methodtol,method) ) (1) 建立函数文件建立函数文件fxy.m：function f=fxy(x,y)global ki;ki=ki+

59、1; %ki用于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);(2) 调用调用dblquad函数函数global ki;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.5745dblquad(exp(-x.2/2).*sin(x.2+y),-2,2,-1,1)例例triplequa

60、d(fun,xmin,xmax,ymin,ymax,zmin,zmaxtriplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax, ,tol,methodtol,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)