MATLAB程序设计与应用下new.docx

1、MATLAB程序设计与应用下newMATLAB程序设计与实例应用(下)晋中学院物理与电子工程学院第六章MATLAB 的数值计算第六章 MATLAB 数值计算6-1 多项式的运算6-1-1 多项式的生成和表达1. 多项式的表达在MATLAB环境下多项式是用向量的形式表达的. 向量最右边的元素表示多项式的0阶，向左数依次表示多项式的第1阶、第2阶、第3阶。例如多项式表示为：5 0 3 2 1。2. 多项式的生成语法：P=ploy(MA)说明：1. 若MA为方阵，则生成的多项式P为方阵MA的特征多项式。2. 若MA为向量，则向量和多项式满足这样一种关系：，生成的多项式为：3. 直接输入的方式生成多项

2、式。例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); %生

3、成常见的多项式表示形式P,Px运行结果：P = 1 -14 71 -154 120Px = x4 - 14 x3 + 71 x2 - 154 x + 1206-1-2 多项式的乘除语法：A. c=conv(a,b)B. q,r=decony(c,a)说明：1. a、b和c分别是多项式的向量表示形式。A表示两个多项式的乘积运算，B表示两个多项式的除法运算。2. q表示除运算的商，r表示除运算的余数。例6-3求多项式和的乘积M(x)。程序设计：cleara=1 5 0; %第一个多项式F(x)b=2 1; %第二个多项式G(x)c=conv(a,b); %求两个多项式的乘积Mx=poly2str(

4、c, x); %用常用的方式表示多项式的积c,Mx%end运行结果：c = 2 11 5 0Mx = 2 x3 + 11 x2 + 5 x例6-4求多项式和的除运算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则是除不尽的

5、余数。2. 运行结果变量Dx表示的商没有加上余数。6-1-3多项式的求导语法：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

6、 = 2 5Dg = 2Dfx = 2 x + 5Dgx = 26-1-4 多项式的求根语法：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.50

7、00程序说明：从运行结果我们可以看到F(x)的根为0和-5，G(x)的根为-0.50，这和我们手工验算的结果完全一致。例6-7求多项式和的根。程序设计：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. 我们可以看到例6-6和例6-7求得的根和我们手工计算的结果是一致的，其中例6-7多

8、项式的根出现了虚根。2. 我们用求伴随矩阵的方法对例6-7再做一次求根运算：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

9、 0.2500 - 0.6614i我们看到两种方法的结果是一致的。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可以是数组也可以是矩阵。如果是数组则对整个数组进行计算，如果是矩阵则是分别对矩阵的每个列向

10、量进行运算。例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

11、.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的最小

12、值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.5

13、937Pmean = 0.3519 -0.1406 0.3607Pstd = 0.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 22程

14、序设计：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 (fontsize14

15、rm气温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所示：程序说明：1. d表示时间第一天、第二天，t表示每天的气温。2. 本例我们选取的拟合多项式的阶次是3。3. p为向量形式的拟合多项式，px为字符形式的拟合多项式。4. plot为绘图函数，它的两个参数

16、分别为图像的横坐标和纵坐标。4个参数则分别为两条曲线的横坐标和纵坐标。5. 图中的箭头可在绘图窗口中打开insert菜单，选择Arrow然后用鼠标从起始位置拖动到终止位置即可。例6-11根据上例中的数据建立的天气温度，拟合方程求5天后的天气气温是多少？程序设计及运行结果：clearx=20; T=0.0010474*x3 - 0.053823*x2 + 0.96436*x + 16.4623T = 22.5995程序设计及运行结果：1. 依照上例7-10天的时间序列d，5天后应该是20天，所以令x=20。2. T即为上例中我们求得的拟合公式，结果求得5天后的日平均气温为22.5995。6-3

17、数值积分和数值微分6-3-1 微分和积分的数学表达式对函数f(x)在区间x1,x2上求积分在数学上常常表示为对函数F(x)求导数运算在数学上常常表示为：6-3-2 函数数值积分命令函数：1. 连续被积函数quadl(f,a,b,t) 自适应递推牛顿科西(NewtonCotes)法求积分2. 离散被积函数trapz(X,Y) 梯形法求数值积分说明：1. 1) 参数f是被积函数表达式字符串或者函数文件名。2) 参数a、b定义函数积分的上下限。3) 参数t定义积分的精度(默认值1.0e-6)。 2. trapz积分中给出Y相对于X的积分值。当Y是(m*n)矩阵时，积分对Y的列向量分别进行，得到一个(

18、1*n)矩阵是Y的列向量对应于X的积分结果。例6-12设，用牛顿科西法求。这个积分可以用手工算出来是2，我们看MATLAB的计算结果。程序设计及运行结果：clears=quadl(sin(x),0,pi)s = 2.0000例6-13设，求。程序设计及运行结果：cleary=quadl(sin(x).2+x.3)./(x-2.*cos(3.*x),0,2)%endans = 11.1943例6-14用离散数据表示，积分区间为分别用离散积分方法对其进行积分。程序设计及运行结果(如图6-3)：clearX=0:0.001:pi; %离散化积分变量Y=sin(X); %生成被积函数的离散序列Z=tr

