第四讲 MATLAB 语言 与现代科学计算

1. 第四讲 MATLAB语言 与现代科学计算 4.1 数值线性代数问题 4.2 数值微积分 4.3 数据插值与统计分析 4.4 数学问题的解析运算与高精度运算 4.5 例子与习题

2. 4.1 数值线性代数问题 一些 MATLAB 特殊矩阵 零矩阵：A=zeros(m, n);其中 （m, n） 定义零矩阵维数大小 全1矩阵：A=ones(m, n);其中 （m, n） 定义矩阵维数大小 单位矩阵：A=eye(m, n);其中 （m, n） 定义零矩阵维数大小 随机元素矩阵：A=rand(m, n);［0，1］上均匀分布 A=randn(m, n);正态分布 对角矩阵：A=diag(v); V 为对角向量

3. 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 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0

4. 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 = hankel(c,r) h = 1 2 3 8 2 3 8 9 3 8 9 10 p = [1 2 3 8 9 10]

5. 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 1 256 64 16 4 1 625 125 25 5 1

6. 4.1 数值线性代数问题 矩阵行列式： det(A) 矩阵的迹： trace(A) (即对角线元素之和） 矩阵的秩： rank(A) 即线性无关的列数或行数 矩阵的特征多项式： poly(A) 即 矩阵的特征根 roots(poly(A))

7. 4.1 数值线性代数问题 矩阵的范数 norm(A, p) 当A为向量时，范数定义稍有不同，如下表所示

8. 4.1 数值线性代数问题 矩阵的特征值与特征向量 ［V, D]= eig(A) 矩阵V的各列为特征向量，矩阵D的对角元素为特征值 >> 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

9. 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.0000 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

10. 4.1 数值线性代数问题 面向矩阵各个元素的函数 B= 函数名(A) 用命令 help elfun可以查看这些命令列表，主要有 sin sinh asin asinh cos cosh acos acosh tan tanh atan atan2 atanh sec sech 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

11. 4.1 数值线性代数问题 对矩阵进行数值分析的函数 B= 函数名(A) 用命令 help datafun可以查看这些命令列表，主要有 max min mean median std var sort sortrows sum 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

12. 数值差分运算 dy=diff(y) (按列运算） 4.2 数值微积分 >> y=magic(6) y = 35 1 6 26 19 24 3 32 7 21 23 25 31 9 2 22 27 20 8 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

13. 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 为被积分函数的调用次数。 tol 的缺省值对 quad 函数为 1e-3 对 quad8 和 quadl 函数为 1e-6 quad8 和 quadl 函数使用的算法比 quad 高得多

14. 4.2 数值微积分 例子：求无穷定积分 该无穷定积分的理论值为 1 >> 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.99999920879563 kk = 89 >> [y1,kk1]=quad8(f,-15,15) y1 = 0.99999999999999 kk1 = 769

15. 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

16. 4.3 数据插值与统计分析 一维插值用函数 y=interp1(x, y, x1, method) x 和 y 两个向量分别表示给定的一组自变量和函数值数据，x1 为一组新的插值点，而得出的 y1 是在这一组插值点处的插值结果，插值方法有 ‘linear’ 线性插值，‘cubic’三次插值，‘spline’ 样条插值等。 多项式拟合函数 p=polyfit(x, y,n) x 和 y 两个向量分别表示给定的一组自变量和函数值数据，n 为预期的多项式的阶次，返回的 p 为插值多项式系数。

17. 4.3 数据插值与统计分析 例：假设已知的数据点来自下面的函数 >> x=0:0.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); >> x1=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

18. 4.3 数据插值与统计分析

19. 4.3 数据插值与统计分析 数据统计分析 ［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) 求相关系数矩阵

