MATLAB语言与现代科学计算.ppt

1、 电子发烧友 电子技术论坛,MATLAB程序语言设计,应自炉 五邑大学信息学院2007年, , 电子发烧友 电子技术论坛,第四讲 MATLAB语言 与现代科学计算,4.1 数值线性代数问题 4.2 数值微积分 4.3 数据插值与统计分析 4.4 数学问题的解析运算与高精度运算 4.5 例子与习题, 电子发烧友 电子技术论坛,4.1 数值线性代数问题,一些 MATLAB 特殊矩阵,零矩阵： A=zeros(m, n); 其中 （m, n） 定义零矩阵维数大小,全1矩阵： A=ones(m, n); 其中 （m, n） 定义矩阵维数大小,单位矩阵： A=eye(m, n); 其中 （m, n） 定

2、义零矩阵维数大小,随机元素矩阵： A=rand(m, n); 0，1上均匀分布 A=randn(m, n); 正态分布,对角矩阵： A=diag(v); V 为对角向量, 电子发烧友 电子技术论坛,4.1 数值线性代数问题,Hilbert 矩阵：A=hilb(n);,The elements of the Hilbert matrices are: H(i, j) = 1/(i+j-1),伴随矩阵：A=compan(p); 其中 p 为多项式系数向量,多项式 对应的向量为 p= 1 2 7 6 9 8 ，它的伴随矩阵为 compan(p) ans = -2 -7 -6 -9 -8 1 0 0

3、0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0, 电子发烧友 电子技术论坛,4.1 数值线性代数问题,Hankel 矩阵：A=hankel(c, r);,A Hankel matrix is a matrix that is symmetric and constant across the anti-diagonals, and has elements h(i,j) = p(i+j-1), where vector p = c r(2:end) completely determines the Hankel matrix.,c = 1:3; r = 7:10; h =

4、 hankel(c,r) h = 1 2 3 8 2 3 8 9 3 8 9 10 p = 1 2 3 8 9 10, 电子发烧友 电子技术论坛,4.1 数值线性代数问题,Vandermonde 矩阵：A=vander(c, r);,A = vander(v) returns the Vandermonde matrix whose columns are powers of the vector V, that is A(i,j) = v(i)(n-j). c=1:5, v=vander(c) c = 1 2 3 4 5 v = 1 1 1 1 1 16 8 4 2 1 81 27 9 3

