MATLAB应用基础第四章Matlab的数值计算

1、第4章MATLAB的数值计算4.1求解线性代数方程组已知线性代数方程组Ax = b可用下面三种方法求解1）直接求逆法求解x = inv(A)*b此法当A不可逆时失效2）左除法求解x = Ab左除法的基础是高斯消元法，由消元法对系数矩阵A进行LU分解。进而得到方程组的解。该法运算量少，运算速度快，而且数值稳定性好，解的精度高。当方程个数大于未知变量个数时，该法可得到问题的最小二乘解。当方程个数小于未知变量个数时，该法可求得有多个0元素的解。3）使用伪逆函数求最小范数解：x = pinv(A)b例如：%3个未知量4个方程A=1 2 3;4 5 6;7 8 0;2 5 8 b=366,804,351

2、,514%计算最小二乘解x=Ab %该馀向量具有最小范数res=A*x-b %生成4个未知量的3个方程A=Ab=b(1:3)%具有最多0元素x=Ab %计算最小范数解xn=pinv(A)*b 运行结果如下：A = 1 2 3 4 5 6 7 8 0 2 5 8b = 366 804 351 514x = 247.9818 -173.1091 114.9273res = -119.4545 11.9455 0 35.8364A = 1 4 7 2 2 5 8 5 3 6 0 8b = 366 804 351x = 0 -165.9000 99.0000 168.3000xn = 30.8182

3、-168.9818 99.0000 159.0545 4.2多项式计算1、多项式的表示方法多项式使用按降幂排列的多项式系数所构成的行向量描述。即C1Xn+C2Xn1+CnX+Cn+1可描述为：C1,C2,Cn,Cn+12、多项式相乘conv(a,b) 其中：a,b为多项式系数向量3、多项式相除q,r = deconv(a,b)其中：q为商，r为余项4、多项式相加两个相同规模的多项式可以相加：a + b注意：若规模不同，则要把低阶多项式补足若干个0。5、多项式微分1） 求多项式p的导数polyder(p) 2） 求多项式之积的导数polyder(a,b) 3）求多项式之比的导数q,d = pol

4、yder(a,b) 6、多项式求值polyval(p,x) 功能：求多项式p对应于x的值7、多项式求根roots(p)8、求一组根所对应的多项式若已知一多项式的根为行向量r时，则poly(r) 为r所对应的多项式。例如：a=1 2 3 4;b=1 4 9 16;c=conv(a,b)d=a+be=c+0 0 0 dq,r=deconv(e,a)h=polyder(c)g=polyder(a,b)q,d=polyder(a,b)p=1,-12,0,25,116;v=polyval(p,2.5)r=roots(p) 运行结果如下：c = 1 6 20 50 75 84 64d = 2 6 12 2

5、0e = 1 6 20 52 81 96 84q = 1 4 9 18r = 0 0 0 0 2 6 12h = 6 30 80 150 150 84g = 6 30 80 150 150 84q = 2 12 42 32 12d = 1 8 34 104 209 288 256v = 30.0625r = 11.7473 2.7028 -1.2251 + 1.4672i -1.2251 - 1.4672i 4.3数值逼近1、曲线拟合若已知离散数据向量x,y，则polyfit(x,y,N)将采用最小二乘法构造一个N阶多项式例如：x=0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1;

6、y=-.447,1.978,3.28,6.16,7.08,7.34, 7.66,9.56,9.48,9.3,11.2;p=polyfit(x,y,2)x1=0:.01:1;z=polyval(p,x1);p10=polyfit(x,y,10)z10=polyval(p10,x1);plot(x,y,bo-,x1,z,k:,x1,z10,r)legend(原始数据,2阶曲线拟和,10阶曲线拟和) 运行结果如下：p = -9.8108 20.1293 -0.0317p10 = 1.0e+006 * Columns 1 through 7 -0.4644 2.2965 -4.8773 5.8233

