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

1、2022/8/19第1页第4章 MATLAB 的数值计算 4.1 数值微积分4.2 矩阵和代数方程4.3 概率分布和统计分析4.4 多项式运算和卷积2022/8/19第2页4.1数值微积分4.1.1 近似数值极限及导数4.1.2 数值求和与近似数值积分4.1.3 计算精度可控的数值积分4.1.4 函数极值的数值求解4.1.5 常微分方程的数值解2022/8/19第3页在MATLAB数值计算中，既没有专门的求极限指令，也没有专门的求导指令。但MATLAB提供了与“求导”概念有关的“求差分”指令。dx=diff(X) %计算向量X的前向差分FX=gradient(F) %求一元(函数)梯度FX,

2、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 近似数值极限及导数2022/8/19第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:end)-F(1:end-2)/2 , FX(end

3、) = F(end)-F(end-1) ; FX长度与F的长度相同当F是矩阵时， FX, FY是与F同样大小的矩阵。 FX的每行对应F相应行元素间的梯度 ; FY的每列对应F相应列元素间的梯度 ; 4.1.1 近似数值极限及导数2022/8/19第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;

4、 %可得到正确结果2022/8/19第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)的数值十分接近，高位有效数字完全相同，

5、df =f(t+dt)-f(t) 造成高位有效数字消失，精度急剧变差。2022/8/19第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/8/19第8页d=pi/100; t=0:d:2*pi;x=sin(t);dxdt_diff=diff(x)/d;dxdt_grad=gradien

6、t(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):

7、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/8/19第9页2022/8/19第10页Sx=sum(X) % su

8、m按列向求和得(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/8/19第11页clear; d=pi/8; t=0:d:pi/2;

9、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

10、求得积分 trapz求得积分 1.5762 1.30132022/8/19第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:可为字符串,内联对象,匿名函数,函数句柄,注意乘除法运算

11、一定加点.用数组运算 tol: 默认积分的绝对精度为10-62022/8/19第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.

12、746824132854452 %实际精度控制到1e-102022/8/19第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) %字符串表示被

13、积函数s_n=dblquad(inline(x.y),0,1,1,2) %内联函数表示被积函数% 一定为数组乘2022/8/19第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)对应的目标函数极小值

14、点由fval给出 fun: 字符串,内联函数,匿名函数,函数句柄,注意乘除法运算一定加点.用数组运算 options: 配置优化参数,可略 exitflag: 给出大于0的数,则成功搜索到极值点 output: 给出具体的优化算法和迭代次数2022/8/19第16页【例4.1-7 】已知 ，在-10 x10区间，求函数的最小值。（1）用“导数为零法”求极值点syms xy=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

