第六章MATLAB数值计算.doc

第六章 MATLAB 数值计算6-1 多项式的运算61-1 多项式的生成和表达1. 多项式的表达在MATLAB环境下多项式是用向量的形式表达的。 向量最右边的元素表示多项式的0阶,向左数依次表示多项式的第1阶、第2阶、第3阶。例如多项式表示为：5 0 3 2 1。2. 多项式的生成语法:P=ploy（MA)说明:1. 若MA为方阵，则生成的多项式P为方阵MA的特征多项式。2. 若MA为向量，则向量和多项式满足这样一种关系:，生成的多项式为:3. 直接输入的方式生成多项式。例6-1利用方阵M=5 6 7;8 9 1;11 12 13生成一个多项式（为方阵M的特征多项式）.程序设计：clearM=5 6 7；8 9 1;11 12 13；P=poly(M）; %产生多项式的向量表达式Px=poly2str（P, x）; 生成常见的多项式表示形式P,Px运行结果:P = 1.0000 27.0000 90。0000 54。0000Px = x3 - 27 x2 + 90 x + 54例6-2利用向量A=2 3 4 5生成一个多项式。程序设计：clearA=2 3 4 5；P=poly(A); 产生多项式的向量表达式Px=poly2str（P， x）; %生成常见的多项式表示形式P，Px运行结果：P = 1 14 71 154 120Px = x4 14 x3 + 71 x2 154 x + 12061-2 多项式的乘除语法:A. c=conv（a，b）B. q,r=decony（c，a）说明：1. a、b和c分别是多项式的向量表示形式。A表示两个多项式的乘积运算，B表示两个多项式的除法运算.2. q表示除运算的商，r表示除运算的余数.例63求多项式和的乘积M(x）。程序设计:cleara=1 5 0; 第一个多项式F(x）b=2 1; 第二个多项式G（x)c=conv(a，b）； 求两个多项式的乘积Mx=poly2str(c, x)； %用常用的方式表示多项式的积c，Mxend运行结果:c = 2 11 5 0Mx = 2 x3 + 11 x2 + 5 x例64求多项式和的除运算D（x)。程序设计：clearc=1 5 0； %第一个多项式F（x)a=2 1; 第二个多项式G(x）q，r=deconv（c,a）; %求F（x)/ G(x)Dx=poly2str(q, x); 用常用的方式表示多项式的积q，r,Dx%end运行结果：q = 0。5000 2.2500r = 0 0 2。2500Dx = 0.5 x + 2.25程序说明：1. 在运行结果中变量q是F(x)除以G（x)的商，而r则是除不尽的余数。2. 运行结果变量Dx表示的商没有加上余数.6-13多项式的求导语法：Dp=polyder(p）说明：p为向量表示的多项式.例6-5求多项式和的一阶导数.我们容易知道以上两个方程的导数手工验算结果为：F (x) =2x+5和G （x）=2我们看MATLAB的计算结果。程序设计：clearf=1 5 0； %第一个多项式F（x)g=2 1; 第二个多项式G（x）Df=polyder（f)； 求F（x）的导数Dg= polyder(g); 求G (x）的导数Dfx=poly2str(Df, x)；Dgx=poly2str(Dg, x);Df，Dg，Dfx，Dgx%end运行结果:Df = 2 5Dg = 2Dfx = 2 x + 5Dgx = 26-14 多项式的求根语法：A.r=roots(p)说明:另外还有一种通过先求多项式的伴随矩阵,然后再求特征值的方法也可以求得多项式的根。例6-6求多项式和的根。我们容易知道方程的根为：F（x）为:x1=0;x2=5G（x）为： x2=1/2我们看MATLAB的计算结果。程序设计：clearf=1 5 0; %第一个多项式F(x）g=2 1; %第二个多项式G(x）rf=roots（f)； %求F(x)的根rg= roots (g）； 求G （x）的根rf, rg， %end运行结果：rf = 0 5rg = 0.5000程序说明：从运行结果我们可以看到F（x）的根为0和5,G（x）的根为0.50，这和我们手工验算的结果完全一致。例67求多项式和的根。程序设计:clearf=1 5 -3; 第一个多项式F(x)g=2 1 0 1； %第二个多项式G（x)rf=roots（f); 求F（x)的根rg= roots (g)； 求G （x）的根rf， rg, end运行结果：rf = 5。5414 0。5414rg = -1。0000 0.2500 + 0.6614i 0.2500 - 0.6614i程序说明：1. 我们可以看到例66和例6-7求得的根和我们手工计算的结果是一致的,其中例67多项式的根出现了虚根。2. 我们用求伴随矩阵的方法对例67再做一次求根运算:clearf=1 5 -3； 第一个多项式F（x）g=2 1 0 1; %第二个多项式G（x)cf=compan（f); 求F(x）的伴随矩阵cg=compan(g）; 求G (x)的伴随矩阵cf， cg, %end运行结果：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.5414 0.5414ccg = 1。0000 0。2500 + 0.6614i 0.2500 - 0.6614i我们看到两种方法的结果是一致的.6-2 数据分析6-11 极值、均值、标准差和中位值的计算语法：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-8对一组随机数组进行均值、方差和中位值的计算。程序设计：clearx=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，answers运行结果：x = Columns 1 through 6 0.4326 -1.6656 0.1253 0.2877 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.1500程序说明:1. 注意方差为标准差的平方。2. 需要注意的是MATLAB的变量定义是区分大小写的。例6-9对一个随机矩阵进行最大值、最小值、均值、标准差、中位值的计算。程序设计：clearX=randn（6,3); 生成一个63的随机矩阵Pmax=max(X）； 计算矩阵X的最大值Pmin=min(X）； %计算矩阵X的最小值Pmean=mean(X)； 计算矩阵X的平均值Pstd=std(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。9962 0。8482 1。240462-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构成了一系列数据点的坐标。例610下列数据为一段时间的日气温平均实际数据,求该数据的拟合方程。设这组气温数据为：18 19 17 18 20 22 22 23 21 21 20 22 23 23 22程序设计:cleard=1:15； 时间t=18 19 17 18 20 22 22 23 21 21 20 22 23 23 22; %和时间相对应的日平均气温p=polyfit(d，t，3)； %生成拟合多项式（向量形式），阶次为3px=poly2str（p, x）; %生成拟合多项式（字符形式）pv=polyval(p,d）； %生成拟合值p,pxplot（d,t， ：r，d,pv,-k， LineWidth, 2）title(fontsize18bf数据拟合图像, Color,r)xlabel (fontsize14rm时间d轴, Color,k)ylabel （fontsize14rm气温t轴, Color,k）gtext（fontsize16fontname宋体实际曲线， Color,r）gtext(fontsize16fontname宋体拟合曲线， Color，k）end运行结果:p = 0.0010 0。0538 0.9644 16。4623px = 0。0010474 x3 - 0。053823 x2 + 0。96436 x + 16。4623图6-1 拟合曲线图图形如6-1所示：程序说明：1. d表示时间第一天、第二天,t表示每天的气温。2. 本例我们选取的拟合多项式的阶次是3。3. p为向量形式的拟合多项式,px为字符形式的拟合多项式。4. plot为绘图函数,它的两个参数分别为图像的横坐标和纵坐标。4个参数则分别为两条曲线的横坐标和纵坐标。5. 图中的箭头可在绘图窗口中打开insert菜单，选择Arrow然后用鼠标从起始位置拖动到终止位置即可。例6-11根据上例中的数据建立的天气温度，拟合方程求5天后的天气气温是多少？程序设计及运行结果:clearx=20; T=0。0010474x3 0。053823x2 + 0.96436*x + 16.4623T = 22.5995程序设计及运行结果：1. 依照上例710天的时间序列d，5天后应该是20天，所以令x=20。2. T即为上例中我们求得的拟合公式，结果求得5天后的日平均气温为22.5995。6-3 数值积分和数值微分63-1 微分和积分的数学表达式对函数f（x)在区间x1,x2上求积分在数学上常常表示为对函数F（x)求导数运算在数学上常常表示为：6-3-2 函数数值积分命令函数：1. 连续被积函数quadl(f，a,b，t，trace) 自适应递推牛顿科西（NewtonCotes)法求积分2. 离散被积函数trapz(X,Y) 梯形法求数值积分说明：1.1) 参数f是被积函数表达式字符串或者函数文件名。2) 参数a、b定义函数积分的上下限。3) 参数t定义积分的精度(默认值1。0e-6)；trace设置是否用图形展示积分过程,1为展示、0不展示.2.trapz积分中给出Y相对于X的积分值。当Y是（m*n)矩阵时，积分对Y的列向量分别进行，得到一个(1*n）矩阵是Y的列向量对应于X的积分结果.例612设，用牛顿科西法求.这个积分可以用手工算出来是2，我们看MATLAB的计算结果。程序设计及运行结果：clears=quadl（sin（x),0,pi)s = 2。0000例6-13设，求。程序设计及运行结果:cleary=quadl（(sin（x)。2+x.3）./(x2.cos（3。x），0,2)%endans = 11.1943例6-14用离散数据表示，积分区间为分别用离散积分方法对其进行积分。程序设计及运行结果（如图6-3）：clearX=0:0。001:pi; 离散化积分变量Y=sin(X)； 生成被积函数的离散序列Z=trapz（X,Y)； 梯形法求积分Z %endZ = 2。0000例615设，求。程序设计及运行结果：clearX=0:0。001:2； %离散化积分变量Y=(sin（X）.2+X。3）./(X2.cos(3.*X)）；Z=trapz(X,Y）；ZendZ = 10.9607程序说明:输入公式时一定要注意算符都是对数组进行运算的，所以是“.*“.”“。/“。63-3 数值微分命令函数：D=diff（Y）说明：Y为一组离散序列，D为与Y同维的离散矩阵，则微分运算是对矩阵的每个列向量微分。例616设函数和进行微分运算,并画出计算结果的图像.程序设计：cleardx=0。001；x=3：dx:3； 生成一个从-3到3间隔为0。001的序列x1=3：dx:3dx；y=sin(x)；z=2.*x.3+3。*x。2+x1；Y=diff(y); 对y进行微分运算DY= diff(y）/dx； 对y进行求导运算Z=diff（z)； %对z进行微分运算DZ=diff（z)/dx; %对z进行求导运算subplot（2,2，1）； 产生两行一列的绘图区间的第一个绘图区间plot(Y,k， LineWidth,1） title(fontsize12fontnameTime New Roman Deltay（x） , Color,k)grid onsubplot（2，2,2）; %产生两行一列的绘图区间的第二个绘图区间plot（Z ，-k， LineWidth，1） title（fontsize12fontnameTime New Roman Delta z （x), Color，k)grid onsubplot（2，2,3）; 产生两行一列的绘图区间的第三个绘图区间plot(x1，DY，-k， LineWidth， 1） title（fontsize12fontnameTime New Romany(x）=cosx, Color，k)xlabel（fontsize12fontnameTime New Romanx, Color,k)ylabel(fontsize12fontnameTime New Romandy/dx, Color,k）axis （-3 3 -1 1）grid onsubplot（2,2,4）; %产生两行一列的绘图区间的第四个绘图区间plot(x1,DZ ,k, LineWidth, 1) title（fontsize12fontnameTime New Roman z(x）=6x2+6x+1, Color,k）xlabel(fontsize12fontnameTime New Romanx， Color,k）ylabel（fontsize12fontnameTime New Romandz/dx, Color，k）axis （-3 3 20 80）grid onend图6-2 函数的微分与导数图运行结果（如图62)：程序说明：1. 输入公式时一定要注意算符都是对数组进行运算的,所以是“。*“.“./”。2. Y=diff（y）和Z=diff（z）和依此画出的图像，纵坐标是两函数的微分,横坐标是数据序列的值并不是实际图像的横坐标。由于纵坐标是,所以纵坐标在数值上是导数的倍。3. 导数图的纵坐标是，n的取值范围应该是。因此取值是到，而微分相应取值是到，所以两序列的数目相差1，所以画图时不能用plot（x，diff(y） /dx)。64 一般非线性方程组的数值解命令函数：x=fsolve（fun,x0，options）说明:1. fun是我们要求解的方程组，可以直接输入，不过由于方程组比较复杂,一般都先生成一个m文件，然后再进行调用。2. x0是给出的这个方程组的初值解,我们可以随意地给出。3. options是命令函数fsolve的参数设置项。因为fsolve求解的过程本身就是一个优化的过程,而options则是设置优化过程的参数的.这里对于一般的非线性方程组求解我们常设置为optiomset（Display, off），意思为不显示优化过程。例：617求方程组的数值解。这个方程组的解经过手工演算我们知道解为：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)=3x+y-z4;2. 在MATLAB的工作区里求解:X0=1 1 1;X=fsolve(xyz,X0, optimset （Display， off)X = 1。3333 -0。6667 0.6667程序说明：1. 在建立一个函数文件时一定要以function作为文件头说明.2. 把建立的函数文件保存到MATLAB默认的文件夹里，如work文件夹。3. zeros(3，1）用来产生一个3行1列的矩阵，矩阵的所有元素都为0.4. 我们可以看到计算的结果和我们手工计算的结果非常近似，差别是计算精度引起的，fsolve的默认求解精度是0.001。5. 事实上对方程组的求解过程是一个优化过程，optimset （Display， off)设置不显示求解过程。例：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)=4x+2y+log(z)1；2. 在MATLAB的工作区里求解:clearX=fsolve(fun, 1 1 1, optimset (Display, off）X = -0.4392 1.4548 1.0156程序说明:fsolve的默认求解精度是0.001，方程成立的精度为0.001。65 微分方程的数值解65-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所示：clearx，y=ode45（funx，0 2pi，0); %可以用 x，y=ode45(funx，0：0。01：2*pi,0)；plot(x，y，k, LineWidth, 2） title（fontsize18fontnameTime New Roman微分方程y（x)=cosx求解图像, Color，r）xlabel（fontsize16fontnameTime New Romanx， Color，k）ylabel（fontsize16fontnameTime New Romany=sinx， Color,k)axis （0 2pi 1 1)图6-3 一阶微分方程解的图像1grid on程序说明：1. 注意函数文件建立的格式,一定要以function作文件头，函数名为funx,两个变量为x和y。这里还要注意两个变量x和y的顺序，否则绘出的图像横坐标和纵坐标就颠倒了。2. 要把建立的函数文件保存在MATLAB默认的文件夹里,否则在调用时要指明路径。3. ode45（funx,0 2*pi，0),注意函数文件的调用格式。本例中计算时间（x变量)是从0开始,到2pi结束，初始值为0。4. 本例的微分方程是我们熟悉的，它的解和绘出的图像与手工计算结果是吻合的.例：6-20求解微分方程并画出解在区间0 1上的图像。程序设计与运行结果：1. 在程序编辑器中建立待解微分方程表达式的函数文件funxx.m，并存入work工作目录：function Y=funxx（x，y）Y=5*y+(x1）*sin（x）+（x+1)cos(x)；2. 在MATLAB的命令窗口求解，结果如图66所示：clearx,y=ode45（funxx,0 1,1)；plot(x,y，k， LineWidth, 2） title（fontsize18fontnameTime New Roman微分方程y(x）=5y+(x1)sinx+（x+1）cosx求解图像， Color，r）xlabel（fontsize16fontnameTime New Romanx， Color,k)ylabel(fontsize16fontnameTime New Romany, Color,k）axis (0 1 0 1。8e+2)图6-4 微分方程解的图像2grid on程序说明：1. 在本例中我们起的函数名为funxx，两个变量为x和y。2. 本例就比较复杂了，手工求解比较困难，但是我们看到用MATLAB很容易就给出了解的图像。653 二阶常微分方程的求解二阶常微分方程的求解使用的仍然是上节介绍的函数，只是在微分方程式的表达上做了一些变换，转换成一阶微分微分方程的表达形式，然后求解。具体的运用，我们在下面的例子中介绍.例：6-21求解二阶微分方程并画出解在区间0 5上的图像.程序设计与运行结果:1. 我们先对微分方程进行变换。设：，于是微分方程可以表达成为一个微分方程组：2. 在程序编辑器中建立待解微分方程表达式的函数文件funerx。m,并存入work工作目录:function Q=funerx(t，x） %x有两个元素x（1）和x(2）Q=x（2);4 cos（t）； Q有两个列向量3. 在MATLAB的命令窗口求解，结果如图6-7所示：cleart，x=ode45（funerx,0 2pi,4 0）; %t,x=ode45(funerx，0：0.01:2pi,4 0）；plot（t,x（:，1），k， t，x(:,2），:k，LineWidth， 2） title(fontsize18fontnameTime New Roman微分方程x+4cost=0求解图像, Color，k）xlabel（fontsize16fontnameTime New Romant， Color，k）ylabel（fontsize16fontnameTime New Romanx, Color，k)legend （x=4cos(t),v=-4sin(t)axis (0 2pi 4 4)图6-5 二阶微分方程解的图像1grid on程序说明：1. t,x=ode45（funerx，0 2*pi,4 0），指定计算区间（t的范围）为0 2pi，指定初始值（x的两个元素)分别为4和0，实际上相当于给定了初始位移是4,而初始速度为0。t，x=ode45(funerx,0:0.01：2pi,4 0）, 指定计算区间（t的范围）为0 2*pi，计算的步长为0。01。而在t，x=ode45（funerx，0 2*pi,4 0）中步长采用默认值。2. plot(t，x(:,1）,k, t,x（：,2），：k,LineWidth， 2），分别绘出t与x的第一列和第二列数据的图像，第一列第二列的图像分别用实线和点线绘出.3. 图67中实线线表示解函数x=4cost在区间0 2pi上的图像，而另一条点线则是x的导数，即v=-4sint在区间0 2*pi上的图像。例：6-22求解二阶微分方程并画出解在区间0 5上的图像。程序设计与运行结果：1. 我们先对微分方程进行变换.设：,于是微分方程可以表达成为一个微分方程组：2. 在程序编辑器中建立待解微分方程表达式的函数文件funerxx.m，并存入work工作目录：function Q=funerxx（t，x） %x有两个元素x（1）和x（2）Q=x（2)； x（1)- （1+ x(1)2)*x（2); %Q有两个列向量3. 在MATLAB的命令窗口求解，结果如图66所示：cleart,x=ode45（funerxx,0 12,1 1）；plot(t，x(：,1），k， t,x（：，2），:k,LineWidth, 2) title（fontsize18fontnameTime New Roman微分方程x+（x2+1)x+x=0求解图像， Color，k）xlabel(fontsize16fontnameTime New Romant, Color，k)ylabel（fontsize16fontnameTime New Romanx， Color,k)legend (x，x）axis （0 12 0.6 1.4)grid on图6-6 二阶微分方程解的图像2程序说明：1. t,x=ode45（funerxx,0 12，1 1）,指定计算区间(t的范围）为0 12