7、-4.2948 2.0211 -0.6032 Columns 8 through 11 0.1090 -0.0106 0.0004 -0.0000 2、一元插值若已知离散数据向量x,y，则yi = interp1(x,y,xi,method)将计算在xi处的函数值yi。其中的参数method可为：nearest 最邻近插值linear 线性插值（缺省）spline 三次样条插值cubic 三次多项式插值（要求x等距）xi为向量时，yi也为向量。例如：x=0:8;y=sin(x);x1=0:.1:8;y1=interp1(x,y,x1,nearest);y2=interp1(x,y,x1,lin

8、ear);y3=interp1(x,y,x1,spline);y4=interp1(x,y,x1,cubic);plot(x,y,o-,x1,y1,b-, x1,y2,k:,x1,y3,r-,x1,y4,g-)legend(原始数据,nearest, linear,spline,cubic) 运行结果如下： 3、二元插值若已知离散数据x,y,z，则zi=interp2(x,y,z,xi,yi,method)功能：计算在xi,yi处的zi值。其中:method可为：linear 线性插值(x,y单增)bilinear 双线性插值(x,y单增)cubic 三次插值(要求x,y单 增等距)bicub

9、ic 双三次插值(要求x,y 单增等距)nearest 最近邻域插值(要求x, y单增)zi=griddata(x,y,z,xi,yi)功能：griddata使用距离倒数的方法,计算非等距节点x,y,z在xi,yi处的zi值。若xi,yi为向量，则zi将为矩阵。(xi为行向量，yi为列向量)。xi和yi也可为由meshgrid函数生成矩阵。例如：%在500米间距正方形网格系统上测得的海底深度(米)%绘制海底深度图x=0:0.5:4; %0.5km为单位y=0:0.5:6; %0.5km为单位z=100 99 100 99 100 99 99 99 100 100 99 99 99 100 99

10、 100 99 99 99 99 98 98 100 99 100 100 100 100 98 97 97 99 100 100 100 99 101 100 98 98 100 102 103 100 100 102 103 101 100 102 106 104 101 100 99 102 100 100 103 108 106 101 99 97 99 100 100 102 105 103 101 100 100 102 103 101 102 103 102 100 99 100 102 103 102 101 101 100 99 99 100 100 101 101 100

