《matlab数值计算数值计算》由会员分享，可在线阅读，更多相关《matlab数值计算数值计算（63页珍藏版）》请在金锄头文库上搜索。
1、数学建模暑期培训班,Matlab 初步,滨州学院数学与信息科学系,王 磊,,Matlab,第4讲 基本的数值计算方法,基本的数值计算问题,数值微分、插值与数值积分 线性方程组的求解(包括矩阵的各种相关计算、线性拟合） 多项式的处理 非线性方程（组）的求解 微分方程（组）的求解,数值微分(导数),用离散方法近似计算函数y=f(x)在某点x=a的导数值,前差公式,后差公式,中点公式,Matlab的（数值）差分计算程序：y=diff(x) 因此，可计算数值导数=diff(y)./diff(x) 如果是符号计算，diff(y,x)直接得到y对x的导函数,插值问题,g表达式复杂,甚至无表达式,1.分段线
2、性插值,实用插值方法,2. 三次样条插值,细木条：样条,输入:节点x0,y0, 插值点x (均为数组，长度自定义)； 输出:插值y (与x同长度数组)。,1. 分段线性插值:已有程序 y=interp1(x0,y0,x) y=interp1(x0,y0,x,linear),2. 三次样条插值:已有程序 y=interp1(x0,y0,x,spline) 或 y=spline(x0,y0,x),用MATLAB作插值计算,注：MATLAB有样条工具箱（Spline Toolbox）,插值的应用,加工时需要x每改变0.05时的y值,chazhi2,图1 零件的轮廓线 （x间隔0.2）,表1 x间隔0
3、.2的加工坐标x,y（图1右半部的数据）,数控机床加工零件,模型 将图1逆时针方向转90度，轮廓线上下对称，只需对上半部计算一个函数在插值点的值。,图2 逆时针方向转90度的结果,令v=x, u= -y,用MATLAB 作数值积分,trapz(x),输入数组x,输出按梯形公式x的积分(单位步长),trapz(x,y),输入同长度数组 x,y,输出按梯形公式 y对x的积分(步长不一定相等),n充分大时In就是I的数值积分,回 忆 定 积 分 的 定 义,用MATLAB 作数值积分,辛普森公式,quad(fun,a,b,tol,trace) I,fn=quad(),用自适应辛普森公式计算 tol为
4、绝对误差，缺省时为10-6,Gauss-Lobatto公式,quadl(fun,a,b,tol,trace) I,fn=quadl(),用自适应Gauss-Lobatto公式计算 tol为绝对误差，缺省时为10-6,注意：fun.m中应以自变量为矩阵的形式输入(点运算),矩形域上计算二重积分的命令：,dblquad(fun,xmin,xmax,ymin,ymax,tol),广义积分、二重和三重积分,长方体上计算三重积分的命令：,triplequad(fun,xmin,xmax,ymin,ymax, zmin,zmax,tol),注：fun是被积函数，本身可以有自己的参数,广义积分：,通过分析和
5、控制误差，转换成普通积分,quadv(fun,a,b,tol,trace),向量值积分：,用MATLAB 作数值积分,例. 计算,1）矩形公式和梯形公式,将(0, /4)100等分,2）辛普森公式和Gauss-Lobatto公式,精确、方便,无法计算用数值给出的函数的积分,Jifen1a.m Jifen1b.m,精确值为,用随机模拟计算数值积分,定积分的计算 重积分的计算 MATLAB实现,方法的直观解释随机投石,向单位正方形里随机投n块小石头,*,*,*,*,*,*,*,若有k块小石头落在1/4单位圆内，当n很大时,1/4单位圆的面积,（计算的一种方法）,1）随机投点法,目的：计算1/4单位
6、圆的面积,大数定律（贝努利定理）,随机变量(X,Y)在单位正方形内均匀分布,点(xi, yi)落在1/4单位圆内概率,一般地,投点坐标(xi, yi), xi, yi是相互独立、(0,1)内均匀分布的随机变量（(0,1)随机数）,设k是n次独立重复试验中事件A发生的次数。p是事件A在每次试验中发生的概率，则对任意的正数, 有,产生 n组(0,1)随机数 (xi,yi),其中 k 组满足,随机投点法(续),大数定律（辛钦定理） 设随机变量,相互独立，服从同一个分布，且具有数学期望,则对任意的正数有,随机变量 X 的概率密度为,的期望为,2) 均值估计法,产生(0,1)随机数xi (i=1,2,n
7、), n很大,用随机模拟方法计算任意区间上的积分,其中ui为(0,1)随机数,不要产生yi, 不用比较,均值估计法 的优点,均值估计法(续),随机模拟法计算重积分,产生相互独立(0,1)随机数 xi,yi,i =1, n；,落在内 m个点记作(xk,yk),k=1, m,可用于任意的 f, , 且可推广至高维,结果的精度和收敛速度与维数无关,计算量大，精度低，结果具有随机性,一般区间重积分的计算,分别为a, b和c, d 区间上的均匀分布随机数，判断每个点是否落在域内，将落在域内的m个点记作,则,MATLAB实现,随机数的产生：unifrnd(a,b,m,n),产生m行n列a,b区间上的均匀分
8、布随机数。当a=0, b=1时，可用 rand(m,n),随机投点法计算 n=10000; x=rand(2,n); k=0; for i=1:n if x(1,i)2+x(2,i)2=1 k=k+1; end end p=4*k/n,积分域和被积函数的对称性,蒙特卡罗方法: x取(0,a)随机数, y取(0,b)随机数,例:炮弹命中概率,1是椭圆在第1象限的部分,线性方程组的一般形式、两类解法,直接法 经过有限次算术运算求出精确解（实际上由于有舍入误差只能得到近似解）- 高斯（Gauss）消元法及与它密切相关的矩阵LU分解等 迭代法 从初始解出发，根据设计好的步骤用逐次求出的近似解逼近精确解
9、 - 雅可比（Jacobi）迭代法和高斯塞德尔（GaussSeidel）迭代法等,线性方程组数值解法的MATLAB实现,若A为可逆方阵，输出原方程的唯一解x,若A为nm矩阵（nm），且ATA可逆，输出原方程的最小二乘解x (即“线性拟合”),求解xA=b, 用右除：x=b/A,直接法一般用LU分解：LUx=b Ux=y; Ly=b,注：如果A的条件数 很大(病态)，将给出警告信息,注：除法计算不仅比 x=inv(A)*b 快，而且精度高,线性方程组数值解法的MATLAB实现,x,y,p=lu(A) 若A可逆，输出x为单位下三角阵L， y为上三角阵U，p为一交换阵P，使PA=LU.,u =cho
10、l(A) 对正定对称矩阵A的Cholesky分解，输出u为上三角阵U，使A=UTU,2. 矩阵LU分解,x,y=lu(A) 若A可逆且顺序主子式不为零， 输出x为单位下三角阵L，y为上三角阵U，使A=LU； 若A可逆，x为一交换阵与单位下三角阵之积.,例. 解,A=10 3 1;2 -10 3;1 3 10, b=14 -5 14, x=Ab, L1,U1=lu(A); L1,U1, A1=L1*U1, L2,U2,P=lu(A); L2,U2,P, A2=L2*U2, A3=inv(P)*A2,并对系数矩阵作LU分解,shiyan51,若第1个方程改为 3x2+x3=14 结果如何,向量和矩
11、阵的范数,度量向量、矩阵大小的数量指标,向量范数,矩阵范数,向量和矩阵范数的相容性条件,线性方程组数值解法的MATLAB实现,3. 范数 条件数 特征值 n=norm(x) 输入x为向量或矩阵，输出为x的2-范数 c=cond(x) 输入x为矩阵, 输出为x的2-条件数 r=rcond(x) 输入x为方阵, 输出为x条件数倒数 e=eig(x) 输入x为矩阵，输出x的全部特征值,方阵x的条件数定义： cond(x)=norm(x)* norm(inv(x) V,D=eig(x) 输入x为矩阵，输出x的特征向量阵V 和特征值对角阵D,使得x*V=V*D,当n很大时Hilbert矩阵呈病态,线性方
12、程组数值解法的MATLAB实现,H=hilb(5), h=rats(H), b=ones(5,1); x=Hb; b(5)=1.1; x1=Hb; x,x1, n1=cond(H), n2=rcond(H),例：观察Hilbert矩阵的病态性,例. Hx=b, 其中 H=hilb(5), b=1,1T,shiyan52,cond(H)=4.7661e+005,1. 提取（产生）对角阵,v=diag(x) 输入向量x，输出v是以x为对角元素的对角阵；输入矩阵x，输出v是x的对角元素构成的向量；,例：v=diag(diag(x) 输入矩阵x，输出v是x的对角元素构成的对角阵，可用于迭代法中从A中提
13、取D。,2. 提取上（下）三角阵,其他相关的MATLAB函数,y=triu(x) 输入矩阵x，输出v是x的上三角阵； v=tril(x) 输入矩阵x，输出v是x的下三角阵；,v=triu(x,1) 同上，但对角元素为0，可从A中提取U; v=tril(x,-1) 同上，但对角元素为0，可从A中提取L。,例. 用迭代法解,shiyan53,MATLAB对稀疏矩阵的处理: 进行大规模计算的优点,a=sparse(r,c,v,m,n) 在第r行、第c列输入数值v，矩阵共m行n列，输出a为稀疏矩阵，只给出(r,c)及v,aa=full(a) 输入稀疏矩阵a，输出aa为满矩阵（包含零元素）,a=spar
14、se(2,2:3,8,2,4), aa=full(a),n=500;b=1:n; a1=sparse(1:n,1:n,4,n,n); a2=sparse(2:n,1:n-1,1,n,n); a=a1+a2+a2; tic;x=ab;t1=toc aa=full(a); tic;xx=aab;t2=toc y=sum(x) yy=sum(xx),例. 分别用稀疏矩阵和满矩阵求解Ax=b, 比较计算时间,t1, t2相差巨大，说明用稀疏矩阵计算的优点 （y=yy 用于简单地验证两种方法结果的一致）,shiyan54,y=f(x),使点(xi,yi) 与曲线 y=f(x)的距离i尽量小,i=1,n,
15、曲线拟合与最小二乘准则,用MATLAB作线性最小二乘拟合,1. 作多项式 f(x)=a1xm+ +amx+am+1拟合,可利用已有程序:,a=polyfit(x,y,m),输入:数据x,y (同长度数组);m (拟合多项式次数) 输出:系数a=a1, am , am+1 (数组)。,2. 对超定方程组,仍用,可得最小二乘意义下的解,多项式在x点的值: y=polyval(a,x),例 汽车刹车距离,刹车距离d与车速v的关系：,数据拟合,(#),车速与刹车距离的实际数据记作(vi, di), i=1,2,7,结果,shiyan57,用MATLAB作多项式运算,1.多项式乘法：p=conv(p1,
16、p2),3.多项式生成: p=poly(x),2.多项式除法,q,r=deconv(p1,p2),生成方阵x的特征多项式，或以向量x为根,4.多项式求根: r=roots(p),8.部分分式展开: r,p,k=residue(b,a),9.矩阵多项式求值: PM=polyvalm(p,M),5.多项式显示: PS=ploy2str(p,s),一般“卷积”运算,6.多项式求导: dp=ployder(p),7.多项式积分: ip=ployint(p),例：多项式求根 （也可说明问题的“病态性”）,考虑如下的问题 f(x)=(x-1)(x-2).(x-20) 显然方程 f(x)=0 的解是 1 2 3 4 19 20 请问: 如下方程的解是什么？,p=p
