《matlab与数值分析精品资料.doc》由会员分享，可在线阅读，更多相关《matlab与数值分析精品资料.doc（41页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、实验指导书：Matlab 软件与数值计算Matlab 的含义是矩阵实验室(Matrix Laboratory),是由美国MathWork公司于1984年推出的一套高性能的数值计算软件,由于它具有使用方便,用户界面友好的特征,能集数值分析,矩阵计算,信号处理和图形显示于一身,更有多种用户工具包可选用,深受计算工作者和工程技术人员的亲睐,目前使用十分广泛.本书特选用Matlab作为数值分析的演示及实验的软件环境.本章只介绍与本书演示及实验有关的Matlab知识,更详细的Matlab知识可参考相关书籍.10.1 矩阵与数组1. 行向量 x=2 -3 1 %中间用空格将数据分开,也可用逗号分开x =2
2、 -3 1 %显示输出结果2. 列向量 y=3 -1 5 % 表示转置y = 3 -153. 矩阵输入A=1 3 2;5 4 6;7 9 8 % ;表示换行A = 1 3 2 5 4 6 7 9 84. 矩阵转置 B=AB = 1 5 7 3 4 92 6 85. 单位阵 I=eye(3) % 3 表示矩阵是3阶方阵 I = 1 0 0 0 1 00 0 16. 零矩阵 Z=zeros(3)Z = 0 0 0 0 0 0 0 0 07. 矩阵加减法 C=A+B,C = 2 8 9 8 8 15 9 15 16D=A-BD = 0 -2 -5 2 0 -35 3 08. 矩阵乘法 E=A*B,E
3、 = 14 29 50 29 77 11950 119 1949. 数组输入除上述的矩阵输入法外,数组还可以有如下输入法:(1) 等差输入 数组=初值:增量:终值a=1:2:10a =1 3 5 7 9增量为1时,可以省略 a=1:5a = 1 2 3 4 5(2) 等分输入法 数组=linspace(初值,终值,等分点数) a=linspace(0,10,5)a = 0 2.5000 5.0000 7.5000 10.000010. 点运算数组的点乘运算是对应分量相乘;矩阵的点乘是对应元素相乘，类似的运算还有点除（./）和点乘方(.).如: x=1:5x = 1 2 3 4 5 y=-2:2
4、y = -2 -1 0 1 2z=x.*yz = -2 -2 0 4 1010.2 函数运算和作图10.2.1 基本初等函数1. 三角函数: sin(x), cos(x), tan(x);2. 反三角函数: asin(x), acos(x), atan(x);3. 指数函数: exp(x) ; 4. 对数函数:(1) 自然对数:log(x); (2) 常用对数: log10(x); (3) 以2为底的对数: log2(x);5. 开平方函数: sqrt(x);6. 绝对值函数:abs(x)7. 计算一般函数值: eval(f )其中f是一个函数表达式的字符串 f=x.*sin(x); x=1:
5、10;y= eval(f)y = Columns 1 through 9 0.8415 1.8186 0.4234 -3.0272 -4.7946 -1.6765 4.5989 7.9149 3.7091 Column 10 -5.440210.2.2 多项式函数 多项式的表示 多项式可用一个向量表示，其分量为其降幂排列的系数如可表为p=1 3 0 5 多项式求值polyval(p,x) p=1 3 0 5; x=1:10; y=polyval(p,x)y =3 1 5 21 55 113 201 325 491 7053. 多项式的零点roots(p) roots(p)ans = 2.051
6、9 + 0.5652i 2.0519 - 0.5652i -1.1038 4. 方阵的特征多项式poly(A) A=3 1 2;1 3 2;1 2 3 p=poly(A)p =1.0000 -9.0000 20.0000 -12.0000 roots(ans)ans = 6.0000 2.0000 1.000010.2.3 矩阵函数1. 矩阵的行列式A=1 3 2;5 4 6;7 9 8 A = 1 3 2 5 4 6 7 9 8d=det(A),d =182. 矩阵的逆 B=inv(A),B =-1.2222 -0.3333 0.55560.1111 -0.3333 0.22220.9444
7、 0.6667 -0.61113. LU分解命令格式: L U P=lu(A)说明:这里L是单位下三角矩阵,U是上三角矩阵,P是置换阵,这是因为本命令中默认加上了按列选主元素.如命令中省略P, 则L不一定是单位下三角矩阵,而PL才是单位下三角矩阵.例: L U P=lu(A)L = 1.0000 0 0 0.7143 1.0000 0 0.1429 -0.7059 1.0000U = 7.0000 9.0000 8.0000 0 -2.4286 0.2857 0 0 1.0588P = 0 0 1 0 1 01 0 0 L U =lu(A)L = 0.1429 -0.7059 1.0000 0
8、.7143 1.0000 0 1.0000 0 0U = 7.0000 9.0000 8.0000 0 -2.4286 0.2857 0 0 1.05884. 矩阵范数命令格式: n1=norm(A,1), n2=norm(A,2), ninf=norm(A,inf)说明:这三个命令分别求A的1-范数,2-范数和-范数,其中n2=norm(A,2)也可写为n2=norm(A)例: A=3 1 2;1 3 2;2 2 3; n1=norm(A,1)n1 =7 n2=norm(A,2)n2 = 6.3723 ninf=norm(A,inf)ninf =75. 特征值分解命令格式: Q,D=eig(
9、A),D=eig(A)说明:此命令可以求出方阵A的特征值和相应的特征向量.例: A=3 1 2;1 3 2;2 2 3A = 3 1 2 1 3 2 2 2 3 D=eig(A)D = 0.6277 2.0000 6.3723 Q,D=eig(A)Q = 0.4544 0.7071 0.5418 0.4544 -0.7071 0.5418 -0.7662 0 0.6426D = 0.6277 0 0 0 2.0000 0 0 0 6.372310.2.4 绘图命令1. 二维绘图函数plot ,常用和最简单的绘图工具(1) 格式: plot(x,y,-r)(2) 说明:将x和y对应分量定义的点依
10、次用红色实线联结(x,y的维数必须一致).(3) 其中线型点型和色号可按下两表选择:表10.1线型符号点型符号实线(减号)实心圆点虚线（双减号）加号点线：星号*虚线间点（减号加点）空心圆点o（小写字母o）叉号x（小写字母x）表10.2兰色红色黄色绿色bryg(4) 绘图辅助命令利用以下命令可以为图形加上不同效果,注意以下命令必须放在相应的plot命令之后表10.3title()在图形上方加标题xlabel()为轴加说明ylabel()为轴加说明grid在图上显示虚线格标度text(x,y,)在图中x,y坐标处显示中的内容gtext( )在图中用光标确定位置处显示中的内容axis(xl,xu,y
11、l,yu)用中四个实数定义x,y方向显示范围hold on后面的plot图形将叠加在一起hold off解除hold on命令，后面plot产生的图形将冲去原有图形例10.1 作出-p,p间函数的图形x=-pi:0.05:pi;y=x.*x.*sin(2*x);plot(x,y,-)grid2. 一元函数绘制函数命令fplot ,可绘制已定义函数在指定区间上的图象,与plot命令类似,但能自适应地选择取值点.plot的线型和色号选项依然适用.格式: fplot(f,a,b),其中f是一个函数表达式的字符串,a,b是x取值的下界和上界.例10.2 作出-p,p间函数的图形f=x.*x.*sin(
12、2*x)fplot(f,-pi,pi)也可以作出与例10.1一样的图形.图1013. 多窗口绘图函数subplot格式:subplot(p,q,r)说明:将图形窗口分为p行q列共pq个格子,在第r个格子上画图,格子是从上到下,从左到右依次计数的.例10.3 作出-1,1间前四个Legendre多项式的图形Legendre多项式是定义在-1,1上的正交多项式系列,它可由下列递推公式生成:x=-1:0.05:1p0=1+0*x;p1=x;p2=(3*x.2-1)/2;p3=(5*x.3-3*x)/2subplot(2,2,1);plot(x,p0)subplot(2,2,2);plot(x,p1)
13、subplot(2,2,3);plot(x,p2)subplot(2,2,4);plot(x,p3)图10210.2.5 Matlab编程将一个完整的命令集合写成M文件便是一段Matlab程序. Matlab程序具有一般程序语言的基本结构和功能.有顺序结构,分支结构和循环结构三种结构.1. 顺序结构: 如程序中没有控制语句,则程序执行时按语句顺序逐句执行.2. 分支结构(条件语句):(1) if end(2) if elseend3. 循环结构 Matlab提供两种循环结构:(1) for 循环for =:end(2) while 循环while end4. 以上两种控制结构中的由逻辑表达式的
14、值决定,1为真,表示条件成立,0为假,表示条件不成立.逻辑表达式由关系表达式和逻辑运算符组合而成, 关系表达式由变量,常量,表达式和关系运算符组合而成,5. 关系表达式中的关系运算符有表10.4相等大于小于等于小于大于等于不等于6. 逻辑表达式中的逻辑运算符有表10.5逻辑与ABA与B同时成立逻辑或ABA或B成立逻辑非A与条件A相反的条件成立例10.5 找出1到100间的完全平方数for i=1:100 if sqrt(i)=fix(sqrt(i) i endend例10.6 找出2000年到2100年间的闰年.i=2000;while i0) %年份能整除4且逢100年能整除400 i en
15、d i=i+1;end7. 自定义函数如M文件的第一行为:function =(自变量表),则这个文件是一个自定义函数,可以如标准函数一样调用.例10.7 求方阵A的Jacobi迭代法的迭代矩阵Bfunction B=bj(A)n=max(size(A);D=zeros(n);for i=1:n D(i,i)=A(i,i);endBJ=-inv(D)*(A-D);以上程序存为名为bj.m的M文件,调用如下 A=3 1 2;1 3 2;1 2 3A = 3 1 2 1 3 2 1 2 3 BJ=bj(A)BJ= 0 -0.3333 -0.6667 -0.3333 0 -0.6667 -0.333
16、3 -0.6667 010.3 线性方程组的数值解10.3.1 直接法,A是一个nn方阵,b是n维列向量1. Gauss消元法 x=Ab 用带选主元素的Gauss 消元法求解.2. LU分解和Doolittle方法L U=lu(A),得到A的LU分解, Doolittle方法: L U=lu(A)y=Lbx=Uy10.3.2 迭代法L是A的下三角部分,U是上三角部分,D是对角部分.1. Jacobi迭代法原理:编程公式对例10.9 用Jacobi迭代法求解当退出运算解 编制程序,存入M文件jacobi.m% 输入数据A=10 -2 -1;-2 10 -1;-1 -2 5;b=3 15 10;e
17、=0.000001; %控制误差%n=max(size(A); %测定维数for i=1:n if A(i,i)=0 对角元为零,不能求解 return endendx=zeros(n,1) %设置初始解k=0; % 预设迭代次数为0kend=50 % 最大迭代次数为50r=1; % 前后项之差的无穷范数,初始值设为1while ke %达到预定精度或迭代超过50次退出计算 x0=x;% 记下前次近似解 for i=1:n s=0; for j=1:i-1 s=s+A(i,j)*x0(j); end for j=i+1:n s=s+A(i,j)*x0(j); end x(i)=b(i)/A(i
18、,i)-s/A(i,i); end r=norm(x-x0,inf);%重新计算前后项之差的无穷范数 k=k+1;endif kkend 迭代不收敛,失败else 求解成功 x kend此程序是一个通用程序，更改A,b,e后可求其他线性方程组的解当迭代次数超过50次(也可以更改)，还不能达到精度要求时，认为迭代失败，输出迭代失败信息迭代收敛时，给出符合精度要求的近似解和迭代次数如出现对角元为零时自动停止计算，并给出失败信息运行结果jacobix = 0.99999984084662 1.99999984084006 2.99999973808525k =162. Gauss-Seidel迭代法
19、原理:编程公式对编程时只要在jacobi.m中作一处更改便可,请读者作为习题完成,将相应文件存为g_s.m.运行结果g_s求解成功x = 0.99999989453805 1.99999994851584 2.99999995831395k = 9可见对此题Gauss-Seidel法比Jacobi法迭代收敛快.10.3.3 迭代法收敛理论L是A的下三角部分,U是上三角部分,D是对角部分.若能将转化为等价的，有，则此迭代格式收敛时，迭代不收敛越接近于零，收敛越快对Jacobi迭代法对Gauss-Seidel迭代法只要算出对应的则很容易判断迭代收敛或比较收敛的快慢.例10.11判别用Jacobi迭
20、代法和Gauss-Seidel迭代法的收敛性1. 前面已有函数bj.m，继续编制函数bs.m如下：function B=bs(A)n=max(size(A);D=zeros(n);for i=1:n for j=1:n if jrho=max(abs(eig(B)9. 判断A1 A1=10 -2 -1;-2 10 -1;-1 -2 5A1 = 10 -2 -1 -2 10 -1 -1 -2 5 BJ=bj(A1)BJ = 0 0.2000 0.1000 0.2000 0 0.1000 0.2000 0.4000 0 rho_j=max(abs(eig(BJ)rho_j = 0.3646 BS=
21、bs(A1)BS = 0 0.2000 0.1000 0 0.0400 0.1200 0 0.0560 0.0680BS = 0 0.2000 0.1000 0 0.0400 0.1200 0 0.0560 0.0680 rho_s=max(abs(eig(BS)rho_s = 0.1372可见两种迭代法均收敛，且Gauss-Seidel迭代法比Jacobi迭代法快2. 判断A2 A2=4 2 1;2 2 1;1 1 1A2 = 4 2 1 2 2 1 1 1 1 BJ=bj(A2)BJ = 0 -0.5000 -0.2500 -1.0000 0 -0.5000 -1.0000 -1.0000
22、 0 rho_j=max(abs(eig(BJ)rho_j = 1.2808 BS=bs(A2)BS = 0 -0.5000 -0.2500 0 0.5000 -0.2500 0 0 0.5000 rho_s=max(abs(eig(BS)rho_s = 0.5000可见Jacobi迭代法不收敛，Gauss-Seidel迭代法收敛10.5 方程和方程组求根10.5.1 二分法例10.17 求 在0,1间的根，要求编制程序equation_B如下：fun=2*x-exp(-x)%输入f(x)a=0;b=1;%输入求根区间a,be1=10(-7);%要求误差1e2=0.0001;%要求误差2r=1
23、;%当前误差,预设为1k=1;%迭代次数m=20;%最大迭代次数,预设为20x=a;fa=eval(fun);x=b;fb=eval(fun);if fa*fb0 不能求出实根 returnendwhile r=e2 & km x=a;fa=eval(fun); x=(a+b)/2; fx=eval(fun); if abs(fx)0 a=x; else b=x; end end r=b-a; k=k+1;endkfroot=x运行结果k =15froot = 0.351710.5.2 Newton法例10.18 求 在0.5附近的根，要求编制程序equation_N如下：fun=2*x-ex
24、p(-x)；%输入f(x)dfun=2+exp(-x)；%输入f (x)x0=0.5；%输入初始解e=10(-7);%要求误差r=1;%当前误差,预设为1k=1;%迭代次数m=20;%最大迭代次数,预设为10while r=e & k p=1 -2 3 -4；x= roots(p)ans = 1.6506 0.1747 + 1.5469i 0.1747 - 1.5469i2. 求函数的零点命令格式：x=fzero(fun,x0) 找出函数在x0 附近的零点. 这个命令是基于二分法的，它要求函数在x0 附近变号，否则可能失败并给出NaN结果例：求在1附近的根 fun=exp(-x/2)-x x=
25、fzero(fun,1)x = 0.703510.6 插值方法10.6.1 Lagrange插值在互异节点有函数值，求作n次Lagrange插值多项式并求在x_star处的插值结果例10.19 已知，用Lagrange插值求的近似值编制程序存入int_lagr.m,并运行x=4 9 16;%输入x值y=sqrt(x);%输入y值x_star=7;%输入x_star值n=max(size(x);%测定x的维数y_star=0;for i=0:n-1 lj=1; for j=0:n-1 if j=i lj=lj*(x_star-x(j+1)/(x(i+1)-x(j+1); end end y_st
26、ar=y_star+y(i+1)*lj;endy_star运行结果：y_star = 2.628610.6.2 Newton插值在互异节点有函数值，求作1到n次Lagrange插值多项式并求在x_star处的插值结果例10.20 已知的值，用Newton插值求的近似值编制程序存入int_newton.m,并运行x=pi/6 pi/4 pi/3 pi/2;%输入x值y=sin(x);%输入y值x_star=2*pi/9;%输入x_star值n=max(size(x);%测定x的维数y_star=y(1)xf=1;%一次因子的乘积预置为1dx=y;%各阶差商,预置为y值for i=1:n-1 %计
27、算各阶差商 dx0=dx;for j=1:n-i dx(j)=(dx0(j+1)-dx0(j)/(x(i+j)-x(j);end df=dx(1); % xf=xf*(x_star-x(i);%计算一次因子的乘积 y_star=y_star+xf*df%计算各次Newton插值的值end运行结果：y_star = 0.5000y_star = 0.6381y_star = 0.6434y_star = 0.642910.6.3 用拟合函数polyfit作插值当拟合多项式次数只比数据点个数少1时，它就是插值多项式所以可以用polyfit( )命令作插值它使用简单，可直接得到插值多项式的降幂标准形
28、式，插值结果与Lagrange插值和Newton插值完全一样命令格式：p=polyfit(x,y,n)这里 x,y是n+1 维数据向量，p得到的是一个n+1维向量，它的分量就是插值多项式按降幂排列的系数例：已知的值，用polyfit插值求的近似值 x=pi/6 pi/4 pi/3 pi/2; y=sin(x); p=polyfit(x,y,3)p = -0.0913 -0.1365 1.0886 -0.0195 polyval(p,2*pi/9)ans = 0.642910.7 数据拟合与函数逼近10.7.1 多项式数据拟合1. 多项式拟合命令polyfit( )命令格式:p=polyfit(
29、x,y,n) 求出n次最小二乘多项式的系数p ,这里x,y是数据向量,它们的维数应大于 n,当它们的维数是n+1时,这个命令求出的是插值多项式例：已知数据：表10.6xi-1.00-0.75-0.5-0.2500.250.50.751.00yi-0.22090.32950.88261.43922.00032.56453.13343.70614.2836求一次,二次,三次拟合多项式x=-1:0.25:1;y=-0.22090.32950.88261.43922.00032.56453.13343.70614.2836;p1=polyfit(x,y,1)p2=polyfit(x,y,2)p3=po
30、lyfit(x,y,3)运行结果:p1 = 2.2516 2.0131p2 = 0.0313 2.2516 2.0001p3 = 0.0021 0.0313 2.2501 2.00012. 构造正规方程求解多项式拟合polyfit命令并不使用正规方程，而是使用更复杂的奇异值分解方法，它可以避免正规方程组的病态如果想构造正规方程组求解可进行如下编程.例10.22 就上面的表10.6数据求次多项式拟合：x=-1:0.25:1;y=-0.22090.32950.88261.43922.00032.56453.13343.70614.2836;x0=x.0x1=xx2=x.2A=x0;x1;x2N=A
31、*Ab=A*yp=Nb运行结果：N = 9.0000 0 3.7500 0 3.7500 0 3.7500 0 2.7656b = 18.1183 8.4437 7.5870p = 2.0001 2.25160.0313注意这里p是升幂排列，与上次的结果是一致的10.7.2非线性拟合1. 非线性拟合如能事先化为线性拟合，最好事先转化用线性拟合处理，这样计算量小且精度高例10.23 已知数据：表10.7xi12345678yi15.320.527.436.649.165.687.8117.6求形如的经验公式(a,b为待定常数)解 编程如下:x=1:8;y=15.320.5 27.4 36.6 4
32、9.165.6 87.8 117.6;Y=log(y);p=polyfit(x,Y,1)p = 0.2912 2.4369a=exp(p(2)b=p(1)运行结果:a = 11.4371b= 0.2912有还可以画出拟合效果图xi=1:0.1:8；yi=a*exp(b*xi)；plot(x,y,+,xi,yi,-)图10810.8 数值积分10.8.1 非复化的数值积分实际上是一个线性组合计算，用Matlab编程计算很容易例10.25 分别用Cotes公式和三点GaussLegendre公式计算（准确值为1.7416）%输入被积函数和积分区间fun=x.*sin(x)a=0;b=2;%四阶Ne
33、wtonCotes公式A=7 32 12 32 7/90;n=max(size(A)-1;h=(b-a)/n;x=a:h:b;IC=eval(fun)*A*(b-a)%三点GaussLegendre求积公式x=1-0.7746,1,1.7746A=5 8 5/9IL=eval(fun)*A计算结果：IC = 1.7417IL = 1.7414其他积分的计算只要变更参数即可10.8.2 复化数值积分计算1 分n段的复化梯形公式这里n是等分段数，是步长，节点为例10.26 用分10段的复化梯形公式计算% 复化梯形积分公式fun=inline(x.*sin(x) % 输入函数a=0;b=2;% 输入
34、区间n=10;% 输入等分段数h=(b-a)/n;% 计算步长% 以下计算积分值I=feval(fun,a)+feval(fun,b);for k=2:n I=I+2*feval(fun,a+(k-1)*h);endI=I*h/2运行结果：I = 1.74192 逐次分半的复化梯形公式这里n是等分段数，是步长，节点为例10.27 用逐次分半的复化梯形公式计算,精确到.%逐次分半的复化梯形积分公式fun=inline(x.*sin(x);% 输入函数a=0;b=2;% 输入区间h=b-a; % 步长预设为区间长e=1e-4; % 要求精度% 以下计算积分值T=(feval(fun,a)+feva
35、l(fun,b)/2*hn=1;r=1; % 当前误差预设为1while re % 当前误差小于要求精度时退出计算 h=h/2;n=2*n; s=0; for k=2:2:n s=s+feval(fun,a+(k-1)*h); end T0=T; T=T0/2+s*h % 逐次分半公式 r=abs(T0-T); % 计算当前误差endT运行结果：T =1.74163 复化Simpson公式这里m是等分段数，是步长，节点为例10.28 用分10段的复化Simpson公式计算,精确到.%quad_4 复化Simpson积分公式fun=inline(x.*sin(x); % 输入函数a=0;b=2;% 输入区间m=10;% 输入等分段数h=(b-a)/m/2;% 计算步长% 以下计算积分值S=feval(fun,a)+feval(fun,b);for k=2:m S=S+2*feval(fun,a+2*(k-1)*h);endfor k=1:m