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 符串

2、形式符串形式,x多项式变量 y=poly2str(P, x) (2) 字符串表示法：字符串表示法： y= poly2sym(P, x) (3) 符号多项式表示法：符号多项式表示法： 完整形式完整形式 p=1,3,0,4,5 y=poly2str(p,x) y = x4 + 3 x3 + 4 x + 5 syms x y=poly2sym(p, x) y=x4 + 3*x3 + 4*x + 5 p=1,3,0,4,5 p = 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.

3、1980i 0.5594 - 1.1980i -0.8843 例例 2.多项式的运算函数多项式的运算函数 多项式的值；根和微分；拟合曲线；部分分式多项式的值；根和微分；拟合曲线；部分分式 polyval polyval polyvalm covolution 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的

4、根的根 P=1,8,3,4,-10; r=roots(P) 多项式根的r向量： 创建多项式创建多项式 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,4 a =1 2 3 4 y=polyval(p

5、,a) y =6 22 56 114 y=polyval(p,2) y = 22 例例 例例 Y = polyvalm(P,X) 矩阵多项式求值矩阵多项式求值 Y=a0*Xn+a1*Xn-1an+1 p X必为方阵，作必为方阵，作自变量自变量代入多项式求值；代入多项式求值； p 结果阶数还是保持方阵阶数。结果阶数还是保持方阵阶数。 X 设A为方阵， P代表多项式x3-5x2-2： pp=polyval(P,A)的含义 pp=A.*A.*A -5*A. *A-2 Pm=polyvalm(P,A)的含义： Pm=A*A*A-5*A*A-2 p=1,-5,0,-2; a=2,4;1,0 a = 2

6、4 1 0 pp=polyval(p,a) pp = -14 -18 -6 -2 pm=polyvalm(p,a) pm =-18 -8 -2 -14 例例 例例 14 p 相同次数多项式，加减其系数向量相同次数多项式，加减其系数向量, , 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,-1

7、0; p2=2,-1,3; p12=conv(p1,p2) p12 = Columns 1 through 5 2 15 -5 24 -20 Columns 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; 卷积是两个变量在某范围内相乘后

8、求和的结果。卷积的变 量是序列x(n)和h(n)，y=conv(x,h) 来实现卷级的，输出的来实现卷级的，输出的 结果个数等于结果个数等于x x的长度与的长度与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=d

9、econv(p12,p2) q = 1 8 0 0 -10 r = Columns 1 through 5 0 0 0 0 0 Columns 6 through 7 0 0 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) 例例 例例 21 I=polyint(2,-1,0,3,5); 例：例： p(x) = 2x3 - x2 + 3求求 常数项常数项 5 (

10、) dp xx 不定积分，常数项取不定积分，常数项取 c 不定积分，常数项不定积分，常数项取取 零零 例例 曲线拟合和数据插值曲线拟合和数据插值 4.2.1 曲线拟合曲线拟合 4.2.2 数据插值数据插值 p 最小二乘法(误差平方和最小) ； 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：标准偏差、残

11、差范数，normr=norm(y-V*p)， 此处的p为求解之后的数值。 年份年份1900191019201930194019501960197019801990 人口人口75.9991.97105.71123.20131.66150.69179.32203.21226.50249.63 对不同年份人口数据分别进行1次、2次和9次多项式拟合 （polyfit），用poly2str表示多项式完整形式法； 1次、2次和9次多项式估计2000年人口(polyval)，结合美国 2000年人口普查截至2000年4月1日美国人口2.81421906亿 数据 绘制1900 2000年间的时间-人口数(po

12、lyval)曲线，要求用 plot不同线型（LineSpec）绘制原始数据点、拟合的1次、 2次和9次多项式曲线，标注图例 ，说明是否阶数越到高 越好 美国从1900年1990年每隔10年进行人口普查的数据如下表所示：（百万） 例例 线型线型说明说明数据点数据点标记符标记符说明说明 颜色颜色 说明说明 LineSpec 例如r-.*、-.r*、*-.r等形式是等效的 plot(x,y,-gh) clear n=1900:10:1990; r=7599,9197,10571,12320,13166,15069,17932,20321,22650,24963; nrf1=polyfit(n,r,1

13、) %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;

14、% 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); 例例 hold on; plot(x,y, *); plot(xx,y1,r-); plot(xx,y2,g-); plot(xx,y3,b-)