20. 4.3 数据插值与统计分析 例子 >> s=randn(10000, 5); M=mean(s), D=std(s), MD=median(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.0079 0.0009 0.0081 -0.0028 1.0000

21. 4.4 数学问题的解析运算与高精度运算 MATLAB本身没有提供解析运算，即符号运算的能力，但借助于著名的符号运算软件 MAPLE， MATLAB提供了一个符号运算工具箱，使得MATAB可以进行符号运算及进行高精度运算 syms 变量列表 定义符号变量 s=sym(A)从A中定义符号对象 s。

22. 4.4 数学问题的解析运算与高精度运算 观察下例 >> a=hilb(5), b=sym(hilb(5)),whos a = 1.0000 0.5000 0.3333 0.2500 0.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 double array b 5x5 2478 sym object

23. 4.4 数学问题的解析运算与高精度运算 可以用命令 help symbolic查看MATLAB提供的符号运算命令，可以看出有很多，我们只介绍少数几个命令。 矩阵求逆 inv(A) >> 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/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]

24. 4.4 数学问题的解析运算与高精度运算 极限、微分、积分和 Taylor级数展开 limit, diff, int, taylor >> syms x a t h; limit(sin(x)/x) ans =1 >> limit((x-2)/(x^2-4),2) ans =1/4 >> limit((1+2*t/x)^(3*x),x,inf) ans =exp(6*t) >> limit(1/x,x,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]

25. 4.4 数学问题的解析运算与高精度运算 极限、微分、积分和 Taylor级数展开 limit, diff, int, taylor >> syms x t; diff(sin(x^2)) ans =2*cos(x^2)*x >> [diff(t^6,1), diff(t^6, 2), diff(t^6, 3),diff(t^6, 4),diff(t^6, 5),diff(t^6, 6)] ans =[ 6*t^5, 30*t^4, 120*t^3, 360*t^2, 720*t, 720] >> syms x x1 alpha u t; >> int(1/(1+x^2)) ans =atan(x) >> int(sin(alpha*u),alpha) ans =-1/u*cos(alpha*u) >> int(x1*log(1+x1),0,1) ans =1/4

26. 4.4 数学问题的解析运算与高精度运算 极限、微分、积分和 Taylor级数展开 limit, diff, int, taylor >> taylor(exp(-x)) ans =1-x+1/2*x^2-1/6*x^3+1/24*x^4-1/120*x^5 >> 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(x^t,3,t) ans =1+log(x)*t+1/2*log(x)^2*t^2

27. 4.4 数学问题的解析运算与高精度运算 无穷级数求和 symsum >> syms n k; symsum(2/2^n+2/3^n, n, 1, inf) ans =3 >> symsum(k, 0, n-1) ans =1/2*n^2-1/2*n >> symsum(k^2,0,10) ans =385 >> symsum(1/k^2,1,Inf) ans =1/6*pi^2

28. 4.4 数学问题的解析运算与高精度运算 积分变换 laplace, ilaplace, fourier, ifourier, ztrans, iztrans >> syms a s t w x >> laplace(t^5) ans =120/s^6 >> laplace(exp(a*s)) ans =1/(t-a) >> laplace(sin(w*x),t) ans =w/(t^2+w^2) >> syms t v w x >> fourier(exp(-x^2),x,t) ans =pi^(1/2)*exp(-1/4*t^2) >> fourier(exp(-t)*sym('Heaviside(t)'),v) ans =1/(1+i*v) >> fourier(diff(sym('F(x)')),x,w) ans =i*w*fourier(F(x),x,w)

29. 4.4 数学问题的解析运算与高精度运算 变精度运算 vpa(S, D)。函数 vpa(S, D)将符号结果 S 的数值答案求到任意指定的有效位数D， D的缺省值为16。 >> vpa(pi, 60) ans =3.14159265358979323846264338327950288419716939937510582097494 >> vpa(sqrt(2), 80) ans = 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070 >> vpa(hilb(2),30) ans = [ 1., .500000000000000000000000000000] [ .500000000000000000000000000000, .333333333333333333333333333333]

30. 4.5 例子与习题 例一、构造一个11阶的 Hilbert矩阵，并得出该矩阵的条件数，再用 inv 函数和 invhilb 函数求该矩阵的逆，并验证得出的逆矩阵是否真的满足逆矩阵的条件 例二、求解下面的线性代数方程，并验证得出的解满足原方程

31. 4.5 例子与习题 例三、求定积分 例四、调用 sort 函数按列对 magic(10)的结果矩阵进行正序排序与逆序排序，按行进行上述同样的排序 例五、演示 symbolic 工具箱的例子。