11、100 100 99 99 100 100 100 100 100 99 99 99 99 100 100 100 99 99 100 99 100 99;mesh(x,y,z)xlabel(X-axis(km)ylabel(Y-axis(km)zlabel(海底深度(m)title(海底深度图)pausexi=linspace(0,4,30);yi=linspace(0,6,40);xxi,yyi=meshgrid(xi,yi);zzi=interp2(x,y,z,xxi,yyi,linear);mesh(xxi,yyi,zzi)title(海底深度图(linear)hold onxx,yy

12、=meshgrid(x,y);plot3(xx,yy,z+0.1,ob)hold offpausezzi=interp2(x,y,z,xxi,yyi,bilinear);mesh(xxi,yyi,zzi)title(海底深度图(bilinear)hold onxx,yy=meshgrid(x,y);plot3(xx,yy,z+0.1,ob)hold offpausezzi=interp2(x,y,z,xxi,yyi,cubic);mesh(xxi,yyi,zzi)title(海底深度图(cubic)hold onxx,yy=meshgrid(x,y);plot3(xx,yy,z+0.1,ob)

13、hold offpausezzi=interp2(x,y,z,xxi,yyi,bicubic);mesh(xxi,yyi,zzi)title(海底深度图(bicubic)hold onxx,yy=meshgrid(x,y);plot3(xx,yy,z+0.1,ob)hold offpausezzi=interp2(x,y,z,xxi,yyi,nearest);mesh(xxi,yyi,zzi)title(海底深度图(nearest)hold onxx,yy=meshgrid(x,y);plot3(xx,yy,z+0.1,ob)hold offpausezzi=interp2(x,y,z,xxi

14、,yyi,cubic);pcolor(xxi,yyi,zzi)shading interphold oncontour(xxi,yyi,zzi,15,k)colormap(cool)colorbar(vert)hold offtitle(海底深度等值线图)运行结果如下：4.4常微分方程的数值解由于任一高阶常微分方程均可化为与之等价的一阶常微分方程组(即状态方程)，所以MATLAB中总是求解一阶常微分方程组，而且需要用m文件定义方程组。例如: 若常微分方程为：设,则原方程可表为故可定义xp.m文件描述为：function xdot = xp(t,y)global MU %0MU=2 method

15、=method-1; l = 0;endelsemethod=-1;endwhile (method 6) method = ;while length(method) method = input(Select a method number: );endendif (method = 0) hold offreturnendOPTIONS=0;if method=2, OPTIONS(6)=1;elseif method=3, OPTIONS(6)=2;elseif method=5, OPTIONS(5)=1;endif test_longl = l + 1;elsel = ;endif

16、 method=4disp()disp( Choose any of the following line search methods)disp()disp( 1) Mixed Polynomial Interpolation)disp( 2) Cubic Interpolation)disp()while length(l) l = input(Select a line search number: );endif l=2, OPTIONS(7)=1; endendif method=4OPTIONS(14)=200;elseOPTIONS(14)=300;endx=-1.9,2;dis

17、p( )if method 1e-8, error(Optimization Toolbox in datdemo), endenddisp(sprintf(nValue of the function at the solution: %g, options(8) ); disp(sprintf(Number of function evaluations: %g, options(10) ); disp(sprintf(Number of gradient evaluations: %gnn, options(11) ); disp(Strike any key for menu)paus

18、eclfhold off contour(xx,yy,meshd,conts,w-) xlabel(x1)ylabel(x2)title(Minimization of the Banana function)hold onplot(-1.9,2,o)text(-1.9,2,Start Point)plot(1,1,o)text(1,1,Solution)end运行结果如下： Choose any of the following methods to minimize the banana function UNCONSTRAINED: 1) Broyden-Fletcher-Golfarb

19、-Shanno 2) Davidon-Fletcher-Powell 3) Steepest Descent 4) Simplex Search LEAST SQUARES: 5) Gauss-Newton 6) Levenberg-Marquardt 0) QuitSelect a method number: 04.8 求解线性规划问题MATLAB中求解的线性规划问题一般形式为： min fTx S.t Axb 或 min fTx S.t Axb x1xx2其中：x为n维未知向量fT=f1,f2,fn为目标函数系数A为约束系数矩阵b为列向量x1,x2为n维常数向量1、求解形式的线性规划问题

20、x,lambda,how = lp(f,A,b)其中: 返回x为最优解，lambda为相应最优拉格朗日因子，how为描述求解状况的字符串。2、求解形式的线性规划问题x = lp(f,A,b,x1,x2)或x = lp(f,A,b,x1,x2,x0) x0为初始点或x = lp(f,A,b,x1,x2,xo,n) n表示由A与b所给出的约束条件中前n个为等式例如:求解线性规划问题 min 2x1+2x2 s.t x1-x2=1 -x1+2x2=0f=2;2;A=-1,1;-1,2;b=-1;0;x=lp(f,A,b,0,0,)运行结果如下:x = 1 0又如:求解线性规划问题 min -3x1+

21、x2 s.t x1+2x2=4 -x1+x2=0f=-3;1;A=1,2;-1,1;b=4;2;x,la,how=lp(f,A,b,0,0,) 运行结果如下:x = 4 0la = 3.0000 0 0 7.0000how =ok 4.9求解二次规划MATLAB求解二次规划问题的一般形式为：例如：求解二次规划问题 s.t 显然：，关键是求H，因为式中的二次项：故求解二次规划的MATLAB函数为x = qp(H,f,A,b)或x = qp (H,f,A,b,x1,x2)其中x1与x2为解的范围，即x1xx2或x = qp(H,f,A,b,x1,x2,x0)其中x0为初始开始点。或x = qp(H,f,A,b,x1,x2,xo,N)其中N表示由A和b所定义的约束中，前N个为等式。例如，求解前述的二次规划问题可使用下面的命令：f=-2;-6;A=1,1;-1,2;b=2;2;H=2,-2;-2,4;disp(使用函数x=qp(H,f,A,b)x=qp(H,f,A,b)pausedisp(使用函数x=qp(H,f,A,b,X1,X2)x=qp(H,f,A,b,0,0,)pausedisp(使用函数x=qp(H,f,A,b,X1,X2,X0)x=qp(H,f,A,b,0,0,0,0) 运行结果如下：使用函数x=qp(H,f,A,b)x =