第三讲-1 MATLAB的数值计算.ppt

1、Matlab程序设计,第四章 数值计算,1 数值微分 2 数值积分 3 函数零、极值求解 4 微分方程求解 5 矩阵运算 6 线性方程求解 7 多项式运算,Matlab程序设计,1 数值微分,任意函数f(x)在x点的导数：,当h0时,引入差分 ：,则用f(x)在点x处的某种差商作为其导数。,MATLAB中用函数diff(x)表示差分：,FX=gradient(F) 求一元梯度 F为向量 FX,FY=gradient(F) 求二元梯度 F为矩阵,Matlab程序设计,1 数值微分, x=0:4;y=x.2 y = 0 1 4 9 16 dy=diff(y) dy = 1 3 5 7 gy=gra

2、dient(y) gy = 1 2 4 6 7,Matlab程序设计,1 数值微分, yy=rand(3)*10 yy = 4.2176 9.5949 8.4913 9.1574 6.5574 9.3399 7.9221 0.3571 6.7874 yy=fix(yy) yy = 4 9 8 9 6 9 7 0 6, fx,fy=gradient(yy) fx = 5.0000 2.0000 -1.0000 -3.0000 0 3.0000 -7.0000 -0.5000 6.0000 fy = 5.0000 -3.0000 1.0000 1.5000 -4.5000 -1.0000 -2.0

3、000 -6.0000 -3.0000,Matlab程序设计,1 数值微分,例：x=sin(t)，采用diff和gradient计算该函数在区间0，2*pi的 近似导数 d=pi/100; t=0:d:2*pi; x=sin(t); dxdt_diff=diff(x)/d; dxdt_grad=gradient(x)/d; 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(0, 2pi) legend

4、(x(t),dxdt_grad,dxdt_diff,Location,North) xlabel(t) hold off,Matlab程序设计,2 数值积分（数值求和、近似数值积分）,Sx=sum(X) 沿列求和,Scs=cumsum(X) 沿列累计求和,St=trapz(x,y) 采用梯形法沿列求函数y 关于自变量x的积分,St=cumtrapz(x,y)采用梯形法沿列求函数y 关于自变量x的累计积分,Matlab程序设计, x=fix(rand(2)*5) x = 1 1 4 4 sum(x) ans = 5 5 x=1:3 x = 1 2 3 sum(x) ans = 6, x=fix(

5、rand(4,2)*5) x = 1 2 0 1 1 4 3 2 cumsum(x) ans = 1 2 1 3 2 7 5 9,2 数值积分（数值求和、近似数值积分）,Matlab程序设计, x=1:4,y=x.2 x = 1 2 3 4 y = 1 4 9 16 trapz(x,y) ans = 21.5000, x=1:4,y=x.2 x = 1 2 3 4 y = 1 4 9 16 cumtrapz(x,y) ans = 0 2.5000 9.0000 21.5000,2 数值积分（数值求和、近似数值积分）,Matlab程序设计,例：x=sin(t)，采用sum和trapz计算该函数在

