第六章MATLAB 数值计算Word格式.docx

1、例6-3求多项式和的乘积M（x）。a=1 5 0; %第一个多项式F（x）b=2 1; %第二个多项式G（x）c=conv（a,b）; %求两个多项式的乘积Mx=poly2str（c, %用常用的方式表示多项式的积c,Mx%endc = 2 11 5 0Mx = 2 x3 + 11 x2 + 5 x例6-4的除运算D（x）。c=1 5 0; %第一个多项式F（x）a=2 1; %第二个多项式G（x）q,r=deconv（c,a）; %求F（x）/ G（x）Dx=poly2str（q, q,r,Dxq = 0.5000 2.2500r = 0 0 -2.2500Dx = 0.5 x + 2.25

2、程序说明：1. 在运行结果中变量q是F（x）除以G（x）的商，而r则是除不尽的余数。2. 运行结果变量Dx表示的商没有加上余数。6-1-3多项式的求导Dp=polyder（p）p为向量表示的多项式。例6-5的一阶导数。我们容易知道以上两个方程的导数手工验算结果为：F （x） =2x+5和G （x）=2我们看MATLAB的计算结果。f=1 5 0;g=2 1;Df=polyder（f）; %求F（x）的导数Dg= polyder（g）; %求G （x）的导数Dfx=poly2str（Df, Dgx=poly2str（Dg, Df,Dg,Dfx,DgxDf = 2 5Dg = 2Dfx = 2 x

3、 + 5Dgx =6-1-4 多项式的求根A.r=roots（p）另外还有一种通过先求多项式的伴随矩阵，然后再求特征值的方法也可以求得多项式的根。例6-6的根。我们容易知道方程的根为：F（x）为：x1=0；x2=-5G（x）为： x2=-1/2rf=roots（f）; %求F（x）的根rg= roots （g）; %求G （x）的根rf, rg, rf = 0 -5rg = -0.5000从运行结果我们可以看到F（x）的根为0和-5，G（x）的根为-0.50，这和我们手工验算的结果完全一致。例6-7f=1 5 -3;g=2 1 0 1; %求F（x）的根 %求G （x）的根 -5.5414 0

4、.5414 -1.0000 0.2500 + 0.6614i 0.2500 - 0.6614i1. 我们可以看到例6-6和例6-7求得的根和我们手工计算的结果是一致的，其中例6-7多项式的根出现了虚根。2. 我们用求伴随矩阵的方法对例6-7再做一次求根运算： %第一个多项式F（x）cf=compan（f）; %求F（x）的伴随矩阵cg=compan（g）; %求G （x）的伴随矩阵cf, cg, cf = -5 3 1 0cg = -0.5000 0 -0.5000 1.0000 0 0 0 1.0000 0ccf=eig（cf） %求特征根ccg=eig（cg） %求特征根ccf = -5.

5、5414ccg = -1.0000 我们看到两种方法的结果是一致的。6-2 数据分析6-1-1 极值、均值、标准差和中位值的计算A. Pmax=max（X）B. Pmin=min（X）C. Pmean=mean（X）D. Pstd=std（X）E. Pmedian=median （X）1. max（X） 、min（X）、mean（X）、std（X）、median （X）分别是用来求数组或矩阵的最大值、最小值、均值、标准差、中位值。 注意std（X）的定义是，std（X,1）的定义是。2. 这里X可以是数组也可以是矩阵。如果是数组则对整个数组进行计算，如果是矩阵则是分别对矩阵的每个列向量进行运算

6、。例6-8对一组随机数组进行均值、方差和中位值的计算。x=randn（1,10）; %生成一个10元素的数组Pmax=max（x）; %计算数组x的最大值Pmin=min（x）; %计算数组x的最小值Pmean=mean（x）; %计算数组x的平均值Pstd=std（x）; %计算数组x的标准差Psqu= Pstd2; %计算数组x的方差Pmed=median （x）; %计算数组x的中位值answers= Pmax Pmin Pmean Pstd Psqu Pmed;x,answersx = Columns 1 through 6 -0.4326 -1.6656 0.1253 0.2877

7、-1.1465 1.1909 Columns 7 through 10 1.1892 -0.0376 0.3273 0.1746answers = 1.1909 -1.6656 0.0013 0.9034 0.8162 0.15001. 注意方差为标准差的平方。2. 需要注意的是MATLAB的变量定义是区分大小写的。例6-9对一个随机矩阵进行最大值、最小值、均值、标准差、中位值的计算。X=randn（6,3）; %生成一个63的随机矩阵Pmax=max（X）; %计算矩阵X的最大值Pmin=min（X）; %计算矩阵X的最小值Pmean=mean（X）; %计算矩阵X的平均值Pstd=std（