15、) %验证极值点处导函数为零y_xs0=vpa(subs(y,x,xs0) %求极值点处极值xs0 = matrix(0.050838341410271656880659496266968)yd_xs0 = 2.2958874039497802890014385492622*10(-41)y_xs0 = -0.001263317776974196724544154344118 无法判断是否最小值2022/8/19第17页【例4.1-7 】已知 ，在-10 x10区间，求函数的最小值。（2）采用优化算法求极小值x1=-10;x2=10;yx=(x)(sin(x)2*exp(-0.1*x)-0.5

16、*sin(x)*(x+0.1); xn0,fval,exitflag,output=fminbnd(yx,x1,x2)xn0 = 2.514797840754235fval = %比“导数为零法”求得的极值更小, 更可能是最小值 -0.499312445280039exitflag = 1output = iterations: 13 funcCount: 14 algorithm: golden section search, parabolic interpolation message: 1x112 char 2022/8/19第18页【例4.1-7 】已知 ，在-10 x10区间，求函

17、数的最小值。（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.023562824723015fval = -3.568014059128578 %最小值exitflag = 1output = iterations: 9 funcCount: 10 algorithm: golden section search, parabolic interpolation message

18、: 1x112 char y=sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);（3）绘图观察最小值xx=-10:pi/200:10;yxx=subs(y,x,xx);plot(xx,yxx)xlabel(x),grid on 2022/8/19第19页2022/8/19第20页例4.1-8: f(x,y)=100(y-x2)2+(1-x)2在区间-5,5的极小值ff=(x)100*(x(2)-x(1).2).2+(1-x(1).2; %匿名函数 x0=-5,-2,2,5;-5,-2,2,5; %4个搜索起点 sx,sfval,sexit,soutput=fmins

19、earch(ff,x0)其理论极小值点为x=1, y=1%sx给出一组使优化函数值非减的局部极小点 sx = 0.99998 -0.68971 0.41507 8.0886 0.99997 -1.9168 4.9643 7.8004sfval = 2.4112e-0102.4112e-010 5.7525e+002 2.2967e+003 3.3211e+005 format short e %取5位科学计数法 disp(ff(sx(:,1),ff(sx(:,2),ff(sx(:,3),ff(sx(:,4) 用x的二元向量表示x,yx0 =-5 -2 2 5 -5 -2 2 52022/8/1

20、9第21页4.1.5 常微分方程Ordinary Differential Equation的数值解t,Y=ode45(odefun,tspan,y0) % 4阶龙格库塔 数值法odefun: 待求解一阶微分方程组的函数文件句柄tspan: 自变量微分二元区间t0, tfy0: 一阶微分方程组的(n*1)初值列向量matlab为解常微分方程初值问题提供了一组配套齐全,结构严整的指令, 包括: ode45, ode23, ode113, ode23t, ode15s, ode23s, ode23tb.在此只介绍最常用的ode45的基本使用方法.ode45使用方法:t二元区间的点系列Y原函数在微分

21、区间点系列上的函数值2022/8/19第22页例4.1-9求解：%解算微分方程tspan=0,30; y0=1;0;tt,yy=ode45(DyDt,tspan,y0); figure(1)plot(tt,yy(:,1)xlabel(t),title(x(t)解:令，上式写成一阶微分方程组形式%画相平面图(函数和其导数勾画的曲线称为相轨迹)figure(2)plot(yy(:,1),yy(:,2)xlabel(位移),ylabel(速度) function ydot=DyDt(t,y)mu=2; ydot=y(2);mu*(1-y(1)2)*y(2)-y(1);据以上方程组，编写M函数文件Dy

22、Dt.m2022/8/19第23页2022/8/19第24页4.2矩阵和代数方程4.2.1 矩阵运算和特征参数4.2.2 矩阵的变换和特征值分解4.2.3 线性方程的解4.2.4 一般代数方程的解2022/8/19第25页4.2.1 矩阵运算和特征参数矩阵与标量之间的四则运算与数组运算相同矩阵和矩阵之间的四则运算矩阵和矩阵之间的加减运算与数组运算相同设 A 是一个 mn 矩阵，B 是一个 pq 矩阵，当 np 时，两个矩阵可以相乘，乘积为 mq 矩阵。矩阵乘法不可逆。在 MATLAB 中，矩阵乘法由“*”实现。矩阵除法在实际中主要用于求解线性方程组矩阵转置：符号“ ”实现矩阵的转置操作。对于实

23、数矩阵， “ ”表示矩阵转置，对于复数矩阵，“ ”实现共轭转置。对于复数矩阵，如果想要实现非共轭转置，可以使用符号“ . ”。2022/8/19第26页format rat %有理格式显示A=magic(2) + j*pascal(2) A = 1 + 1i 3 + 1i 4 + 1i 2 + 2iA1=A A1=1 - 1i 4 - 1i %共轭转置 3 - 1i 2 - 2iA2=A. A2=1 + 1i 4 + 1i %非共轭转置，数组运算操作 3 + 1i 2 + 2i例4.2-2 矩阵和数组转置操作的差别2022/8/19第27页计算矩阵标量特征参数秩,迹,行列式的指令rank(A)

24、 %求秩（Rank）det(A) %求行列式（Determinant）trace(A) %求迹（Trace），即矩阵主对角元素的和2022/8/19第28页 A=reshape(1:9,3,3); r=rank(A) %求秩 d3=det(A) %非满秩矩阵的行列式一定为零 d2=det(A(1:2,1:2) %求子式的行列式 t=trace(A)【例4.2-3】矩阵标量特征参数计算示例。A = 1 4 7 2 5 8 3 6 9r = 2d3 =0d2 = -3t = 15 2022/8/19第29页【例4.2-4】矩阵标量特征参数的性质。format short; rand(twister

25、,0)A=rand(3,3); B=rand(3,3);C=rand(3,4); D=rand(4,3);tAB=trace(A*B) %任何符合矩阵乘法规则的两个矩阵乘积tBA=trace(B*A) %的“迹”不变。同阶乘积“迹”不变tCD=trace(C*D) %两个“內维”相等矩阵的乘积“迹”不变tDC=trace(D*C)tAB = 2.6030tBA = 2.6030tCD = 4.1191tDC = 4.1191 dCD=det(C*D)dDC=det(D*C) dCD = 0.0424 dDC = -2.6800e-018d_A_B=det(A)*det(B)dAB=det(A*

26、B) dBA=det(B*A) 同阶矩阵乘积行列式等于各矩阵行列式之乘积d_A_B = 0.0094dAB = 0.0094dBA = 0.0094 非同阶矩阵乘积行列式不等于各矩阵行列式之乘积2022/8/19第30页dx=diff(X) %计算向量X的前向差分FX=gradient(F) %求一元(函数)梯度FX, FY =gradient(F) %求二元（函数）梯度4.1.1 近似数值极限及导数S1=sum(X,1) % sum按列向求和得(1n)数组, =sum(X) (X多行)S2=sum(X,2) % sum按行向求和得(n1)数组Scs=cumsum(X)%沿X列向求累计和, 仍

27、是数组, 最后一行等于Sx4.1.2 数值求和与近似数值积分St=trapz(t,X) 或 St=dt*trapz(X) %梯形法求積分Sct=cumtrapz(t,X) 或 Sct=dt*cumtrapz(X) %梯形法沿列向求X关于x的累计积分,最后一个值等于StX= 1 2; 3 4 S1= 4 6 S2= 3 7sum(S1) = 10sum(S1,1) = 4 6Review2022/8/19第31页4.1.3 计算精度可控的数值积分一重积分(quadrature精度可控)：S1=quad(fun,a,b,tol) %自适应Simpson法S2=quadl(fun,a,b,tol)

28、%自适应罗巴托 Lobatlo法二重积分(精度可控) ：S3=dblquad(fun,xmin,xmax,ymin,ymax,tol)三重积分(精度可控) ： S4=triplequad(fun,xmin,xmax,ymin,ymax, zmix,zmax,tol) fun:可为字符串,内联对象,匿名函数,函数句柄,注意乘除法运算一定加点.用数组运算 tol: 默认积分的绝对精度为10-6Review2022/8/19第32页4.1.4 函数极值的数值求解 x,fval,exitflag,output=fminbnd(fun,x1,x2,options) %一元函数在区间bound(x1,x2

29、)中极小值 x,fval: 极值点坐标和对应目标函数极值x,fval,exitflag,output=fminsearch(fun,x0,options) %单纯形法求搜索起点x0附近多元函数极值点 % x每列代表一个候选极值点，各列按目标函数极小值递增顺序 %x(:,1)对应的目标函数极小值点由fval给出 fun: 字符串,内联函数,匿名函数,函数句柄, 注意乘除法运算一定加点.用数组运算 options: 配置优化参数,可略 exitflag: 给出大于0的数,则成功搜索到极值点 output: 给出具体的优化算法和迭代次数Review2022/8/19第33页rank(A) %求秩（R

30、ank）det(A) %求行列式（Determinant）trace(A) %求迹(Trace)，即矩阵主对角元素的和Reviewt,Y=ode45(odefun,tspan,y0) % 4阶龙格库塔 数值法odefun: 待求解一阶微分方程组的函数文件句柄tspan: 自变量微分二元区间t0, tfy0: 一阶微分方程组的(n*1)初值列向量ode45使用方法:t二元区间的点系列Y原函数在微分区间点系列上的函数值4.1.5 常微分方程Ordinary Differential Equation的数值解4.2.1 矩阵运算和特征参数2022/8/19第34页4.2.2 矩阵的变换和特征值分解R

31、, ci=rref(A) %借助初等变换将A缩减行变成行阶梯矩阵R。B=orth(A)可得到矩阵A的正交基，B的列与A的列可张成相同的空间，而且B的列是正交的，因此B*B=eye(rank(A)，B的列数正好是A的秩。 X=null(A) %A矩阵零空间的全部正交基，满足AX=0, X*X=I。B=orth(A)V, D=eig(A) %A矩阵的特征值、特征向量分解，使AV=VD。ci 是行数组, 其元素表示A中线性独立列的序号。length(ci)=rank(A)2022/8/19第35页【例4.2-5】行阶梯阵简化指令rref计算结果的含义。A=magic(4)R,ci=rref(A) %

32、行阶梯分解A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1R = 1 0 0 1 0 1 0 3 0 0 1 -3 0 0 0 0ci = 1 2 3 r_A=length(ci) % =rank(A)r_A = 3 aa=A(:,1:3)*R(1:3,4) %A前三列线形组合err=norm(A(:,4)-aa) aa = 13 8 12 1err = 0 2022/8/19第36页A=reshape(1:15,5,3) X=null(A)S=A*XS =1.0e-014 * 0 -0.1776 -0.2665 -0.3553 -0.5329A =1 6

33、11 2 7 12 3 8 13 4 9 14 5 10 15n= 3l= 1Rank_A=2X =0.4082 -0.8165 0.4082ans = 1 例4.2-6 矩阵零空间及含义。设 是矩阵 的零空间,即n=size(A,2)l=size(X,2)Rank_A=rank(A)n-l=rank(A)2022/8/19第37页【例4.2-7】简单实阵的特征分解。eig，cdf2rdf Complex Diagonal Form to Real-block Diagonal Form%（1）A=1,-3;2,2/3V,D=eig(A) A =1.0000 -3.0000 2.0000 0.

34、6667V = 0.7746 0.7746 0.0430 - 0.6310i 0.0430 + 0.6310iD =0.8333 + 2.4438i 0 0 0.8333 - 2.4438i %（2）VR,DR=cdf2rdf(V,D) VR =0.7746 0 0.0430 -0.6310DR =0.8333 2.4438 -2.4438 0.8333 %（3）A1=V*D/V %因计算误差有很小虚部A1_1=real(A1) %去除虚部后等于AA2=VR*DR/VRerr1=norm(A-A1,fro)err2=norm(A-A2,fro) A1 =1.0000 + 0.0000i -3.

35、0000 2.0000 - 0.0000i 0.6667 A1_1 =1.0000 -3.0000 2.0000 0.6667A2 =1.0000 -3.0000 2.0000 0.6667err1 =6.7532e-016err2 =4.4409e-016 Frobenius矩阵范数 实数块对角阵实阵矩阵Frobenius范数,类似向量2范数，不同于矩阵2范数2022/8/19第38页4.2.2 线性方程的解对于方程组Amnx=b (m个方程，n个未知数)， 当向量b在矩阵A列向量所张空间中, rank(A,b) =rank(A)= r a) 若n =r，则解唯一。 b) 若n r，则解不唯

36、一。 当向量b不在矩阵A列向量所张空间中, 则无准确解但存在最小二乘解。1. 线性方程解的一般讨论matlab定义的左除运算可以很方便地解上述方程组： x=Ab (x=A-1b只适用于A非奇异时)2.除法运算解方程解 当m=n时, “恰定”方程 当mn时, “超定”方程 当mr, 解不唯一% （3）求特解和通解，并对由之构成的全解进行验算xs=Ab; xg=null(A); % xg是齐次方程Ax=0的解c=rand(1); ba=A*(xs+c*xg) % ba是A乘 “一个随机的全解”Warning: Rank deficient, rank = 2, tol = 1.8757e-014.

37、ba =13.0000 14.0000 15.0000 16.0000 例4.2-8: 求方程 的解norm(ba-b) ans = 1.1374e-014 2022/8/19第40页例4.2-9: “逆阵法”和”左除法”解恰定方程的性能对比。测试阵产生指定异常值和特殊带宽的随机阵3. 矩阵的逆A_1=inv(A) %求非奇异阵A的逆randn(state,0); A=gallery(randsvd,300,2e13); %产生条件数 为2e13的300阶随机矩阵;x=ones(300,1); %定义真解b=A*x;cond(A) %验算矩阵条件数, 结果a1.9978e+013, 值越大,

38、阵越病态,%求逆法tic; xi=inv(A)*b; ti=toceri=norm(x-xi) rei=norm(A*xi-b)/norm(b) %左除法tic; xd=Ab;td=tocerd=norm(x-xd)red=norm(A*xd-b)/norm(b)ti =0.0185eri =0.0883rei =0.0051td =0.0066erd =0.0298red =8.7810e-015矩阵计算对于误差越敏感, 数值稳定性差 2022/8/19第41页求解任意函数f(x)=0(可能无解,单解,多解)的步骤:作图获取初步近似解 观察f(x)与横轴的交点坐标，用zoom放大，用ginp

39、ut得较精确些的交点坐标值。4.3 一般代数方程的解 用”泛函”指令求精确解 x,favl = fzero(fun,x0) %一元函数求零点 x,fval = fsolve(fun,x0) %解非线性方程组 fun: 内联对象,匿名函数,函数句柄,字符串;被解函数自变量一般用x, 注意乘除法运算一定加点.用数组运算 x0: 零点初始猜测值. x0为标量时取与之最靠近的零点; x0取a,b时,在此区间内寻找一个零点. x: 所求零点的自变量值 fval: 函数值2022/8/19第42页例4.2-10: 求f(t)=(sin2t)e-0.1t0.5|t|的零点y_C=inline(sin(t).

40、2.*exp(-0.1*t)-0.5*abs(t),t); t=-10:0.01:10;Y=y_C(t); clf, plot(t,Y,r); hold onplot(t,zeros(size(t),k);xlabel(t);ylabel(y(t) t4 = 0.5993y4 = 0tt =-2.0039 -0.5184 -0.0042 0.6052 1.6717 t3 = -0.5198y3 = 5.5511e-017zoom ontt,yy=ginput(5),zoom offt4,y4=fzero(y_C,tt(4) t3,y3=fzero(y_C,tt(3) %t=0处没穿越横轴202

41、2/8/19第43页4.3 概率分布和统计分析4.3.1 概率函数、分布函数、逆分布函数和随机数的发生4.3.2 随机数发生器和 统计分析指令2022/8/19第44页4.3 概率分布和统计分析1. 二项分布（Binomial distribution ）4.3.1 概率函数、分布函数、逆分布函数和随机数的发生pk=binopdf(k,N,p) 事件A发生k次的概率Fk=binocdf(k,N,p) 事件A发生次数不大于k的概率R=binornd(N,p) 产生符合二项分布B(N,p)的(mn) 随机数组（元素值为事件A可能发生的次数）。Binomial Cumulative Distribu

42、tion FunctionBinomial Probability Density Function2022/8/19第45页N=100;p=0.5;k=0:N;pdf=binopdf(k,N,p); cdf=binocdf(k,N,p);h=plotyy(k,pdf,k,cdf);set(get(h(1),Children),Color,b,Marker,.,MarkerSize,13)set(get(h(1),Ylabel),String,pdf)set(h(2),Ycolor,1,0,0)set(get(h(2),Children),Color,r,Marker,+,MarkerSize,4)set(get(h(2),Ylabel),String,cdf)xlabel(k) grid on 例4.3-1. 画出N=100，p=0.5情况下的二项分布概率特性曲线。4.3.1 概率函数、分布函数、逆分布函数和随机数的发生2022/8/19第46页px=normpdf(x,M