6、区间0，pi/2的数值积分 d=pi/10; t=0:d:pi/2; x=sin(t); sx=sum(x)*d %矩形法 st=trapz(t,x) %梯形法,2 数值积分（数值求和、近似数值积分）,sx = 1.1488 st = 0.9918,Matlab程序设计,2 数值积分（精度可控）,一重积分： q1=quad(fun,a,b) %自适应Simpson法 q2=quadl(fun,a,b) %自适应Lobatlo法 二重积分： SS=dblquad(fun,xmin,xmax,ymin,ymax,tol) 三重积分： SSS= triplequad(fun,xmin,xmax,ym

7、in,ymax, zmin,zmax,tol),Matlab程序设计,例：求 ,其精确值为0.74684 %fun.m function y=fun(x) y=exp(-x.*x); %主程序 Hf=fun; Isin1=quad(Hf,0,1) simpson法 Isin2=quadl(Hf,0,1) Lobatlo法 dt=0.01; t=0:dt:1; Isin3=dt*sum(fun(t) 矩形法 Isin4=trapz(t, fun(t) % 梯形法,Isin1 Isin2 Isin3 Isin4 0.74682418, 0.74682413, 0.75365739, 0.74681

8、800,Matlab程序设计,凡以函数为输入宗量的指令，被称为泛函指令。,3 函数零、极值求解,Matlab程序设计,例1(L4.5-1): 求f(t)=(sin2t)e-at-b|t|的零点 P1=0.1;P2=0.5; y_C=sin(x).2.*exp(-P1*x)-P2*abs(x); x=-10:0.01:10; Y=eval(y_C); clf,plot(x,Y,r);hold on plot(x,zeros(size(x),k); xlabel(t);ylabel(y(t) % yc=inline(sin(x).2.*exp(-P1*x)-P2*abs(x) % plot(x,y

9、c(P1,P2,x),zoom on tt,yy=ginput(5);zoom off t4,y4,exitflag=fzero(y_C,tt(4),P1,P2),Matlab程序设计,例2：f(x)=x2+3x+2在-5 5区间的极小值,并绘图 y=inline(x.2+3.*x+2), x=-5:0.05:5; plot(x,y(x), hold on f=fminbnd(y,-5,5) plot(f,y(f),or) 例3：f(x)=100(x2-x12)2+(1-x1)2在区间0，2的极小值 function f=xun(x) f=100*(x(2)-x(1).2).2+(1-x(1)

10、.2; x=fminsearch(xun,0,2),Matlab程序设计,数值求导 dxdt_diff=diff(x)/d; dxdt_grad=gradient(x)/d; 近似数值积分 S=dt*sum(x) 矩形法 S=trapz(t, x) % 梯形法 精度可控数值积分 q1=quad(fun,a,b) %自适应Simpson法 q2=quadl(fun,a,b) %自适应Lobatlo法 SS=dblquad(fun,xmin,xmax,ymin,ymax,tol) 函数零、极值求解 z=fzero(fun,x0) z=fminbnd(fun,x1,x2) z=fminsearch(

11、fun,x0),Matlab程序设计,4 微分方程求解,微分方程求解的仿真算法有多种，常用的有Euler（欧拉法）、Runge Kutta(龙 格-库塔法)。 Euler法称单步法，用于一阶微分方程,Matlab程序设计,当给定仿真步长时： 所以 yn+1 = yn + hf (xn,yn) n=0,1,2 y(x0)=y0,Matlab程序设计,Runge Kutta法 龙格-库塔法：实际上取两点斜率 的平均值来计算的，其精度高于欧拉算法 。 龙格-库塔法：ode23 ode45,k1=hf(xn,yn) k2=hf(xn+h,yn+k),Matlab程序设计,例： 为方便令 分别对y1,y

12、2求一阶导数，整理后写成一阶微分方程组形式 y1=y2 , y1(0)=0.25 y2=y2(1-y12)-y1 y2(0)=0 命令格式: T,Y = ode23(odefun,tspan,y0) 1. 建立m文件 2. 解微分方程,Matlab程序设计,方法一：建立函数m文件 function ydot=DyDt(t,y) ydot=zeros(2,1); ydot(1)=y(2); ydot(2)=y(2)*(1-y(1)2)-y(1); 给定区间、初始值;求解微分方程 t0=0; tf=20; y0=0.25 0; t,y=ode23(DyDt, t0, tf, y0) plot(t,

13、y), figure(2),plot(y(:,1),y(:,2),Matlab程序设计,方法二：利用句柄函数求解： 建立m文件 function dxdt=DyDt(t,y) dxdt=y(2);y(2)*(1-y(1)2)-y(1); 求解微分方程 t,y=ode23(DyDt,0 30,0 0.25); plot(t,y); figure(2) plot(y(:,1),y(:,2),Matlab程序设计,Matlab程序设计,5 矩阵运算,矩阵加、减、乘、除 A+B；A-B；a+B；a-B；a*B;A*B; (A，B为矩阵，a为标量) 矩阵转置B=A A=magic(2)+j*pascal

14、(2) A = 1.0000 + 1.0000i 3.0000 + 1.0000i 4.0000 + 1.0000i 2.0000 + 2.0000i A1=A , A2=A. A1 = 1.0000 - 1.0000i 4.0000 - 1.0000i 3.0000 - 1.0000i 2.0000 - 2.0000i A2 = 1.0000 + 1.0000i 4.0000 + 1.0000i 3.0000 + 1.0000i 2.0000 + 2.0000i,Matlab程序设计,5 矩阵运算,矩阵求逆 X=inv(A) 求方阵的逆 X*A为单位矩阵 inv(magic(2)*magic

15、(2) ans = 1.0000 -0.0000 0 1.0000 X=pinv(A) 求矩阵的伪逆 A*X*A=A，A*X和X*A都是Hermitian矩阵, a=rand(2,3),a*pinv(a)*a a = 0.6787 0.7431 0.6555 0.7577 0.3922 0.1712 ans = 0.6787 0.7431 0.6555 0.7577 0.3922 0.1712,Matlab程序设计,5 矩阵运算,矩阵标量特征参数 rank(A)：秩 矩阵A中线性无关列（行）向量组中最大向量数 矩阵A中最高非零子行列式的阶数 矩阵A中最高非奇异子矩阵的维数 trace(A)：迹

16、 矩阵A主对角元素之和 det(A)：行列式,Matlab程序设计,5 矩阵运算,矩阵标量特征参数 A=fix(rand(3)*10) A = 0 3 4 8 9 3 6 0 7 r=rank(A),t=trace(A),d=det(A) r = 3 t = 16 d = -330,Matlab程序设计,5 矩阵运算,矩阵变换和特征值分解 R,ci=rref(A)：将矩阵A变换成行阶梯矩阵R ci是行数组。它的元素表示A中线性独立“列”的序号。 X=null(A)：矩阵A零空间的全部正交基，AX=0 Z=orth(A)：矩阵A值空间的全部正交基。 V,D=eig(A)：矩阵A的特征值、特征向量

17、分解。 AV=VD。该指令只用特征值不同的情况。如果特征值相同，采用指令jordan，进行Jordan分解。,Matlab程序设计,5 矩阵运算, A=magic(4),R,ci=rref(A) A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 R = 1 0 0 1 0 1 0 3 0 0 1 -3 0 0 0 0 ci = 1 2 3,Matlab程序设计,5 矩阵运算, A=reshape(1:15,5,3);X=null(A),S=A*X, rank(A),size(A,2)-size(X,2) X = 0.4082 -0.8165 0.4082 S

18、 = 1.0e-014 * -0.2665 -0.1776 -0.0888 -0.0888 -0.0888 ans = 2 ans = 2,Matlab程序设计,6 线性方程求解,对于方程ax=b，a 为anm矩阵(n方程个数，m为未知数个数)，有三种情况： 当n=m时，此方程成为“恰定”方程 当nm时，此方程成为“超定”方程 当nm时，此方程成为“欠定”方程 matlab定义的除运算可以很方便地解上述三种方程,Matlab程序设计,求解线性方程ax=b三法： (1) x=a-1 b 矩阵逆(奇异值分解法,解可靠), 解：x=inv(a)b 欠定和超定方程时用伪逆求解(最小范数的解)： x=p

19、inv(a)*b (2) x=ab 左除运算(解基本可靠,速度快)， 解：x=ab (3) 正则方程法(精度差)， 解: x=inv(a*a)*a*b,Matlab程序设计,解: a=1 2;2 3;b=8;13; x=inv(a)*b x=ab,a x = b,例: x1+2x2=8 2x1+3x2=13,x = x = 2.00 2.00 3.00 3.00,Matlab程序设计,例: x1+2x2=1 2x1+3x2=2 3x1+4x2=3,a x = b,a=1 2;2 3;3 4;b=1;2;3; 解1: x=ab 解2: x=inv(aa) a b,x = x = 1.00 1.0

20、0 0 0.00,Matlab程序设计,x1+2x2+3x3=1 2x1+3x2+4x3=2,a x = b,a=1 2 3;2 3 4;b=1;2; x=ab x=pinv(a)b,x = x = 1.00 0.83 0 0.33 0 -0.17,Matlab程序设计,7 多项式运算,matlab语言把多项式表达成一个行向量， 该向量中的元素是按多项式降幂排列的。 P(x)=anxn+an-1xn-1+a0 可用行向量 p=an an-1 a1 a0表示,Matlab程序设计,多项式运算函数(P160 表4.4-1),Matlab程序设计,7.1 poly 产生特征多项式系数向量 特征多项式

21、一定是n+1维的 特征多项式第一个元素一定是1 例:a=1 2 3;4 5 6;7 8 0; p=poly(a) 产生特征多项式系数向量 得：p =1.00 -6.00 -72.00 -27.00 p是多项式p(x)=x3-6x2-72x-27的matlab描述方法，我们可用： p1=poly2str(p,x) 显示数学多项式的形式 得：p1 =x3 - 6 x2 - 72 x - 27,Matlab程序设计,7.2 roots 求多项式的根,a=1 2 3;4 5 6;7 8 0; p=poly(a) p = 1.00 -6.00 -72.00 -27.00 r=roots(p) r = 1

22、2.12 -5.73 r是矩阵a的特征值 -0.39,Matlab程序设计,当然我们可用poly令其返回多项式形式 p2=poly(r) p2 = 1.00 -6.00 -72.00 -27.00 matlab规定多项式系数向量用行向量表示，一组根用列向量表示。,Matlab程序设计,7.3 conv多项式乘运算,例:a(x)=x2+2x+3; b(x)=4x2+5x+6; 求：c = (x2+2x+3)(4x2+5x+6) a=1 2 3;b=4 5 6; c1=conv(a,b) c2=conv(1 2 3,4 5 6) 卷积 c2 = 4 13 28 27 18 p=poly2str(c

23、,x) p = 4 x4 + 13 x3 + 28 x2 + 27 x + 18,Matlab程序设计,7.4 deconv多项式除运算,a=1 2 3; c = 4.00 13.00 28.00 27.00 18.00 d=deconv(c,a) 向量解卷积 d =4.00 5.00 6.00,Matlab程序设计,7.5 多项式微分,matlab提供了polyder函数多项式的微分。 命令格式： polyder(p): 求p的微分 polyder(a,b): 求多项式a,b乘积的微分 p,q=polyder(a,b): 求多项式a,b商的微分 例：a=1 2 3 4 5; poly2str(a,x) ans = x4 + 2 x3 + 3 x2 + 4 x + 5 b=polyder(a) b = 4 6 6 4 po