19、apz(X,Y); %梯形法求积分Z %endZ = 2.0000例6-15设，求。程序设计及运行结果：clearX=0:0.001:2; %离散化积分变量Y=(sin(X).2+X.3)./(X-2.*cos(3.*X);Z=trapz(X,Y);Z%endZ = 10.9607程序说明：输入公式时一定要注意算符都是对数组进行运算的，所以是“.*”“.”“./“。6-3-3 数值微分命令函数：D=diff(Y)说明：Y为一组离散序列，D为与Y同维的离散矩阵，则微分运算是对矩阵的每个列向量微分。例6-16设函数和进行微分运算，并画出计算结果的图像。程序设计：cleardx=0.001;x=-3

20、:dx:3; %生成一个从-3到3间隔为0.001的序列x1=-3:dx: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,-k, LineWidth,1) title(fontsize12fontnameTime New Roman Deltay(x) , Color,k)grid onsubplot(2,2,2)

21、; %产生两行一列的绘图区间的第二个绘图区间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(fontsiz

22、e12fontnameTime 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

23、(-3 3 -20 80)grid on%end运行结果（如图6-2）：程序说明：1. 输入公式时一定要注意算符都是对数组进行运算的，所以是“.*”“.”“./”。2. Y=diff(y)和Z=diff(z)和依此画出的图像，纵坐标是两函数的微分，横坐标是数据序列的值并不是实际图像的横坐标。由于纵坐标是，所以纵坐标在数值上是导数的倍。3. 导数图的纵坐标是，n的取值范围应该是-1。因此取值是到，而微分相应取值是到，所以两序列的数目相差1，所以画图时不能用plot(x,diff(y) /dx)而是plot(x1,diff(y) /dx)。6-3 一般非线性方程组的数值解命令函数：x=fsolve

24、(fun,x0,options)说明：1. fun是我们要求解的方程组，可以直接输入，不过由于方程组比较复杂，一般都先生成一个m文件，然后再进行调用。2. x0是给出的这个方程组的初值解，我们可以随意地给出。3. options是命令函数fsolve的参数设置项。因为fsolve求解的过程本身就是一个优化的过程，而options则是设置优化过程的参数的。这里对于一般的非线性方程组求解我们常设置为optiomset(Display, off)，意思为不显示优化过程。例：6-17求方程组的数值解。这个方程组的解经过手工演算我们知道解为：x=4/3，y=-2/3，z=-2/3。程序设计与运行结果：1

25、. 建立函数文件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 (Display, off)X = 1.3333 -0.6667 -0.6667程序说明：1. 在建立一个函数文件时一定要以function作为文件头说明。2. 把建立的函数文件保存到MATLAB默认的文件夹里，如work文件夹。3. zeros(3,1)用来产

26、生一个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)=4*x+2y+log(z)-1

27、;2. 在MATLAB的工作区里求解：clearX=fsolve(fun, 1 1 1, optimset (Display, off)X = -0.4392 1.4548 1.0156程序说明：fsolve的默认求解精度是0.001，方程成立的精度为0.001。6-4 微分方程的数值解6-4-1 微分方程的基本形式的导数等于，即：6-4-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)