第五章_用Matlab进行现代科学运算

第五章用Matlab进行现代科学计算5.1解析解与数值解在实际的工程技术中，一般通过数值解法来获得问题的解。至少有两个原因：1、解析解不存在如圆周率的值就没有解析解，超越积分也没有解析解。在这种情况下，就必须采用数值解技术。22baxedx2、解析解存在但不实用如求解n元一次代数方程组的问题，从理论上讲，总可以把多元一次方程组简化成解析可解的形式。然而当n较大时，需要的基本运算次数是非常惊人的，对于求解实际问题所需要的计算成本是根本无法接受的，只能采用数值解法。5.2数值线性代数问题及求解5.2.1特殊矩阵的MATLAB输入零矩阵：zeros(m,n);(m,n)定义矩阵维数幺矩阵：ones(m,n);单位矩阵：eye(m,n);随机元素矩阵：rand(m,n);[0,1]上均匀分布randn(m,n);正态分布对角阵：diag(v);v为对角元素组成的向量5.2.2矩阵的特征参数运算矩阵的行列式：det(A)矩阵的迹：trace(A)，即对角线元素之和矩阵的秩：rank(A)，即线性无关的列数或行数矩阵的范数norm(A,p)11-范数，A的最大列和，max(sum(abs((A)))22-范数，最大奇异值，norm(A)Inf无穷范数，A的最大行和，max(sum(abs(A')))'fro'Frobenius范数，sqrt(sum(diag(A'*A)))当A为向量时，范数定义稍有不同：norm(A,p)sum(abs(A).^p)^(1/p)norm(A)norm(A,2)norm(A,inf)max(abs(A))norm(A,-inf)min(abs(A))p1矩阵的特征多项式：c=poly(A)111()detnnnnCssIAscscsc矩阵的特征值：特征方程的根roots(poly(A))或eig(A)若A为矩阵，则返回的c就是向量[1,c1,…,cn-1,cn]降幂排列矩阵的特征值与特征向量xxA[v,d]=eig(A)矩阵v的各列为特征向量，矩阵d的对角元素为特征值a=[123;456;780];[v,d]=eig(a)v=-0.2998-0.7471-0.2763-0.70750.6582-0.3884-0.6400-0.09310.8791d=12.1229000-0.3884000-5.7345矩阵求逆B=inv(A)为求矩阵A的逆矩阵a=[123;456;780];b=inv(a);c=a*b;d=b*a;[bcd]ans=-1.77780.8889-0.11111.00000-0.00001.00000.000001.5556-0.77780.2222-0.00001.00000-0.00001.00000-0.11110.2222-0.11110.0000-0.00001.00000.00000.00001.000伴随阵法，初等变换法。。。多项式及多项式矩阵求值•多项式求值B=polyval(aa,x)•多项式矩阵求值B=polyvalm(aa,A)•其中aa为多项式系数降幂排列构成的向量；x为指定的标量，A为给定的方阵面向矩阵各个元素的函数B=函数名(A)用命令helpelfun可以查看这些命令列表，主要有sinsinhasinasinhcoscoshacosacoshtantanhatanatan2atanhsecsechasecasechcsccschacscacschcotcothacotacothexploglog10log2pow2sqrtabsanglecomplexconjimagrealunwrapisrealcplxpairfixfloorceilroundmodremsign对矩阵进行数值分析的函数B=函数名(A)用命令helpdatafun可以查看这些命令列表，主要有maxminmeanmedianstdvarsortsortrowssumprodhisthistctrapzcumsumcumprodcumtrapzdiffgradientdel2corrcoefcovsubspacefilterfilter2convconv2convndeconvdetrendfftfft2fftnifftifft2ifftnfftshiftifftshift5.2.3矩阵的相似变换与分解1、三角形分解•最基本的分解方法。可把任何一个方阵表示成两个三角矩阵的乘积，其中一个方阵是换位的下三角矩阵，另一个是上三角矩阵。这种分解法常称为LU分解或LR分解，其算法多为高斯消去法。•是求解线性方程组、inv()和det()的基础。•用lu()函数可得到分解后的两个三角矩阵，其调用形式为[L,U]=lu(A)2、正交分解•也称为“QR”分解，将一个矩阵表示成一个正交矩阵和一个上三角矩阵的乘积，A=Q*R•用qr()函数可得到分解后的两个矩阵，其调用形式为[Q,R]=qr(A)3、奇异值分解[U,S,V]=svd(A)其中U和V是正交矩阵，而S是对角矩阵。即A=U*S*V’奇异值分解是矩阵分析的有力工具。Matlab的一些函数也是基于此实现的，包括求伪逆函数pinv(A)，求秩函数rank(A)，求欧式范数norm(A,2)，求条件数cond(A)5.3数值微积分问题5.3.1数值差分运算MATLAB提供了计算给定向量差分的函数diff()，其调用方法是dy=diff(y)(按列运算，行数减1）y=magic(6)y=351626192433272123253192222720828331710153053412141643629131811diff(y)ans=-32311-54128-23-514-5-231931-5-17-522-231-541-2631-514-55.3.2数值积分求解函数定积分的数值方法是多种多样的，如梯度法，Simpson法、Romberg法等都是经常采用的方法。基本思想是将整个积分空间分割成若干个自空间，而每个小空间上的函数积分可求，因而整个空间函数积分可求。[y,n]=sum(F,a,b);等宽矩形法求定积分[y,n]=trapz(F,a,b);梯形法求定积分[y,n]=quad(F,a,b,tol);自适应辛普森积分[y,n]=quadl(F,a,b,tol);自适应Lobatto积分其中F为描述被积函数的字符串变量，一般为一个F.m函数文件名（加引号）；还可用inline()来定义一个函数。a,b为积分上下限，tol为变步长积分用误差限，y为积分结果，n为被积分函数的调用次数。常见的一元函数数值积分指令：例子：求无穷定积分dxex2/221f=inline('1/sqrt(2*pi)*exp(-x.^2/2)','x');[y,kk]=quad(f,-8,8)y=1.00000197533430kk=81[y1,kk1]=quad8(f,-8,8)y1=1.00000000000023kk1=161[y,kk]=quad(f,-15,15)y=0.99999920879563kk=89[y1,kk1]=quad8(f,-15,15)y1=0.99999999999999kk1=769该无穷定积分的理论值为1例子：dxdyyxex22/1122sin212f=inline(‘1/(2*pi)*exp(-x.^2/2).*sin(x.^2+y)‘,'x','y');y=dblquad(f,-2,2,-1,1)y=0.2506用双重积分函数dblquad(F,x_m,x_M,y_m,y_M,tol)三重积分函数triplequad(F,x_m,x_M,y_m,y_M,z_m,z_M,tol)5.4常微分方程的数值解法5.4.1一般常微分方程的数值解法•通常把含有自变量t、未知的一元函数x(t)及其导数或微分的方程，叫做常微分方程（ODE）。如一阶常微分方程的初值问题：0(,),1,,,(0)iixftinxxx•求解常微分方程组的数值方法是多种多样的，如常用的Euler法、Runge-Kutta法、Adams线性多步法、Gear法等。MATLAB求解常微分方程的函数如下：函数说明函数说明ode23低阶法解非刚性微分方程ode45中阶法解非刚性微分方程ode113变阶法解非刚性微分方程ode23t梯形法解适度刚性ODE和DAEode15s变阶法解刚性微分方程和DAE（微分代数方程）ode23sode23tb低阶法解刚性微分方程ode23()和ode45()最常用，采用自适应变步长求解方法[t,x]=ode23(方程函数名,[t0,tf],x0,'选项')[t,x]=ode45(方程函数名,[t0,tf],x0,'选项')方程函数名为描述系统状态方程的M函数的名称，该函数名应该用引号括起来。t0和tf分别为起始和终止计算时间，x0为系统的初始状态变量的值。返回值为求解的时间变量和相应的状态列向量构成的矩阵转置方程函数名的编写格式是固定的：functionxdot=方程函数名(t,x)其中t为时间变量，x为方程的状态变量，xdot代表状态变量的导数。注意，及时微分方程是非时变的，也应该在函数输入变量列表中写上t占位。在lorenzeq.m函数中：functionxdot=lorenzeq(t,x)xdot=zeros(3,1);xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)]t_final=100;x0=[0;0;1e-10];[t,x]=ode45('lorenzeq',[0,t_final],x0);figure(1);set(gcf,'position',[7225390270]);plot(t,x),figure(2);set(gcf,'position',[405225390270]);plot3(x(:,1),x(:,2),x(:,3));Axis([1045-2020-2030])020406080100-40-20020406010203040-20020-20-100102030谢谢！