第三章MATLAB数值计算

1、12022-5-1第三章第三章 MATLAB数值计算数值计算22022-5-1何时需采用数值计算？n解析解不存在；n解析解存在，但非常复杂；如：由Cramer法则，可将n元一次方程组化简为n个n-1元一次方程； n-1元一次方程组可化简为n-1个n-2元一次方程；结论： n元一次方程组的解析解是可以求出的。但求解需进行方程(n-1)(n+1)!+n次基本运算。一个20元一次方程组需进行方程9.70731020次基本运算。常见的问题：力学中的有限元法；差分方程及微分方程的数值解法；FFT方法等。32022-5-1n3.1 线性方程组线性方程组n3.2 矩阵的非线性运算矩阵的非线性运算n3.3 数

2、值数值积分积分n3.4 数值微分数值微分n3.5 常微分方程常微分方程n3.6 MEX技术简介技术简介42022-5-1 Ax = b 有x = A-1b，但实际上并不显式求A-13.1 线性方程组例子： 7x = 21 x = 21/7=3如果求逆 x = 7-1 21 = .142857 21 = 2.99997这就需要一次除和一次乘，且精度更低52022-5-1运算符 AX = BX = AB 左除 XA = BX = B/A 右除62022-5-13x3的例子6475156230710321xxx65 5 46237710 32132121xxxxxxxx5 . 21 . 6755 .

3、 2061 . 000710321xxx1 . 65 . 2761 . 0055 . 200710321xxx2 . 65 . 272 . 60055 . 200710321xxx72022-5-1算法矩阵表示2 . 60055 . 200710U104. 03 . 0015 . 0001L010100001PPALU 单位下三角阵上三角阵置换阵82022-5-1LU分解例子5156230710A1000100011P104. 000100012M0101000012P105 . 0013 . 00011M103 . 0015 . 00011L104. 000100012L92022-5-1线

4、性方程组求解Ax = b LU = PA Ly = Pb Ux = y x = Ab例： clearA=1 2 3 4 5 6 7 8 1;L,U,P = lu(A)Err=P*A - L*U102022-5-1矩阵的QR分解：nQR分解也称为长方阵的正交分解，把长方阵分解为正交矩阵和上三角矩阵的积。n格式：Q,R=qr(A)例：X=1 2 3 4 5 6 7 8 2 1 6 3;Q,R =qr(X)Q1,R1 =qr(X)112022-5-1矩阵的Cholesky分解：nCholesky分解把对称正定矩阵分解为上三角矩阵和其转置的积，即X=CC，C为上三角阵。n格式：C=chol(A)例：X

5、=6 1 1 2 1 8 3 1 1 3 7 0 2 1 0 4;C=chol(X)C*C %验证X=CC122022-5-1分析在三位精度的算法中，x2是由两个与舍入误差相同阶的数相除得到的，因此x2可以是任意的数。然后完全任意的x2代入第一个方程得到x1 ，所以第一方程的残量是小的。我们也可以期望第二个方程的残量也是小的，因为矩阵是接近奇异的。部分主元的高斯消去法可以确保产生小的残量部分主元的高斯消去法可以确保产生小的残量“小的”舍入误差的阶与三个数量有关：原来矩阵悉数矩阵单元的大小，消去过程中系数矩阵单元的大小和计算解的元素的大小。如果其中一个是“大的”，残量在绝对意义上不必是小的。即使

