《第4MAT六B的数值计算教学提纲》由会员分享，可在线阅读，更多相关《第4MAT六B的数值计算教学提纲（62页珍藏版）》请在金锄头文库上搜索。
1、2020/7/16,第1页,第4章 MATLAB 的数值计算,4.1 数值微积分 4.2 矩阵和代数方程 4.3 概率分布和统计分析 4.4 多项式运算和卷积,2020/7/16,第2页,4.1数值微积分,4.1.1 近似数值极限及导数 4.1.2 数值求和与近似数值积分 4.1.3 计算精度可控的数值积分 4.1.4 函数极值的数值求解 4.1.5 常微分方程的数值解,2020/7/16,第3页,在MATLAB数值计算中，既没有专门的求极限指令，也没有专门的求导指令。但MATLAB提供了与“求导”概念有关的“求差分”指令。 dx=diff(X) %计算向量X的前向差分 FX=gradient
2、(F) %求一元(函数)梯度 FX, FY =gradient(F) %求二元（函数）梯度,对diff而言，当X是向量时，dx= X(2:n)-X(1:n-1) ;当X是矩阵时，dx= X(2:n, :)-X(1:n-1, :) 。 dx的长度比x的长度少1个元素。,4.1.1 近似数值极限及导数,2020/7/16,第4页,dx=diff(X) %计算向量X的前向差分 FX=gradient(F) %求一元(函数)梯度 FX, FY =gradient(F) %求二元（函数）梯度,对gradient而言，当F是向量时，FX(1) = F(2)-F(1), FX(2:end-1) = (F(3
3、:end)-F(1:end-2)/2 , FX(end) = F(end)-F(end-1) ; FX长度与F的长度相同 当F是矩阵时， FX, FY是与F同样大小的矩阵。 FX的每行对应F相应行元素间的梯度 ; FY的每列对应F相应列元素间的梯度 ;,4.1.1 近似数值极限及导数,2020/7/16,第6页,数值极限和导数的应用应十分谨慎,%（1）自变量的增量取得过小（eps数量级） d=pi/100; t=0:d:2*pi; x=sin(t); dt=5*eps; x_eps=sin(t+dt); dxdt_eps=(x_eps-x)/dt; plot(t,x,LineWidth,5)
4、hold on plot(t,dxdt_eps) hold off legend(x(t),dx/dt) xlabel(t),数值导数受计算中有限精度困扰，当增量dt过小时，f(t+dt)与f(t)的数值十分接近，高位有效数字完全相同， df =f(t+dt)-f(t) 造成高位有效数字消失，精度急剧变差。,2020/7/16,第7页,数值极限和导数的应用应十分谨慎,%（2）自变量的增量取得适当 x_d=sin(t+d); dxdt_d=(x_d-x)/d; plot(t,x,LineWidth,5) hold on plot(t,dxdt_d) hold off legend(x(t),dx
5、/dt) xlabel(t),d=pi/100;,2020/7/16,第8页,d=pi/100; t=0:d:2*pi; x=sin(t); dxdt_diff=diff(x)/d; dxdt_grad=gradient(x)/d;,【例4.1-3】已知 采用diff和gradient计算该函数在区间 中的近似导函数。,subplot(1,2,1); plot(t,x,b);hold on plot(t,dxdt_grad,m,LineWidth,8) plot(t(1:end-1),dxdt_diff,.k,MarkerSize,8) axis(0,2*pi,-1.1,1.1); title
6、(0, 2pi) legend(x(t),dxdt_grad,dxdt_diff,Location,North) xlabel(t), box off;hold off,subplot(1,2,2) kk=(length(t)-10):length(t); hold on; plot(t(kk),dxdt_grad(kk),om,MarkerSize,8) plot(t(kk-1),dxdt_diff(kk-1),.k,MarkerSize,8) title(end-10, end) legend(dxdt_grad,dxdt_diff,Location,SouthEast) xlabel(t
7、),box off; hold off,宏观上, diff和gradient结果大致相同,细节上, diff和gradient数值有差异， diff没有给出最后一点导数,2020/7/16,第9页,Sx=sum(X) % sum按列向求和得(1n)数组 Scs=cumsum(X) %沿X列向求累计和, 仍是数组, 第(i, k)个元素是X数组第k列前i个元素的和。最后一行等于Sx,4.1.2 数值求和与近似数值积分,St=trapz(t,X) 或 St=dt*trapz(X) %梯形法求積分 Sct=cumtrapz(t,X) 或 Sct=dt*cumtrapz(X) %梯形法沿列向求X关于x
8、的累计积分,最后一个值等于St,S=dt*sum(X), S=sum(t,X) %近似矩形法求积分,Scs=cumsum(t,X) =dt*cumsum(X),2020/7/16,第10页,clear; d=pi/8; t=0:d:pi/2; y=0.2+sin(t); s=sum(y); s_sa=d*s;% s_sa=sum(t, y), 近似矩形法积分 s_ta=trapz(t,y); %梯形法积分,disp(sum求得积分,blanks(3),trapz求得积分) disp(s_sa, s_ta) t2=t,t(end)+d; y2=y,nan; stairs(t2,y2,:k);ho
9、ld on,plot(t,y,r,LineWidth,3) h=stem(t,y,LineWidth,2); set(h(1),MarkerSize,10) axis(0,pi/2+d,0,1.5) hold off; shg sum求得积分 trapz求得积分 1.5762 1.3013,2020/7/16,第11页,4.1.3 计算精度可控的数值积分,一重积分(quadrature精度可控)： S1=quad(fun,a,b,tol) %自适应Simpson法 S2=quadl(fun,a,b,tol) %自适应罗巴托 Lobatlo法 二重积分(精度可控) ： S3=dblquad(fu
10、n,xmin,xmax,ymin,ymax,tol) 三重积分(精度可控) ： S4=triplequad(fun,xmin,xmax,ymin,ymax, zmix,zmax,tol) fun:可为字符串,内联对象,匿名函数,函数句柄,注意乘除法运算一定加点.用数组运算 tol: 默认积分的绝对精度为10-6,2020/7/16,第12页,（1） syms x Isym=vpa(int(exp(-x2),x,0,1) % x.2数组乘方亦可 Isym = 0.74682413281242702539946743613185,例4.1-5：求 .,% （2） 梯形法积分 format long
11、 d=0.001;x=0:d:1; Itrapz=d*trapz(exp(-x.*x) % x.*x必须为数组乘 Itrapz = 0.746824071499185,（3） fx=exp(-x.2); %一定用数组乘 Ic=quad(fx,0,1,1e-8) Ic = 0.746824132854452 %实际精度控制到1e-10,2020/7/16,第13页,（1）符号计算法 syms x y s=vpa(int(int(xy,x,0,1),y,1,2) Warning: Explicit integral could not be found. s = 0.405465108108164
12、38197801311546435,例4.1-6：求 .,（2）数值积分法 format long s_n=dblquad(x,y)x.y,0,1,1,2) %匿名函数表示被积函数 s_n = 0.405466267243508,s_n=dblquad(x.y,0,1,1,2) %字符串表示被积函数 s_n=dblquad(inline(x.y),0,1,1,2) %内联函数表示被积函数,2020/7/16,第14页,4.1.4 函数极值的数值求解,x,fval,exitflag,output=fminbnd(fun,x1,x2,options) %一元函数在区间bound(x1,x2)中极小
13、值 x,fval: 极值点坐标和对应目标函数极值 x,fval,exitflag,output=fminsearch(fun,x0,options) %单纯形法求搜索起点x0附近多元函数极值点 % x每列代表一个候选极值点，各列按目标函数极小值递增顺序 %x(:,1)对应的目标函数极小值点由fval给出 fun: 字符串,内联函数,匿名函数,函数句柄,注意乘除法运算一定加点.用数组运算 options: 配置优化参数,可略 exitflag: 给出大于0的数,则成功搜索到极值点 output: 给出具体的优化算法和迭代次数,2020/7/16,第15页,【例4.1-7 】已知 ，在-10 x1
14、0区间，求函数的最小值。,（1）用“导数为零法”求极值点 syms x y=sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1); yd=diff(y,x); %求导函数 xs0=solve(yd,x) %求导函数为零的根，即极值点 yd_xs0=vpa(subs(yd,x,xs0) %验证极值点处导函数为零 y_xs0=vpa(subs(y,x,xs0) %求极值点处极值 xs0 = matrix(0.050838341410271656880659496266968) yd_xs0 = 2.2958874039497802890014385492622*10(-41
15、) y_xs0 = -0.001263317776974196724544154344118 无法判断是否最小值,2020/7/16,第16页,【例4.1-7 】已知 ，在-10 x10区间，求函数的最小值。,（2）采用优化算法求极小值 x1=-10;x2=10; yx=(x)(sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1); xn0,fval,exitflag,output=fminbnd(yx,x1,x2) xn0 = 2.514797840754235 fval = %比“导数为零法”求得的极值更小, 更可能是最小值 -0.499312445280039 e
16、xitflag = 1 output = iterations: 13 funcCount: 14 algorithm: golden section search, parabolic interpolation message: 1x112 char,2020/7/16,第17页,【例4.1-7 】已知 ，在-10 x10区间，求函数的最小值。,（4）据图形观察，重设fminbnd的搜索区间 x11=6;x2=10; yx=(x)(sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1); xn00,fval,exitflag,output=fminbnd(yx,x11,x2) xn00 = 8.023562824723015 fval = -3.568014059128578 %最