15、 多项式次数要适当，多项式次数要适当， 过低误差大，过高波过低误差大，过高波 动明显动明显 已知：10个实验数据，分别采用2次、5次、9次和11次拟合，选出最佳拟合次数 3. 3. 函数拟合函数拟合 clear x=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) 例例 * *eb xya ln( )ln( )*yab

16、xly 2 2 1/() 1/ yaxbxc dyyaxbxc 例例 插值构造的函数必须通过已知数据点； 拟合则不要求,只要均方差最小。 都需根据已知数据构造函数； 可使用得到函数计算未知点的函数值。 4.2.2 数据插值数据插值 p构造函数近似表达式的方法。 p常用多项式作插值函数,称多项式插值。 p多项式插值基本思想： 高次多项式或分段的低次多项式为被插值函数f(x)的近似解析表达式。 1. 插值法插值法 拉格朗日插值法 牛顿插值法 埃尔米特插值法 分段插值法 样条插值法等 函数函数y=f(x) 一维插值原理：一维插值原理： 调用格式调用格式： yi = interp1(x,y,xi,me

17、thod) p 已知X,Y两个等长向量，采样点和样本值； p Xi是向量或标量，欲插值点； p Yi是一个与Xi等长的插值结果。 采用差值法估计美国的人口数量 (1)2000年人口及 1900 1996年每隔4的人口数(interp1)； (2)将上题原始数据点、2阶拟合曲线和插值曲线绘制在同 一图，标注坐标轴为时间和人口（xlabel、 ylabel）、图形标题（title）为插值和拟合和图例 legend 年份年份1900191019201930194019501960197019801990 人口人口75.9991.97105.71123.20131.66150.69179.32203.

18、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,n4) nr4=interp1(n,r,n4,spline) % 插值1900 1996年每隔4的人口数 nr2000=interp1(n,r,2000) % 线性插值2000年人口数 nr2000=interp1(n,r,

19、2000,spline) % 插值2000年人口数 plot(n,r,or,n4,nrfv2,-*k, n4,nr4,-dg) %绘图 legend(原始数据,2阶多项式拟合曲线,插值曲线) %图例 例例 插值点样本值落在两相邻数据点的连线上. (缺省). 两点间插值点对应值就是离两点最近那点值。 已知数据求3次多项式，用多项式求插值。 每分段内构造一3次多项式，插值函数除满足插值条件和节点处光滑条 件，再按照样条函数插值。 yy = spline(x,Y,xx) clear x=0:1:10; y=sin(x); xx=0:0.2:10; yy=interp1(x,y,xx) subplot

20、(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) X X1 1取值范围不能超出取值范围不能超出X X给给 定范围，否则会给出定范围，否则会给出 “NaN”“NaN”错误。错误。 例例 38.02

21、5221 38.075267 38.125318 38.175487 38.225815 38.2751502 38.3253162 38.3756207 38.4258414 38.4758156 38.5255975 38.5753473 38.6251601 38.675710 38.725386 38.775325 38.825247 38.875232 38.925207 38.975182 39.025157 例例 3838.238.438.638.83939.239.4 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 plot(xrd

22、1(:,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 38.8750 38.9250 38.9750 39.0250 y=xy(:,2) y = 221 2

23、67 318 487 815 1502 3162 6207 8414 8156 5975 3473 1601 710 386 325 247 232 207 182 157 3838.138.238.338.438.538.638.738.838.939 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 yys=spline(xrd1(:,1),xrd1(:,2),xx) 函数函数z=f(x,y) 二维插值的原理：二维插值的原理： p 向量的采样点向量的采样点X,Y，与采样点对应函数值，与采样点对应函数值Z； p 向量或标量的欲插值点向量或标量的欲插

24、值点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,y,z,xi,yi,method) 输入参量（XI,YI）

25、通常是规则的格点 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、程序退出标志、程序退出标志exitflag 输出结果输出结

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

27、und in the interval -28, 190.51 l 配置优化参数, 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

28、，x2：范围； p option：的选项设定 x = 2.5148 f = -0.4993 e = 1 o = 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) end x,f,e,o=fminbnd(mmfun,-10,10,optims

29、et(display,iter) 0.12 sin( )0.5 (0.1) sin( ) x yexxx 求 在（-10，10） 最小值 例例 x,f,e,o=fminbnd(mmfun,6,10,optimset(display,iter) x = 8.0236 f = -3.5680 e = 1 o = iterations: 9 funcCount: 10 algorithm: golden section search, parabolic interpolation message: 优化已终止: 当前的 x 满足使用 1.000000e-04 的 OPTIONS.Tol. 黄金分割

30、搜索，抛物线插值 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) 最小最小 最小最小 0.12 sin( )0.5 (0.1) sin( ) x yexxx 求 在（6，10）最小值 -10-8-6-4-20246810 -4 -3 -2 -1 0 1 2 3 4 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)*