8、X）; %计算矩阵X的标准差Pmed=median （X）; %计算矩阵X的中位值X,Pmax,Pmin,Pmean,Pstd,Pmed;X = -0.1867 1.0668 0.7143 0.7258 0.0593 1.6236 -0.5883 -0.0956 -0.6918 2.1832 -0.8323 0.8580 -0.1364 0.2944 1.2540 0.1139 -1.3362 -1.5937Pmax = 2.1832 1.0668 1.6236Pmin = -0.5883 -1.3362 -1.5937Pmean = 0.3519 -0.1406 0.3607Pstd = 0

9、.9962 0.8482 1.24046-2-2 曲线拟合A. P=polyfit（X,Y,N）B. Yval=polyval（P,X）1. polyfit（X,Y,N）根据输入数据X和Y生成一个N阶的拟合多项式。2. polyval（P,X）根据数据X，用拟合多项式P生成拟合好的数据。3. 这里X和Y构成了一系列数据点的坐标。例6-10下列数据为一段时间的日气温平均实际数据，求该数据的拟合方程。设这组气温数据为：18 19 17 18 20 22 22 23 21 21 20 22 23 23 22d=1:15; %时间t=18 19 17 18 20 22 22 23 21 21 20 2

10、2 23 23 22; %和时间相对应的日平均气温p=polyfit（d,t,3）; %生成拟合多项式（向量形式），阶次为3px=poly2str（p, %生成拟合多项式（字符形式）pv=polyval（p,d）; %生成拟合值p,pxplot（d,t, :*r,d,pv,-k, LineWidth, 2）title（fontsize18bf数据拟合图像Color,r）xlabel （fontsize14rm时间d轴kylabel （fontsize14rm气温t轴gtext（fontsize16fontname宋体实际曲线fontsize16fontname宋体拟合曲线p = 0.0010

11、-0.0538 0.9644 16.4623px = 0. x3 - 0. x2 + 0.96436 x + 16.4623图形如6-1所示：1. d表示时间第一天、第二天，t表示每天的气温。2. 本例我们选取的拟合多项式的阶次是3。3. p为向量形式的拟合多项式，px为字符形式的拟合多项式。4. plot为绘图函数，它的两个参数分别为图像的横坐标和纵坐标。4个参数则分别为两条曲线的横坐标和纵坐标。5. 图中的箭头可在绘图窗口中打开insert菜单，选择Arrow然后用鼠标从起始位置拖动到终止位置即可。例6-11根据上例中的数据建立的天气温度，拟合方程求5天后的天气气温是多少？程序设计及运行结

12、果：x=20; T=0.*x3 - 0.*x2 + 0.96436*x + 16.4623T = 22.59951. 依照上例7-10天的时间序列d，5天后应该是20天，所以令x=20。2. T即为上例中我们求得的拟合公式，结果求得5天后的日平均气温为22.5995。6-3 数值积分和数值微分6-3-1 微分和积分的数学表达式对函数f（x）在区间x1,x2上求积分在数学上常常表示为对函数F（x）求导数运算在数学上常常表示为：6-3-2 函数数值积分命令函数：1. 连续被积函数quadl（f,a,b,t,trace） 自适应递推牛顿科西（NewtonCotes）法求积分2. 离散被积函数trap

13、z（X,Y） 梯形法求数值积分1. 1） 参数是被积函数表达式字符串或者函数文件名。2） 参数a、b定义函数积分的上下限。3） 参数t定义积分的精度（默认值1.0e-6）；trace设置是否用图形展示积分过程，1为展示、0不展示。2. trapz积分中给出Y相对于X的积分值。当Y是（m*n）矩阵时，积分对Y的列向量分别进行，得到一个（1*n）矩阵是Y的列向量对应于X的积分结果。例6-12设，用牛顿科西法求这个积分可以用手工算出来是2，我们看MATLAB的计算结果。s=quadl（sin（x）,0,pi）s = 2.0000例6-13，求y=quadl（sin（x）.2+x.3）./（x-2.*

14、cos（3.*x）,0,2）ans = 11.1943例6-14用离散数据表示，积分区间为分别用离散积分方法对其进行积分。程序设计及运行结果（如图6-3）：X=0:0.001:pi; %离散化积分变量Y=sin（X）; %生成被积函数的离散序列Z=trapz（X,Y）; %梯形法求积分Z Z =例6-152; %离散化积分变量Y=（sin（X）.2+X.3）./（X-2.*cos（3.*X）;Z 10.9607输入公式时一定要注意算符都是对数组进行运算的，所以是“.*”“.”“./“。6-3-3 数值微分D=diff（Y）Y为一组离散序列，D为与Y同维的离散矩阵，则微分运算是对矩阵的每个列向量