6、残量很小，但不能保证解的误差很小。它们之间的关系可以部分由矩阵的条件数条件数来刻画132022-5-1范数如果Ax = b, 如何测量x关于A和b的敏感性？pxxpnipip1 ,/11iiniiniixxxxxxmax , ,2/121211 x =(1:4)/5 n1 = norm(x,1) n2 = norm(x) ninf = norm(x,inf)142022-5-1扰动分析bAx bbxxAxMb xmbbbAxx)(条件数是相对误差放大因子，也是刻画矩阵奇异性的测度。152022-5-1矩阵范数向量有三种常用范数，相对应的矩阵范数的三种形式为： （A的行范数） （A的列范数） （

7、 是ATA 的最大特征值） （ A的2范数）11maxniji niAa 111maxnijjniAa 12A1162022-5-1A=1 2 3 4 5 6 7 8 9;k2=norm(A) k2_1=norm(A,2)k1=norm(A,1) k_inf=norm(A,inf)A=1 2 3 4 5 6 7 8 1;k2=norm (A) k2_1=norm (A,2)k1=norm (A,1) k_inf=norm (A,inf)例：172022-5-1矩阵条件数xAxmxAxMmin ,max条件数：1minmax)(AAxAxxAxa182022-5-1矩阵条件数 cond(A) 或

8、cond(A,2):使用svd(A) cond(A,1)计算k1(A)：使用inv(A)，工作量比 cond(A,2)小 cond(A,inf): 同cond(A,1) condest(A) 估计k1(A)：用lu(A)用Higham和Tisseur的算法，适合大的稀疏矩阵192022-5-1A=1 2 3 4 5 6 7 8 9;k2=cond(A) k2_1=cond(A,2)k1=cond(A,1) k_inf=cond(A,inf)A=1 2 3 4 5 6 7 8 1;k2=cond(A) k2_1=cond(A,2)k1=cond(A,1) k_inf=cond(A,inf)例：2

9、02022-5-1非齐次线性方程组的求解n对于非齐次线性方程组Ax=b而言，要根据系数矩阵A的秩和增广矩阵B=A b的秩和未知数个数n的关系，才能判断方程组Ax=b的解的情况；1如果系数矩阵的秩=增广矩阵的秩=n，则方程组有唯一解；2如果系数矩阵的秩=增广矩阵的秩n，则方程组有无穷多解；3如果系数矩阵的秩增广矩阵的秩，则方程组无解。212022-5-1n对于非齐次线性方程组Ax=b而言，首先应判断方程组解的情况，其次，若有解，先求出方程组的特解，再次，求出对应齐次方程组的通解，最后，写出非齐次方程组的通解，即特解+对应齐次方程组的通解。n求Ax=b对应的齐次方程组Ax=0的通解，可以利用函数n

10、ull或对系数矩阵A施行行变换；n求Ax=b的特解，方法较多，根据方程组中方程的个数m和未知数的个数n，可以把方程组Ax=b分为： 超定方程组（mn）； 恰定方程组（m=n）； 欠定方程组（mn）。可以根据系数矩阵A的性质选用适当的计算方法，如可利用左除左除 “” 、系数矩阵A的逆或伪逆、或利用lu分解、或cholesky分解等。222022-5-1clearA=4 2 -1;3 -1 2;11 3 0; b=2 10 8; B=A b; %增广矩阵 n=3; RA=rank(A) RB=rank(B)if (RA=RB & RA=n) X=Ab else if (RA=RB &

11、; RA *Inner matrix dimensions must agree.372022-5-12牛顿柯特斯法牛顿柯特斯法基于牛顿柯特斯法，基于牛顿柯特斯法，MATLAB给出了给出了quadl函数来求定积函数来求定积分。分。（在（在6.5以前的版本中为：以前的版本中为： quad8 ） 该函数的调用格式为：该函数的调用格式为：I,n=quadl(fname,a,b,tol,trace)其中参数的含义和其中参数的含义和quad函数相似，只是函数相似，只是tol的缺省值取的缺省值取10-6。该函数可以更精确地求出定积分的值，且一般情况下函数该函数可以更精确地求出定积分的值，且一般情况下函数调

12、用的步数明显小于调用的步数明显小于quad函数，从而保证能以更高的效函数，从而保证能以更高的效率求出所需的定积分值。率求出所需的定积分值。382022-5-1例例3-2 求定积分。求定积分。(1) 被积函数文件被积函数文件fx.m。 function f=fx(x) f=x.*sin(x)./(1+cos(x).*cos(x);(2) 调用函数调用函数quadl求定积分。求定积分。 I=quadl(fx,0,pi) I = 2.4674注：在6.5以上的版本中将提示：Warning: QUAD8 is obsolete. We use QUADL instead.392022-5-1例例3-3

13、 分别用分别用quad函数和函数和quadl函数求定积分的近似函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。值，并在相同的积分精度下，比较函数的调用次数。调用函数调用函数quad求定积分：求定积分：clearfx=inline(exp(-x) ); I,n=quad(fx,1,2.5,1e-10)I = 0.28579444254766n = 65构造构造inline函数，函数，可省略不必要的函可省略不必要的函数文件数文件402022-5-1 调用函数调用函数quadl求定积分：求定积分：clear fx=inline(exp(-x);I,n=quadl(fx,1,2.5,1

14、e-10)I=vpa(I,10)输出：输出：I = 0.2858n = 18I =.2857944425412022-5-13被积函数由一个表格定义被积函数由一个表格定义在在MATLAB中，对由表格形式定义的函数关系的求定积分中，对由表格形式定义的函数关系的求定积分问题用问题用trapz(X,Y)函数。其中向量函数。其中向量X,Y定义函数关系定义函数关系Y=f(X)。例：用例：用trapz函数（梯形法函数（梯形法trapezoidal ）计算定积分。）计算定积分。命令如下：命令如下：clearX=1:0.01:2.5;Y=exp(-X); %生成函数关系数据向量生成函数关系数据向量trapz(

15、X,Y)vpa(ans,11)ans = 0.2858ans = .28579682416422022-5-13.3.3 二重定积分的数值求解二重定积分的数值求解使用使用MATLAB提供的提供的dblquad函数就可以直函数就可以直接求出上述二重定积分的数值解。该函数接求出上述二重定积分的数值解。该函数的调用格式为：的调用格式为：I=dblquad(f,a,b,c,d,tol,trace)该函数求该函数求f(x,y)在在a,bc,d区域上的二重定区域上的二重定积分。参数积分。参数tol，trace的用法与函数的用法与函数quad完完全相同。全相同。432022-5-1例例3-5 计算二重定积分

16、计算二重定积分(1) 建立一个函数文件建立一个函数文件fxy.m：function f=fxy(x,y)global ki;ki=ki+1; %ki用于统计被积函数的调用次数用于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);(2) 调用调用dblquad函数求解。函数求解。global ki; ki=0;I=dblquad(fxy,-2,2,-1,1) %单引号可有可无单引号可有可无kiI = 1.57449318974494ki = 1038442022-5-1clearfxy=inline(exp(-x.2/2).*sin(x.2+y);I=dblquad(fx

17、y,-2,2,-1,1)运用运用inline函数，更简洁函数，更简洁452022-5-13.4 数值微分数值微分3.4.1 数值微分的实现数值微分的实现在在MATLAB中，没有直接提供求数值导数的函数，中，没有直接提供求数值导数的函数，只有计算向前差分的函数只有计算向前差分的函数diff，其调用格式为：，其调用格式为：DX=diff(X)：计算向量计算向量X的向前差分，的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1。DX=diff(X,n)：计算计算X的的n阶向前差分。例如，阶向前差分。例如， diff(X,2)=diff(diff(X)。DX=diff(A,n,dim)：

18、计算矩阵计算矩阵A的的n阶差分，阶差分，dim=1时时(缺缺省状态省状态)，按列计算差分；，按列计算差分；dim=2，按行计算差分。，按行计算差分。462022-5-1例例3-6 用不同的方法求函数用不同的方法求函数f(x)在在-3，3的数值导数，并在的数值导数，并在同一个坐标系中做出同一个坐标系中做出f(x)的图像。的图像。f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);132625326212(5)523411

19、522126(5)fxxxxxdfxxgdxxxxx472022-5-1例例3-6 程序如下：程序如下：clearf=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);x=-3:0.01:3;p=polyfit(x,f(x),5); %用用5次多项式次多项式p拟合拟合f(x)dp=polyder(p); %对拟合多项式对拟合多项式p求导数求导数dpdpx=polyval(dp,x); %求求dp在假设点的函数值在假设点的

20、函数值dx=diff(f(x,3.01)/0.01; %直接对直接对f(x)求数值导数求数值导数gx=g(x); %求函数求函数f的导函数的导函数g在假设点的导数在假设点的导数plot(x,dpx, k.,x,dx,.,x,gx, r-); %作图作图482022-5-13.5 常微分方程初值问题的数值解法常微分方程初值问题的数值解法3.5.1 龙格库塔法简介龙格库塔法简介Runge-Kutta方法：1905（德国）httsssshyyhsyhtfsshyhtfsshyhtfsytfsnnnnnnnnnnnn1432113423121)22(6),()2,2()2,2(),(若f不依赖y,则为

21、simpson公式步长 hRe(max)492022-5-13.5.2 龙格库塔法的实现龙格库塔法的实现 基于龙格库塔法，基于龙格库塔法，MATLAB提供了求提供了求非刚性非刚性常微分方程数值解的函常微分方程数值解的函数，一般调用格式为：数，一般调用格式为： t,y=ode23(fname,tspan,y0) t,y=ode45(fname,tspan,y0)其中其中fname是定义是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。的函数文件名，该函数文件必须返回一个列向量。tspan形式为形式为t0,tf,表示求解区间。表示求解区间。y0是初始状态列向量。是初始状态列向量。t和和

22、y分别给分别给出时间向量和相应的状态向量。出时间向量和相应的状态向量。刚性刚性常微分方程数值解的函数，一般调用格式为：常微分方程数值解的函数，一般调用格式为：ode15s及及ode23s502022-5-1例例3-10 试求初值问题的数值解，并与精确解相比较。试求初值问题的数值解，并与精确解相比较。 (1) 建立函数文件建立函数文件funt.m。function yp=funt(t,y)yp=(y2-t-2)/4/(t+1);(2) 求解微分方程。求解微分方程。cleart0=0; tf=10; y0=2;t,y=ode45(funt,t0,tf,y0); %求数值解求数值解y1=sqrt(t

23、+1)+1; %求精确解求精确解diary e:L.txt %生成日记文件生成日记文件t %输出行向量输出行向量yy1 diary off % y为数值解，为数值解，y1为精确值，显然两者近似。为精确值，显然两者近似。224(1)ytyt 512022-5-1nload-加载数据文件a.dat或b.txt中的数据 例：振动测试的功率谱分析，powerp.mntype-在workspace中显示a.dat或b.txt中的数据 例：type data.dat522022-5-1例例3-11 求解著名的求解著名的Van der Pol方程。方程。2(1)0yyyy变形：令变形：令12,xyxy122

24、2121(1)xxxxxx 532022-5-1function y=vdp_eq(t,x) mu=1;y=x(2); -mu*(x(1).2-1).*x(2)-x(1)clearx0=-0.2; -0.7; t_final=20;t,y=ode45(vdp_eq,0,t_final,x0);plot(t,y)figure; plot(y (:,1),y(:,2), -)1222121(1)xxxxxx 注意自变量必须写成注意自变量必须写成x(1)和和 x(2)542022-5-1例例3-12 Duffing方程，求时间相应，试绘制系方程，求时间相应，试绘制系统相平面图。统相平面图。30 xx

25、xx1232121xxxxxx 552022-5-1function y=duff_eq(t,x) mu=3; alpha=.3; k=-1;y=x(2); -k*x(1)-alpha*x(2)-mu*x(1).3;clear; closex0=-1;1; t_final=50;t,y=ode45(duff_eq,0,t_final,x0);plot(t,y);gridfigure; plot(y (:,1),y(:,2),:)1232121xxxxxx 562022-5-1n例例3-13 Lorenz系统，求时间相应，试绘制系统系统，求时间相应，试绘制系统相平面图。相平面图。()xxyyxz

26、rxyzxybz 1083b572022-5-1Lorenz系统方程：系统方程：function y=lorenzeq(t,x) sigma=10; b=8/3; r=28;y=-sigma*(x(1)-x(2); -x(1)*x(3)+r*x(1)-x(2); x(1)*x(2)-b*x(3);clear; close %混沌现象t_final=100; x0=0.01;0.01;0;t,x=ode45(lorenzeq,0,t_final,x0);plot(t,x),figure; plot3(x(:,1),x(:,2),x(:,3); grid()xxyyxzrxyzxybz 582022-5-13.6 M