山东大学MATLAB编程指导课件第4章 MATLAB 的数值计算

2022/4/20第1页第第4章章 MATLAB 的数值计算的数值计算 4.1 数值微积分数值微积分4.2 矩阵和代数方程矩阵和代数方程4.3 概率分布和统计分析概率分布和统计分析4.4 多项式运算和卷积多项式运算和卷积2022/4/20第2页4.1数值微积分数值微积分4.1.1 近似数值极限及导数近似数值极限及导数4.1.2 数值求和与近似数值积分数值求和与近似数值积分4.1.3 计算精度可控的数值积分计算精度可控的数值积分4.1.4 函数极值的数值求解函数极值的数值求解4.1.5 常微分方程的数值解常微分方程的数值解2022/4/20第3页u在在MATLAB数值计算中，既没有专门的求极限指令，数值计算中，既没有专门的求极限指令，也没有专门的求导指令。但也没有专门的求导指令。但MATLAB提供了与提供了与“求导求导”概念有关的概念有关的“求差分求差分”指令。指令。udx=diff(X) %计算向量计算向量X的前向差分的前向差分uFX=gradient(F) %求一元求一元(函数函数)梯度梯度uFX, FY =gradient(F) %求二元（函数）梯度求二元（函数）梯度u对对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 近似数值极限及导数近似数值极限及导数2022/4/20第4页udx=diff(X) %计算向量计算向量X的前向差分的前向差分uFX=gradient(F) %求一元求一元(函数函数)梯度梯度uFX, FY =gradient(F) %求二元（函数）梯度求二元（函数）梯度u对对gradient而言，当而言，当F是向量时，是向量时，FX(1) = F(2)-F(1), FX(2:end-1) = (F(3:end)-F(1:end-2)/2 , FX(end) = F(end)-F(end-1) ; FX长度与长度与F的长度相同的长度相同u当当F是矩阵时，是矩阵时， FX, FY是与是与F同样大小的矩阵。同样大小的矩阵。 FX的的每行对应每行对应F相应行元素间的梯度相应行元素间的梯度 ; FY的每列对应的每列对应F相应相应列元素间的梯度列元素间的梯度 ; 4.1.1 近似数值极限及导数近似数值极限及导数2022/4/20第5页数值极限和导数的应用应十分谨慎数值极限和导数的应用应十分谨慎x=eps;L1=(1-cos(2*x)/(x*sin(x), L2=sin(x)/x,L1 = 0L2 = 1syms tf1=(1-cos(2*t)/(t*sin(t);f2=sin(t)/t;Ls1=limit(f1,t,0)Ls2=limit(f2,t,0) Ls1 = 2Ls2 = 1 x=pi/1000; %可得到正确结果可得到正确结果2022/4/20第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)hold onplot(t,dxdt_eps)hold offlegend(x(t),dx/dt)xlabel(t) 【例例4.1-2】已知已知, 求该函数在区间求该函数在区间 中的近似导函数。中的近似导函数。数值导数受计算中有限精度困扰，当增量数值导数受计算中有限精度困扰，当增量dt过小时，过小时，f(t+dt)与与f(t)的数值十分接近，高位有效数字完全相的数值十分接近，高位有效数字完全相同，同， df =f(t+dt)-f(t) 造成高位有效数字消失，精度造成高位有效数字消失，精度急剧变差。急剧变差。2022/4/20第7页数值极限和导数的应用应十分谨慎数值极限和导数的应用应十分谨慎%（2）自变量的增量取得适当）自变量的增量取得适当x_d=sin(t+d);dxdt_d=(x_d-x)/d;plot(t,x,LineWidth,5)hold onplot(t,dxdt_d)hold offlegend(x(t),dx/dt)xlabel(t) 【例例4.1-2】已知已知, 求该函数在区间求该函数在区间 中的近似导函数。中的近似导函数。d=pi/100;2022/4/20第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 onplot(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(0, 2pi)legend(x(t),dxdt_grad,dxdt_diff,Location,North)xlabel(t), box off;hold offsubplot(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),box off; hold off 宏观上宏观上, diff和和gradient结果大致相同结果大致相同细节上细节上, diff和和gradient数值有差异，数值有差异， diff没有给出最后一点导数没有给出最后一点导数2022/4/20第9页2022/4/20第10页Sx=sum(X) % sum按列向求和得按列向求和得(1n)数组数组Scs=cumsum(X)%沿沿X列向求列向求累计和累计和, 仍是数组仍是数组, 第第(i, k)个元素是个元素是X数组第数组第k列前列前i个个元素的和。最后一行等于元素的和。最后一行等于Sx4.1.2 数值求和与近似数值积分数值求和与近似数值积分St=trapz(t,X) 或或 St=dt*trapz(X) %梯形法求積分梯形法求積分Sct=cumtrapz(t,X) 或或 Sct=dt*cumtrapz(X) %梯形法沿列向求梯形法沿列向求X关于关于x的的累计积分累计积分,最后一个值等于最后一个值等于StS=dt*sum(X), S=sum(t,X) %近似矩形法求积分近似矩形法求积分Scs=cumsum(t,X) =dt*cumsum(X) 2022/4/20第11页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); %梯形法积分梯形法积分【例例4.1-4】求积分求积分 , 其中其中disp(sum求得积分求得积分,blanks(3),trapz求得积分求得积分)disp(s_sa, s_ta)t2=t,t(end)+d; y2=y,nan;stairs(t2,y2,:k);hold onplot(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.30132022/4/20第12页4.1.3 计算精度可控的数值积分计算精度可控的数值积分一重积分一重积分(quadrature精度可控精度可控)：S1=quad(fun,a,b,tol) %自适应自适应Simpson法法S2=quadl(fun,a,b,tol) %自适应罗巴托自适应罗巴托 Lobatlo法法二重积分二重积分(精度可控精度可控) ：S3=dblquad(fun,xmin,xmax,ymin,ymax,tol)三重积分三重积分(精度可控精度可控) ： S4=triplequad(fun,xmin,xmax,ymin,ymax, zmix,zmax,tol) fun:可为字符串可为字符串,内联对象内联对象,匿名函数匿名函数,函数句柄函数句柄,注意乘除法注意乘除法运算一定加点运算一定加点. .用数组运算用数组运算 tol: 默认积分的绝对精度为默认积分的绝对精度为10-62022/4/20第13页（1）syms xIsym=vpa(int(exp(-x2),x,0,1) % x.2数组乘方亦可数组乘方亦可Isym = 0.74682413281242702539946743613185 例例4.1-5：求：求 . % （2） 梯形法积分梯形法积分format longd=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-102022/4/20第14页（1）符号计算法）符号计算法syms x ys=vpa(int(int(xy,x,0,1),y,1,2)Warning: Explicit integral could not be found. s = 0.40546510810816438197801311546435 例例4.1-6：求：求 . （2）数值积分法）数值积分法format longs_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) %内联函数表示被积函数内联函数表示被积函数% 一定为数组乘一定为数组乘2022/4/20第15页4.1.4 函数极值的数值求解函数极值的数值求解 x,fval,exitflag,output=fminbnd(fun,x1,x2,options) %一元函数一元函数在区间在区间bound(x1,x2)中极小值中极小值 x,fval: 极值点坐标和对应目标函数极值极值点坐标和对应目标函数极值x,fval,exitflag,output=fminsearch(fun,x0,options) %单纯形法求搜索起点单纯形法求搜索起点x0附近附近多元函数多元函数极值点极值点 % x每列代表一个候选极值点，各列按每列代表一个候选极值点，各列按目标函数极小值递增顺序目标函数极小值递增顺序 %x(:,1)对应的目标函数极小值点由对应的目标函数极小值点由fval给出给出 fun: 字符串字符串,内联函数内联函数,匿名函数匿名函数,函数句柄函数句柄,注意乘除法运注意乘除法运算一定加点算一定加点. .用数组运算用数组运算 options: 配置优化参数配置优化参数,可略可略 exitflag: 给出大于给出大于0的数的数,则成功搜索到极值点则成功搜索到极值点 output: 给出具体的优化算法和迭代次数给出具