15、微分。例6-16设函数进行微分运算，并画出计算结果的图像。dx=0.001;x=-3:dx:3; %生成一个从-3到3间隔为0.001的序列x1=-3:3-dx;y=sin（x）;z=2.*x.3+3.*x.2+x-1;Y=diff（y）; %对y进行微分运算DY= diff（y）/dx; %对y进行求导运算Z=diff（z）; %对z进行微分运算DZ=diff（z）/dx; %对z进行求导运算subplot（2,2,1）; %产生两行一列的绘图区间的第一个绘图区间plot（Y,1） fontsize12fontnameTime New Roman Deltay（x） grid onsubpl

16、ot（2,2,2）; %产生两行一列的绘图区间的第二个绘图区间plot（Z ,fontsize12fontnameTime New Roman Delta z （x）subplot（2,2,3）; %产生两行一列的绘图区间的第三个绘图区间plot（x1,DY, 1） fontsize12fontnameTime New Romany（x）=cosxxlabel（fontsize12fontnameTime New Romanxylabel（fontsize12fontnameTime New Romandy/dxaxis （-3 3 -1 1）subplot（2,2,4）; %产生两行一列的绘

17、图区间的第四个绘图区间plot（x1,DZ ,fontsize12fontnameTime New Roman z（x）=6x2+6x+1fontsize12fontnameTime New Romandz/dxaxis （-3 3 -20 80）运行结果（如图6-2）：1. 输入公式时一定要注意算符都是对数组进行运算的，所以是“.*”“.”“./”。2. Y=diff（y）和Z=diff（z）和依此画出的图像，纵坐标是两函数的微分，横坐标是数据序列的值并不是实际图像的横坐标。由于纵坐标是，所以纵坐标在数值上是导数的倍。3. 导数图的纵坐标是，n的取值范围应该是因此取值是到，而微分相应取值是，

18、所以两序列的数目相差1，所以画图时不能用plot（x,diff（y） /dx）。6-4 一般非线性方程组的数值解x=fsolve（fun,x0,options）1. 是我们要求解的方程组，可以直接输入，不过由于方程组比较复杂，一般都先生成一个m文件，然后再进行调用。2. x0是给出的这个方程组的初值解，我们可以随意地给出。3. options是命令函数fsolve的参数设置项。因为fsolve求解的过程本身就是一个优化的过程，而options则是设置优化过程的参数的。这里对于一般的非线性方程组求解我们常设置为optiomset（Displayoff），意思为不显示优化过程。例：6-17求方程组

19、的数值解。这个方程组的解经过手工演算我们知道解为：x=4/3，y=-2/3，z=-2/3。程序设计与运行结果：1. 建立函数文件xyz，并保存：function fun=xyz（X）x=X（1）;y=X（2）;z=X（3）;fun=zeros（3,1）;fun（1）=x+y+z;fun（2）=x-y+3*z;fun（3）=3*x+y-z-4;2. 在MATLAB的工作区里求解：X0=1 1 1;X=fsolve（xyz,X0, optimset （） 1.3333 -0.6667 -0.66671. 在建立一个函数文件时一定要以function作为文件头说明。2. 把建立的函数文件保存到MAT

20、LAB默认的文件夹里，如work文件夹。3. zeros（3,1）用来产生一个3行1列的矩阵，矩阵的所有元素都为0。4. 我们可以看到计算的结果和我们手工计算的结果非常近似，差别是计算精度引起的，fsolve的默认求解精度是0.001。5. 事实上对方程组的求解过程是一个优化过程，optimset （）设置不显示求解过程。6-18求非线性方程组这个方程组手工计算是难以完成的。1. 建立函数文件fun，并保存：function Q=fun（T）x=T（1）;y=T（2）;z=T（3）;Q=zeros（3,1）;Q（1）=x+y-z;Q（2）=cos（x）+y3+z-5;Q（3）=4*x+2y+l

21、og（z）-1;X=fsolve（fun, 1 1 1, optimset （ -0.4392 1.4548 1.0156fsolve的默认求解精度是0.001，方程成立的精度为0.001。6-5 微分方程的数值解6-5-1 微分方程的基本形式的导数等于，即：6-5-2 一阶常微分方程的求解A. T,Y=ode45（odefun,tspan,y0）B. T,Y=ode23（odefun,tspan,y0）1. ode45（）和ode23（）是两个求解微分方程的命令函数，ode45（）应用最广，在容许误差大的情况下ode23（）的效率要比ode45（）高一些。2. T和Y分别表示微分方程的解的两个变量的数值序列。3. odefun是待解微分方程表达式的函数文件名。4. tspan表示运算的起至时间，是行向量，例如0 10表示运算时间从0开始到10结束。5. y0初始状态，用列向量表示，其元素分别表示微分方程组各个表达式的初始状态值。6-19求解微分方程并画出解在区间0 2*pi上的图像。1. 在程序编辑器中建立待解微分方程表达式的函数文件funx.m，并存入work工作目录：function Y=funx（x,y）Y=cos（x）;2. 在MATLAB的命令窗口求解，结果如图6-3所示：x,y=ode45（funx,0 2*pi,0）;