5、1 256 64 16 4 1 625 125 25 5 1, 电子发烧友 电子技术论坛,4.1 数值线性代数问题,矩阵行列式： det(A),矩阵的迹： trace(A) (即对角线元素之和）,矩阵的秩： rank(A) 即线性无关的列数或行数,矩阵的特征多项式： poly(A) 即,矩阵的特征根 roots(poly(A), 电子发烧友 电子技术论坛,4.1 数值线性代数问题,矩阵的范数 norm(A, p),当A为向量时，范数定义稍有不同，如下表所示, 电子发烧友 电子技术论坛,4.1 数值线性代数问题,矩阵的特征值与特征向量,V, D= eig(A) 矩阵V的各列为特征向量，矩阵D的对

6、角元素为特征值, a=1 2 3;4 5 6;7 8 0; v,d=eig(a) v = -0.2998 -0.7471 -0.2763 -0.7075 0.6582 -0.3884 -0.6400 -0.0931 0.8791 d = 12.1229 0 0 0 -0.3884 0 0 0 -5.7345, 电子发烧友 电子技术论坛,4.1 数值线性代数问题,矩阵求逆与线性方程求解,B= inv(A) 求矩阵A的逆矩阵, a=1 2 3;4 5 6; 7 8 0; b=inv(a); c=a*b;d= b*a; b c d ans = -1.7778 0.8889 -0.1111 1.000

7、0 0 -0.0000 1.0000 0.0000 0 1.5556 -0.7778 0.2222 -0.0000 1.0000 0 -0.0000 1.0000 0 -0.1111 0.2222 -0.1111 0.0000 -0.0000 1.0000 0.0000 0.0000 1.0000, 电子发烧友 电子技术论坛,4.1 数值线性代数问题,面向矩阵各个元素的函数,B= 函数名(A),用命令 help elfun 可以查看这些命令列表，主要有,sin sinh asin asinh cos cosh acos acosh tan tanh atan atan2 atanh sec s

8、ech asec asech csc csch acsc acsch cot coth acot acoth exp log log10 log2 pow2 sqrt abs angle complex conj imag real unwrap isreal cplxpair fix floor ceil round mod rem sign, 电子发烧友 电子技术论坛,4.1 数值线性代数问题,对矩阵进行数值分析的函数,B= 函数名(A),用命令 help datafun 可以查看这些命令列表，主要有,max min mean median std var sort sortrows su

9、m prod hist histc trapz cumsum cumprod cumtrapz diff gradient del2 corrcoef cov subspace filter filter2 conv conv2 convn deconv detrend fft fft2 fftn ifft ifft2 ifftn fftshift ifftshift, 电子发烧友 电子技术论坛,4.2 数值微积分,数值差分运算 dy=diff(y) (按列运算）, y=magic(6) y = 35 1 6 26 19 24 3 32 7 21 23 25 31 9 2 22 27 20 8

10、 28 33 17 10 15 30 5 34 12 14 16 4 36 29 13 18 11 diff(y) ans = -32 31 1 -5 4 1 28 -23 -5 1 4 -5 -23 19 31 -5 -17 -5 22 -23 1 -5 4 1 -26 31 -5 1 4 -5, 电子发烧友 电子技术论坛,4.2 数值微积分,数值积分 y, n=quad( F, a, b, tol); y, n=quad8(F, a, b, tol); y, n=quadl(F, a, b, tol);,其中F为被积函数， a, b为积分上下限， tol 为误差限 y 为积分结果，n 为被

11、积分函数的调用次数。,tol 的缺省值对 quad 函数为 1e-3 对 quad8 和 quadl 函数为 1e-6 quad8 和 quadl 函数使用的算法比 quad 高得多, 电子发烧友 电子技术论坛,4.2 数值微积分,例子：求无穷定积分, f=inline(1/sqrt(2*pi)*exp(-x.2/2),x); y,kk=quad(f,-8,8) y = 1.00000197533430 kk = 81 y1,kk1=quad8(f,-8,8) y1 = 1.00000000000023 kk1 = 161, y,kk=quad(f,-15,15) y = 0.99999920

12、879563 kk = 89 y1,kk1=quad8(f,-15,15) y1 = 0.99999999999999 kk1 = 769,该无穷定积分的理论值为 1, 电子发烧友 电子技术论坛,4.2 数值微积分,例子：双重积分,用函数 dblquad(F, xmin, xmax, ymin, ymax, tol) 求双重积分, f=inline(exp(-x.2/2).*sin(x.2+y),x,y); y=dblquad(f,-2,2,-1,1) y = 1.57450259485169, 电子发烧友 电子技术论坛,4.3 数据插值与统计分析,一维插值用函数 y=interp1(x, y

13、, x1, method) x 和 y 两个向量分别表示给定的一组自变量和函数值数据，x1 为一组新的插值点，而得出的 y1 是在这一组插值点处的插值结果，插值方法有 linear 线性插值，cubic三次插值，spline 样条插值等。,多项式拟合函数 p=polyfit(x, y,n) x 和 y 两个向量分别表示给定的一组自变量和函数值数据，n 为预期的多项式的阶次，返回的 p 为插值多项式系数。, 电子发烧友 电子技术论坛,4.3 数据插值与统计分析,例：假设已知的数据点来自下面的函数, x=0:0.12:1; y=(x.2-3*x+5).*exp(-5*x).*sin(x); x1=

14、0:0.02:1; y0=(x1.2-3*x1+5).*exp(-5*x1).*sin(x1); y1=interp1(x,y,x1); y2=interp1(x,y,x1,cubic); y3=interp1(x,y,x1,spline); plot(x,y,x,y,o); hold on; plot(x1,y1,y2,y3,:,x,y,o,x1,y0) max(abs(y0-y2), max(abs(y0-y3) ans = 0.01765138386097 0.00861395506624, 电子发烧友 电子技术论坛,4.3 数据插值与统计分析, 电子发烧友 电子技术论坛,4.3 数据插

15、值与统计分析,数据统计分析,x, i= max(X), 或 x , i= min(X) 返回矩阵 X 各列的最大值或最小值， i 为各列最大值或最小值所在位置的行号构成的向量,x, i= sort(X) 对矩阵X各列排序,y= mean(X) 求矩阵X中各列向量的均值,y= std(X) 求矩阵X中各列向量的标准差,y= median(X) 求矩阵X中各列向量的中间值,y= cov(X) 求协方差,y= corrcoef(X) 求相关系数矩阵, 电子发烧友 电子技术论坛,4.3 数据插值与统计分析,例子, s=randn(10000, 5); M=mean(s), D=std(s), MD=m

16、edian(s) M = 0.0148 -0.0149 -0.0026 -0.0067 -0.0037 D = 0.9940 1.0024 1.0098 1.0044 1.0043 MD = 0.0250 -0.0104 0.0064 -0.0090 -0.0185 V=corrcoef(s) V = 1.0000 0.0047 0.0100 0.0197 -0.0079 0.0047 1.0000 0.0219 -0.0062 0.0009 0.0100 0.0219 1.0000 0.0021 0.0081 0.0197 -0.0062 0.0021 1.0000 -0.0028 -0.0

17、079 0.0009 0.0081 -0.0028 1.0000, 电子发烧友 电子技术论坛,4.4 数学问题的解析运算与高精度运算,MATLAB本身没有提供解析运算，即符号运算的能力，但借助于著名的符号运算软件 MAPLE， MATLAB提供了一个符号运算工具箱，使得MATAB可以进行符号运算及进行高精度运算,syms 变量列表 定义符号变量,s=sym(A) 从A中定义符号对象 s。, 电子发烧友 电子技术论坛,4.4 数学问题的解析运算与高精度运算,观察下例, a=hilb(5), b=sym(hilb(5),whos a = 1.0000 0.5000 0.3333 0.2500 0.

18、2000 0.5000 0.3333 0.2500 0.2000 0.1667 0.3333 0.2500 0.2000 0.1667 0.1429 0.2500 0.2000 0.1667 0.1429 0.1250 0.2000 0.1667 0.1429 0.1250 0.1111 b = 1, 1/2, 1/3, 1/4, 1/5 1/2, 1/3, 1/4, 1/5, 1/6 1/3, 1/4, 1/5, 1/6, 1/7 1/4, 1/5, 1/6, 1/7, 1/8 1/5, 1/6, 1/7, 1/8, 1/9 Name Size Bytes Class a 5x5 200 d

19、ouble array b 5x5 2478 sym object, 电子发烧友 电子技术论坛,4.4 数学问题的解析运算与高精度运算,可以用命令 help symbolic 查看MATLAB提供的符号运算命令，可以看出有很多，我们只介绍少数几个命令。, a=sym(hilb(4); b=inv(a); a b ans = 1, 1/2, 1/3, 1/4, 16, -120, 240, -140 1/2, 1/3, 1/4, 1/5, -120, 1200, -2700, 1680 1/3, 1/4, 1/5, 1/6, 240, -2700, 6480, -4200 1/4, 1/5, 1

20、/6, 1/7, -140, 1680, -4200, 2800 c=a*b c = 1, 0, 0, 0 0, 1, 0, 0 0, 0, 1, 0 0, 0, 0, 1,矩阵求逆 inv(A), 电子发烧友 电子技术论坛,4.4 数学问题的解析运算与高精度运算,极限、微分、积分和 Taylor级数展开 limit, diff, int, taylor, syms x a t h; limit(sin(x)/x) ans =1 limit(x-2)/(x2-4),2) ans =1/4 limit(1+2*t/x)(3*x),x,inf) ans =exp(6*t) limit(1/x,x,

21、0,right) ans =inf limit(1/x,x,0,left) ans =-inf limit(sin(x+h)-sin(x)/h,h,0) ans =cos(x) v = (1 + a/x)x, exp(-x);limit(v,x,inf,left) ans = exp(a), 0, 电子发烧友 电子技术论坛,4.4 数学问题的解析运算与高精度运算,极限、微分、积分和 Taylor级数展开 limit, diff, int, taylor, syms x t; diff(sin(x2) ans =2*cos(x2)*x diff(t6,1), diff(t6, 2), diff(

22、t6, 3),diff(t6, 4),diff(t6, 5),diff(t6, 6) ans = 6*t5, 30*t4, 120*t3, 360*t2, 720*t, 720, syms x x1 alpha u t; int(1/(1+x2) ans =atan(x) int(sin(alpha*u),alpha) ans =-1/u*cos(alpha*u) int(x1*log(1+x1),0,1) ans =1/4, 电子发烧友 电子技术论坛,4.4 数学问题的解析运算与高精度运算,极限、微分、积分和 Taylor级数展开 limit, diff, int, taylor, tayl

23、or(exp(-x) ans =1-x+1/2*x2-1/6*x3+1/24*x4-1/120*x5 taylor(log(x),6,1) ans =x-1-1/2*(x-1)2+1/3*(x-1)3-1/4*(x-1)4+1/5*(x-1)5 taylor(sin(x),pi/2,6) ans =1-1/2*(x-1/2*pi)2+1/24*(x-1/2*pi)4 taylor(xt,3,t) ans =1+log(x)*t+1/2*log(x)2*t2, 电子发烧友 电子技术论坛,4.4 数学问题的解析运算与高精度运算,无穷级数求和 symsum, syms n k; symsum(2/2

24、n+2/3n, n, 1, inf) ans =3 symsum(k, 0, n-1) ans =1/2*n2-1/2*n symsum(k2,0,10) ans =385 symsum(1/k2,1,Inf) ans =1/6*pi2, 电子发烧友 电子技术论坛,4.4 数学问题的解析运算与高精度运算,积分变换 laplace, ilaplace, fourier, ifourier, ztrans, iztrans, syms a s t w x laplace(t5) ans =120/s6 laplace(exp(a*s) ans =1/(t-a) laplace(sin(w*x),t) ans =w/(t2+w2), syms t v w x fo