31、(sin(x)2)-0.5*(x+0.1)*sin(x),-10,10) x = 2.514797840754235 fval = -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.830203748934343 fval = -3.947274022275747 -10-8-6

32、-4-20246810 -4 -3 -2 -1 0 1 2 3 4 最大值求解：最大值求解：-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) 最大最大 例例 p在初始x0附近寻找局部最小值； p使用options选项来指定优化器的参数； 著名著名Rosenbrocks Banana 测试函数测试函数,理论极小值理论极小值x=1,y=1.求求极小值点极小值点 x = 1 1 fval = 0 exitflag =1 iterations: 2

33、4 funcCount: 49 algorithm: Nelder-Mead simplex direct search message: 优化已终止: 当前的 x 满足使用 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).2 end x,fval,exitflag,output=fminsearch(100*(x(2)-x(1).2).2+(1-x(1).2,1;1) 代数方程代数方

34、程 微分方程微分方程 线性方程线性方程 非线性方程非线性方程 常微分方程常微分方程 偏微分方程偏微分方程 4.4 非线性方程数值求解非线性方程数值求解 4.4.1 单变量非线性方程求解单变量非线性方程求解 () 4.4.2 非线性方程组的求解非线性方程组的求解() 迭代解法迭代解法 (fzero) x,fval,exitflag,output= fzero(fun, x1,x2, options) x,fval,exitflag,output=fzero(fun, x0, options) 2.调用格式：调用格式： 函数是否穿越横轴决定零点数值解，求方程函数是否穿越横轴决定零点数值解，求方程f

35、(x)=0根根 1.解方程原理：解方程原理： 3. 求根函数求根函数fun【f(x)】的的调用调用 求f(x)=x-10 x+2=0在x0=0.5附近的根。 例例 最优化的结构体最优化的结构体output，最优化取下列域：，最优化取下列域： algorithm 使用算法使用算法 funcCount 函数个数评估函数个数评估 interval iterations 发现区间的迭代次数发现区间的迭代次数 iterations 发现零值点迭代次数发现零值点迭代次数 message 退出信息退出信息 x,f,e,o=fzero(fun,0.5,optimset(display,iter) x =0.3

36、758；f =0；e = 1 o = interval iterations: 8 iterations: 5 funcCount: 21 algorithm: bisection, interpolation二分法 message: 在区间 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.0012 x,fval=fzero(fun,-1) x = -1.989761447718557 fval =

37、 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 fval=0 x,y=fzero(fun,-3,0.9) x = -1.9898 y = 0 例例 x,f,e,o=fzero(fun,8,optimset(display,iter) x = NaN fval = NaN exitflag =-3 正在退出 fzero: 终止搜索包含符号变化的区间 因为在搜索期间遇到 NaN 或 Inf 函数值。

38、请检查函数或使用其他起始值重试。 o = intervaliter: 22 iterations: 0 funcCount: 44 algorithm: bisection, interpolation message: 正在退出 fzero: 终止搜索包含符号变化的区间 因为在搜索期间遇到 NaN 或 Inf 函数值。 例例 无根为Nan, （Function handle) x,y=fzero(fun,0.9) x = 0.3758 y = 2.2204e-016 x,fval=fzero(x)x-10 x+2,-1,6) x =0.3758 fval = 2.2204e-016 clas

39、s(f) function_handle 待求最小值或根函数f(x)：函数名(少用)；字符串表达式，匿 名对象、函数句柄、内联函数； x,y,e,o=fzero(x-10.x+2,-1,6,optimset(Display,iter) x =-1.9898 y =0 e =1 o = interval iterations: 12 iterations: 6 funcCount: 31 algorithm: bisection, interpolation message: Zero found in the interval 0.28, -2.28 例例例例 inline 本质和本质和fun

40、ction一样，它直接内嵌在命令一样，它直接内嵌在命令 行里，不用单独定义行里，不用单独定义function。 函数表达式函数表达式 inlineinline(expression(expression ) ) f=inline(x-10.x+2); f(3) ans = -995 f=inline(x-10.x+2,x) Z=fzero(f,0.5) 或或z=fzero(inline(x-10.x+2),0.5) hd = fun fa =(x)x-10 x+2 fs =x-10 x+2 fi=inline(x-10 x+2) Name Size Bytes Class fa 1x1 16

41、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 x,fval=fzero(x3-2*x-5,3) x = 2.0946 fval = -8.8818e-016 例例 4.4.2 非线性方程组求解非线性方程组求解 x向量；向量； F(x) 函数向量函数向量。 p x：返回的向量向量解; p fval=F(x)：函数值向

42、量；函数值向量； p Jacobian：解x处的Jacobian阵； p fun：用于定义需求解函数; p x0：初始估计值，列向量列向量； p option：的选项设定 方程组的标准形式：方程组的标准形式：F(x) = 0 x,fval,exitflag,output,jacobian=fsolve(fun,x0,options) function F = 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) 求方程组解求方程组解例例 Function F

43、= myfun (x)。 % x为自变量所构成的数组。为自变量所构成的数组。 F=表达式表达式1；表达式；表达式2；表达式表达式m 1 2 2 120 12 20 x x xxe xxe 在在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(2)=x(2)-0.6*cos(x(1)+0.3*sin(x(2); 1 0.6 sin( 1) 0.3 cos( 2)0 2 0.6 cos( 1) 0.3 sin( 2)0 xxx xxx x,fval=fso

44、lve(myfun,0.5,0.5,optimset(Dis play,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 将求得的解代回原方程，检验结果是否正确，将求得的解代回原方程，检验结果是否正确， 可见得到了较高精度的结果。可见得到了较高精度的结果。 q=myfun(0.6354,0.3734) q = 1.0e-004 * -0.2744 -0.6564 例例

45、 4.5 常微分方程常微分方程 微分方程：y=f(t,y), t0ttn 初始条件：y( t0 )= y0 1. 一阶常微分方程一阶常微分方程 ： 方程的初值问题数值解： 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),或含混合矩阵方程M(t,y), 时间t范围，形式t0,tf, 或 初始条件的值, 用odeset设置可选参数 龙格库塔法龙格库塔法 t,y=ode23(odefun,tsp

46、an,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=1; y0=2; t,y=ode23(funt,t0,tf,y0); %求数值解 y1=sqrt(t+1)+1; %求精确解 t y y1 y为数值解，y1为

47、精确值，显然两者近似。 t=linspace(0,1,11) 例例 求解著名的求解著名的Rossler微分方程组微分方程组 a=b=0.2,c=5.7,x(0)=y(0)=z(0)=0 function 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,a,b,c) plot(t,y) figure plot3(y(:,1),y(:,2),y(:,3) 01020304050607

48、08090100 -15 -10 -5 0 5 10 15 20 25 -10 -5 0 5 10 15 -20 -10 0 10 0 5 10 15 20 25 ( )( )( ) ( )( )( ) ( ) ( ) ( ) x ty tz t y tx tay t z tbx tc z t 例例 例例 2. 二阶常微分方程二阶常微分方程 ： x=f(t,x,x) 求解振荡器的经典求解振荡器的经典Ver der pol的微分方程的微分方程 令y(1)=x, y(2)=d x/dt 222 2 (1) y(2) (1x )x= (1y(1) )y(2)y(1)(2) dydx dtdt dx

49、dyd x dt dtdt 一阶微分方程组一阶微分方程组: 初始条件初始条件： x(t=0)=1，y(1)(t=0)=1 x(t=0)=0, y(2)(t=0)=0 0510152025303540 -3 -2 -1 0 1 2 3 global MU MU=1 t,y=ode23(verderpol1,0,40,1;0); plot(t,y) function yy=verderpol1(t,y) global MU yy=y(2); MU*(1-y(1).2).*y(2)-y(1); 0510152025303540 -6 -4 -2 0 2 4 6 MU=3 2. 高于高于2阶常微分方程

50、阶常微分方程的的求解求解基本过程基本过程 （1）规律、定律、公式）规律、定律、公式的的描述描述形式：形式： ( ) ( ,., )0 n F y y yyt (1) 000101 ( ),( ),.( ) n n y tyy tyyty 初始条件初始条件： 微分方程微分方程： （2） (1) 12 ,., n n yy yyyy 高阶方程高阶方程(组组)转换一阶：转换一阶： ),( ),( ),( 2 1 2 1 ytf ytf ytf y y y y n n 1 1 0 0 02 01 0 )( )( )( )( nn y y y ty ty ty ty 一阶微分方程组： 初始条件： 列向量

51、列向量 (2) (1) (1) ( ) ( ) (2)=z(1) (3 )z(2) - n z(1) (2) k k n n k y y y z z zz zz zz zz yy y y yfz zz (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) = = k k kk y y yz zy z z (k+1) (k) (k+1) (k) n( -1) f,., n yy y yyt （） （3）根据(1)与(2)，编写导

52、数M函数文件； （4）将M文件与初始条件传递给求解器Solver; （5）运行后得到ODE的、在指定时间区间解列向量y （其中包含y及不同阶的导数）。 t,y=solver(odefun,tspan,y0) solversolver 求解器求解器 solver 奇异矩阵奇异矩阵 奇异矩阵奇异矩阵 函数指令含 义函 数含 义 求 解 器 S o l ver ode23 普通2-3阶 法解ODE odefile 包含ODE的 文件 ode23s 低阶法解刚 性ODE 选项 odeset 创 建 、 更 改 Solver选项 ode23t 解适度刚性 ODE odeget 读取Solver的 设置值

53、 ode23t b 低阶法解刚 性ODE 输出 odeplot ODE的时间 序列图 ode45 普通4-5阶 法解ODE odephas2 ODE的二维 相平面图 ode15s 变阶法解刚 性ODE odephas3 ODE的三维 相平面图 ode113 普通变阶法 解ODE odeprint 在命令窗口输 出结果 求 解 器 Solver ODE类型特点说明 ode45非刚性 一步算法；4，5阶Runge- Kutta方程；累计截断误差 达(x)3 大部分场合的首选算法 ode23非刚性 一步算法；2，3阶Runge- Kutta方程；累计截断误差 达(x)3 使用于精度较低的情形 ode

54、113非刚性 多步法；Adams算法；高 低精度均可到10-310-6 计算时间比ode45短 ode23t适度刚性采用梯形算法适度刚性情形 ode15s刚性 多步法；Gears反向数值 微分；精度中等 若ode45失效时，可尝试使 用 ode23s刚性 一步法；2阶Rosebrock算 法；低精度 当精度较低时，计算时间比 ode15s短 ode23tb 刚性梯形算法；低精度 当精度较低时，计算时间比 ode15s短 6. 常微分方程组解法器参数常微分方程组解法器参数 4.6 数值积分和微分数值积分和微分 4.6.1 数值积分数值积分 4.6.2 数值微分数值微分 将区间将区间a,b等分等分

55、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(-X); %生成函数关系数据向量 trapz(X,Y) ans = 0.28579682416393 2. 数值积分的实现方法数值积分的实现方法 z=tr

56、apz(Y) z =28.5797 Z = trapz(Y) Z=trapz(X,Y) q,n=quad(exp(-x),1,2.5) q = 0.2858 n = 13 例例 ：被积函数； ：被积函数调用次数； ：定积分下限和上限； ：控制积分精度，缺省1e-6； q,n=quad(fun,a,b,tol,trace) q=quad(fun,a,b,tol,trace) 非非0展现积分过程；展现积分过程； 0不展现，缺省时不展现，缺省时trace=0； 返回参数返回参数I即定积分值即定积分值 是否展现积分过程是否展现积分过程 P为压力kpa，V为体积，m3；n为摩尔数kmol； R为理想气体

57、常数8.314kpam3/kmolK;T为温 度。气缸内1mol气体，温度为300k，气体在 整个过程恒温，体积由V1=1m3膨胀到V2=5m3 (/)ln( 2/ 1)WPdVnRT V dVnRTVV PV 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).*sin(x+pi/6); (2) 调用数值积分函数调用数值积分函数quad S,n=quad(fesin

58、,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求积分求积分,blanks(4),trapz求积分求积分) disp(s,ss,st) h

59、old on plot(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).*cos(x); (2) 调用函数调用函数quad8求定积分。求定积分。 I

60、=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.28579444254766 n = 65 调用函数quad8求定积分： format long; fx=inline(exp(-x); I,n=quad8(fx,1,2.5,1e-10) I = 0.2857944425475