Matlab70使用说明-数值计算部分

1、代号9- 06A模块计算机辅助设计层次基础型hat lab 7.0使用说明（数值计算部分）华中科技大学国家机械基础课程教学基地2010年9月i第一部分基本操作11.Matlab 的使用11- 1.直接输入命令11- 2.用M 文件开发程序 12.M文件程序的主要语句和主要函数 22- 1.Matlab的数字特征22- 2.主要语句32- 3常用函数42- 4.几个常用的命令53矩阵的有关计算53- 1.矩阵的输入53- 2.矩阵/向量的运算63- 3.矩阵的范数63- 4.向量的范数73-5.矩阵的条件数73- 6.矩阵的特征值和特征向量 84.Matlab 绘图94- 1.绘图的基本命令94

2、- 2.图形的交互编辑11第二部分数值计算 12 .方程求根121- 1.牛顿迭代法121- 2.图解法确定迭代的初始点 132.线性方程组132- 1.迭代法的收敛性 132- 2.线性方程组的病态问题142- 3 .求解线性方程组153.插值和拟合163- 1.Lagra nge 插值163- 2.代数多项式插值173-3.插值误差 173-4.分段线性插值183- 5.数据的曲线拟合 184.数值积分204- 1.复合梯形求积公式204- 2.复合 Simps on 求积公式205.常微分方程的数值解法 215- 1.Euler 方法215- 2.改进的Euler方法235- 3.四阶龙

3、格-库塔方法24习题27一、方程求根 27二、线性方程组 27三、插值与拟合 28四、数值积分29五、 常微分方程30计算方法实验报告 31一、方程求根 31二、线性方程组 31三、插值与拟合 32四、数值积分32五、 常微分方程33III第一部分基本操作1.Matlab的使用Matlab的使用方法有两种：（1）在 Matlab的命令窗口（ Matlab Comma nd Win dows ）中直接输入命 令，即可得到结果；（2）在Matlab的编辑窗口（ Matlab Editor）内编写M文件，然后在命令窗口执行该文件，得到所需的结果。1- 1.直接输入命令在命令窗口中，直接输入命令。例如

4、：键入x=3+5显示x=8若键入x=3+5 ；不能直接显示结果，还须键入x方可显示x=81- 2.用M文件开发程序1. 设置当前目录M文件分为脚本 M文件（相当于主程序）和函数 M文件（相当于子程序或函数）。子程序或函数 调用时，须在当前目录内操作，故建议将用户创建的子目录设置为当前目录。这样，所有程序及命令的 操作都在用户子目录内，较为方便。以下操作只须设置一次，以后再进入Matlab系统时，当前目录即为已设定的目录。设置当前目录步骤为： 将鼠标移至Matlab图标，按右键弹出下拉菜单； 点击属性”弹出Matlab属性窗口； 将该窗口内 起始位置”栏中的路径更改为用户创建的子目录路径； 点击

5、应用”最后点击确定”。2. 如何编写程序对于Matlab命令窗口的上方菜单条，点击File宀New宀M-file，则弹出Matlab的编辑窗口。在编辑窗口内键入M文件，最后点击该窗口上方菜单条中的FiletSave保存文件。3例示2 21)计算函数 f (x, x ) = 2x+7x - 3x - 5x + 3禺-12 编辑窗口中输入该函数 M文件functiony=demof_1(x1,x2) y=2*x1A2+7*x2A2- 3*x1*x2-5*x1+3*x2-12;输入完毕后存盘(默认文件名为demof_1.m) 在命令窗口中调用该函数键入y=demof_1(2,3)显示y=40键入y=

6、demof_1(3,-5)显示 y=1962)计算一组数据x i(i = 1,2?, r)的S =, 和S2 = 1旦n n 在编辑窗口中键入该函数的M文件fun ctio n s1,s2=demof_2(x)n=len gth(x);s1=sum(x)/n;s2=sqrt(sum(x.A2)/n);这是一个具有多个返回值的函数，调用语句的左边须为向量。输入完毕后保存文件(默认文件名为demof_2.m)。 在命令窗口中调用函数输入数组x=1,2,3,4,5输入调用t1,t2=demof_2(x)显示t仁3t2=3.3166 或者，通过主程序调用该函数i) 首先，在编辑窗口输入脚本M文件如下：

7、x=i nput(x= =);s1,s2=demof_2(x);fprin tf(s1=%12.5fs2=%12.5f ,s1,s2);ii) 保存该文件，备下次使用。iii) 运行该脚本M文件，即在命令窗口点击FileOpen,调用该文件(该文件在编辑窗口)；点击Debug t Run,在命令窗口显示待输入的数值变量提示x=输入一组数据，如2,-3,0,9,13,55则显示结果s1=12.66667s2=23.409402. M文件程序的主要语句和主要函数2- 1.Matlab的数字特征1数字类型Matlab中所有变量均为双精度，整型和实型之间没有差别。几个特殊的数字，如pi表示n的值、in

8、f表示a、eps表示计算机精度 2.2204 X10-16等。2. 算术运算符加减乘除乘方+-*/A3. 逻辑比较符大于小于大于等于小于等于andorif语句中的不等于if语句中的等于=5,tap=10;elseifa5,t ap=-10;（如有必要，可多次使用elseif）else tap=1;end2. for-end 语句 for x=1:0.5:9y=xA2-5*x-3;end上述语句一次计算 x=1, 1.5, 2,，9时y的值。若改为for x=-2,0,15,则依次计算x=-2,0,1 5时的值。 for x=0:0.1:10y=si n（ x）;if y0,y=0;（循环中可以

9、插入if-e nd语句）endend3. while-end 语句 i=1;c=0;x=-8,0,12,33,-6,5,-7;while i=le ngth(x)if x(i)xlimit,break;end?endwhile 1-end将可以无限循环下去，当条件xxlimit得到满足时，通过 break语句中止循环。4. 输入、输出语句输入(input)语句举例如下：z=i nput(i nput z=);输入一个数zp=in put(Your name=,s);输入一个字符或字符串输出(fprintf)与C语言的输出语句类似，如fprintf(Volume=%12.5fn,vol);除%1

10、2.5f 格式，还可以输出 %12.5e 格式、%d 格式、%s格式等。2- 3.常用函数1. 常用教学函数 三角函数： sin(x)、cos(x)、tan(x)、asin(x)、acos(x)、atan(x) 初等函数：abs(x)、sqrt(x)、round(x)(四舍五入取整)、fix(x)(去尾数取整)、mod(x,y)(取余数)、exp(x)(以e为底的指数)、log(x)(以e为底的对数)、Iog10(x)(以10为底的对数)2. 常用功能函数sum(x)求向量/数组兀素之和max(x)求向量/数组的取大值min (x)求向量/数组的最小值rand(n)生成随机数，当n=1时返回一

11、个随机数；n1时，返回nxn随机矩阵feval(f_name,x)计算以x为参数，名为f_name的函数，例如，s_n ame=si n则y=feval(s_name,x)等价于 y=sin(x)eval(s)执行字符串s所代表的命令例如s=x=cos(pi/3)贝Ueval(s)等价于执行x=cos(pi2- 4 .几个常用的命令1.cd显示或改变当前目录cd显示当前目录 cdc:matlabrllwork将此目录设置为当前目录2.dir列出当前目录或列出指定目录中的文件dir显示当前目录中的内容dirc:matlabrllbin显示此目录中的内容3.disp显示变量内容4.type列出指定

12、文件的全部内容5.clear清除内容中的变量和函数6.home清屏并将光标移至左上角7.exit 或quit退出 Matlab3.矩阵的有关计算3- 1.矩阵的输入1. 在命令窗口内直接输入?123?对于a=?45?6?,可采用如下输入?789?a=1,2,3;4,5,6;7,8,9或 a=1 2 3;4 5 6;7 8 9或 a=1 2 34 5 67 8 92. 创建M文件对于大型矩阵直接输入不方便，可以采用创建M文件的方式，即在编辑窗口中输入该矩阵(输入格式同3-1的1)，然后保存在当前目录，每次进入Matlab时该矩阵已自动调入内容。3. 生成特殊矩阵的命令 zeros生成0阵例如，

13、a=zeros(3)生成 3X3 的 0 阵 a=zeros(3,2)生成 3 X 的 0 阵 b=zeros(size(a)生成与矩阵a维数相同的0阵17 ones生成1阵格式同zeroseye生成单位阵4.向量的输入格式同zeros行向量x=1,2,3,4,5列向量x=1,2,3,4,53- 2.矩阵/向量的运算1. 加减c=a+bc=a-bc=a+3（矩阵亦可与数量运算）2. 乘c=a*bc=3*a3. 求逆c=i nv（a） （a须为非奇异方阵）4. 除 矩阵的左除c=ab（等效于a-1*b） 矩阵的右除c=a/b（等效于b*a-1）5. 行列式的值b=det（a） （ a须为方阵）6

14、. 转置c=a3- 3.矩阵的范数对于矩阵A=s3 .1. 无穷范数定义:=maxj=1即矩阵元素的绝对值按行相加的最大值输入命令p=n orm(A, inf)显示p=192.1-范数定义:即矩阵元素绝对值按列相加的最大值输入命令显示3.2- 范数p=no rm(A,1)p=22定义：|A 2=代A),入max代A)表示为AA的最大特征值输入命令p=norm(A,2)或 p=norm(A)显示p=15.86873-4.向量的范数对于向量V=1 3 6 -71. 无穷范数即 max(abs(V)即 min(abs(V)定义：V输入命令 显示输入命令显示2.1- 范 数定义：V 1输入命令显示3.

15、2- 范 数定义：V II2=max v1 j jnp=n orm(V,i nf) p=7p=no rm(V,-i nf) p=1n=口i=1p=no rm( V,1)p=17I勺2i=1输入命令p=norm(V,2)或 p=norm(V)显示p=9.74683-5.矩阵的条件数定义：A为非奇异矩阵，cond(A)r = A UA r为矩阵A的条件数。其中r= ,1,2 ;分别称为无穷 范数条件数、1-范数条件数、2-范数条件数。?324?对于矩阵A = ?241?17?1.无穷范数条件数输入命令p=co nd(A,i nf)显示p=29.52632.1- 范数条件数输入命令p=co nd(A

16、,1)显示p=30.31583.2- 范数条件数输入命令p=co nd(A) 或 p=co nd(A,2)显示p=19.05753-6.矩阵的特征值和特征向量定义：设A是nxn的矩阵，满足关系式AXi = AX ,则(i=1,2,n)为A的特征值，向量 X为A对应的特征向量。?3 2?对于矩阵A= ?1 2?1. 特征值输入命令d=eig(A)显示d= 41(即 1=4、2=1)2. 特征向量和特征值输入命令v,d=eig(A)显示v= 0.8944-0.70710.44720.7071d= 400170.8944?以上结果即为1=4、其对应的特征向量为?；?0.4472?- 0.7071?2

17、=1、其对应的特征向量为 ?。? 0.7071 ?4.Matlab 绘图4- 1.绘图的基本命令3时，线变细，反之亦然，此数的默认值是0.5。l. plot画出点数据集合的曲线1) 基本格式plot(x,y) x,y是点的坐标数组 例如 x=0:0.2:10y=x.A3-12*x.A2-7*x+12Plot(x,y)2) 扩展格式 plot(x,y,l in ewidth,3)设置线的粗细linewidth和3控制所绘线的粗细，当用较小的数代替plot(x,y,+) plot(x,y,r)设置线的颜色红黄紫青绿白里八、rymcgbwkr表示画出红色的线，其余颜色设置见下表:设置绘出线的标记上星

18、形占八、十字形叉形菱形圆形五角形正方形三角形*+XDOPSA述命令表示用标记+绘出线，设置标记见下表:颜色与标记可以组合使用,例如plot(x,y, g)2. grid 画网格格式 grid ongrid off绘图形加上网格去掉网格3. hold保持图形格式 hold on保持图形取消此功能使用plot绘出一条曲线后，再绘出另一条曲线之前,hold off需用此命令将原曲线保持下来。4. xlabel和ylabel给坐标轴加标注给x和y坐标轴加上标注例如 xlabel(x)给x轴加上标注y=s in(x)ylabel(y=sin(x)给 y 轴加上标注5. title给图形加上标题例如 ti

19、tle(pressure Rario)6. clf和cla清除图形窗口中的全部内容7. subplot绘制图形阵列调用格式为subplot(m ,n ,k)在图形窗口内绘制 mxn个图形阵列，k是图形窗口的序号。见下述 例2。8.示例例1输入下面的M文件并运行之clf; x=0:0.2:10;y仁sin( x); y2=exp(-0.4*x); y3=y1.*y2; plot(x,y1,r);hold onplot(x,y2,g); plot(x,y3,li newidth,3);grid onxlabel(X);ylabel(Y1=si n(X)Y2=exp(-0.4*X)Y3=Y1*Y2)

20、;g-oa&96 4 2 D 2 4 6 01010 Q-o.-ol 罕世(XJTqf-d益吏塞暑clf; t=0:0.1:30; subplot(2,2,1);plot(t,si n( t); title(SubPlot No:1);SubPlot Nd:1xlabel(t);ylabel(Y=sin (t); subplot(2,2,2);plot(t,t.*si n( t); title(SubPlot No:2);0 6xlabel(t);ylabel(Y=t*sin (t); subplot(2,2,3);plot(t,t .*sin(t).八2); title(SubPlot No

21、:3);10-4C0SubPlot No:23020口=xlabel(t);ylabel(Y=t*sin( t)A2); subplot(2,2,4);plot(t,t.A2.*si n(t).A2); title(SubPlot No:4);xlabel(t);ylabel(Y=tA2*sin( t)A2);垂直叠放的两个图形可通过如下方式绘制:CMi _520SubPIni Nd:330r1000ID 20ISubPlot Nd: 4例2输入下面的M文件并运行之subplot(2,1,1);plot(;subplot(2,1,2);plot(;Figure 1File Edit view

22、Insert Tools： Desktop nfindow4-2.图形的交互编辑图4-2-1所示为图形窗口左上角局部。1.改变曲线的属性用鼠标点击菜单栏中按钮13,移动鼠标10.50使其指向所要编辑曲线上的任意一点。 按下鼠.06标右键，曲线变成连续的黑点， 并弹出下拉菜 单。点击菜单可以改变曲线的 颜色、线宽、标记和曲线类型等属性。2.其他的编辑命令SubPlot No：14-2-1 . I30aw 一严IrA通过编辑窗口内的Tools菜单项亦可元成右干编辑功能。读者可以自行摸索。第二部分数值计算.方程求根1- 1.牛顿迭代法x- 0 01 例1求解方程(0.01x+1)sin x- 0.0

23、096 = 0在xo=4附近的一个根。x +11 .编写牛顿迭代法的函数M文件function x=Newt_ n(f_n ame,x0)x=x0;xb=x+1;k=0;n=50;del_x=0.01;while abs(x-xb)0.000001k=k+1;xb=x;if k=n break;e nd; y=feval(f_ name,x);y_driv=(feval(f_ name,x+del_x)-y)/del_x;x=xb-y/y_driv;fprin tf(k=%3.0f,x=%12.5e,y=%12.5e,yd=%12.5en,k,x,y,y_driv);en d;fprin tf

24、(nFinal an swer=%12.6en,x);用默认文件名 Newt_ n.m存盘。2. 编写待求方程的函数M文件function y=e qn _1(x)y=(0.01*x+1)*si n(x)-(x-0.01)*(xA2+1)A(-1)-0.0096;用默认文件名eqn_1.m存盘。3. 求解步骤在命令窗口中键入 Newt_n(eqn_1,4)显示： k= 1,x=2.36795e+000,y=-1.03138e+000,yd=-6.31954e-001k= 2,x=2.92631e+000,y=3.48815e-001,yd=-6.24716e-001k= 6,x=2.82170

25、e+000,y=-9.55313e-009,yd=-8.88819e-001Final an swer=2.821702e+000 注意：Newt_n.m和eqn_1.m应保存在当前目录。1- 2.图解法确定迭代的初始点X 0 01 方程（0.01x+1）sin x-20.0096 = 0具有多重根，若要求出某一区间（如x|0,1d|）的所x +1有根，可用图解法确定其初始点，然后再调用Newt_ n函数解之，具体步骤为：1. 绘方程曲线该方程的解也就是? % ?y2=(0.01x+1) sin x=x- 0.01f +10.0096的交点，故应绘出曲线y1和y2，观测两条曲线交点的坐标值。输

26、入M文件如下:elf;hold on; x=0:0.01:10;y1=(0.01*x+1).*si n(x);y2=(x-0.01)./(x.A2+1)-0.0096;plot(x,y1,r);plot(x,y2,g);观测图形，两条曲线在 x0,10|有3个交点，交点x坐标值约为3、6.5和9.5，可将这些点作为牛顿迭代法的初始点。2. 求解键入：Newt_ n(eqn_1,3) 显示：an s=2.8217键入：Newt_n(eqn_1,6.5)显示:an s=6.4351 键入：Newt_n(eqn_1,9.5)显示:an s=9.3189线性方程组2- 1.迭代法的收敛性?10x -禺

27、-2x = 7.2? 例2对于? -（xH0x+! 52x=34=23.3建立Joeobi迭代格式为:?x k+1 ?00.10.2?x1 ?0.72?k+1 ?Xi?=?0.100.2?/?+?0.83?，判断其收敛性。? M? ? ?0.20.20 ?k ?0.42? 解法一：?00.10.2?输入迭代矩阵B=?0.100.2?0.20.20 ?1)调用求矩阵特征值命令(eig命令)输入：eig(B)显示：ans= - 0.10000.3372-0.2372矩阵B的谱半径PB)为该矩阵的最大特征值，故PB)= 0.33721，该迭代格式收敛。2) 解法二：调用求矩阵范数的命令(norm命令

28、)若矩阵b的任意一种算子的范数小于1，则迭代收敛。故可计算b |、b 或 b 12中的任意一种，只 要小于1,即可判断出迭代收敛。输入：n orm(B,1)显示：an s=0.4000即|B =0.40001，该方程组是病态方程组。亦可调用con d(A)或con d(A,i nf)命令 来判断。2-3.求解线性方程组1. 调用x=Ab命令求解方程组对于线性方程组 AX=b （其中A为n0的矩阵、x和b均为n维列向量），有X=A -1b,该式与Matlab 的矩阵左除运算等效，即 X=Ab。?3x + 2为=-1 例5求解? x - x? =1?32 ?、 输入系数矩阵A= ?和常数列向量?1

29、 -1? 输入：X=Ab显示：X=0.2000-0.80002. 利用LU分解求解方程组1) LU分解? 2x + x - 3xj = 2 ?对于- x + 3x + 2x = 0?123? 3x +x?-共=11 -3?2?2? ? 输入系数矩阵A=? 1 32 ?和常数项列向量b=刃？ ? 31 -3?1?调用LU分解命令lu命令）输入:l,u,p=lu(A)显示：1=1.000000-0.33331.000000.66670.10001.0000u=3.00001.0000-3.000003.33331.000000-1.1000p=001010100其中I为单位下三角阵（对角线元素均为

30、1），u为上三角；p为置换矩阵，它记录了选主元高斯消去法中行交换的信息。它们之间的关系为p*A=I*u2）利用LU分解求解线性方程组 对于AX=b，用p左乘方程两边，得/ pA=lu 贝U luX=pb24令uX=y(a)则ly=pb(b)通过对原方程组系数矩阵A的LU分解，将求解操作步骤为:求方程组(b)中的未知变量y输入：y=l(p*b) 显示：AX=bW问题转化为求解方程组 (a)和方程组(b)。其y=1.00000.3333 1.3000 求方程组(a)中的未知变量X输入：X=uy显示：X=-1.00000.4545-1.1818亦可将上述 、的命令合 并为 X=u(l(p*b)3.插

31、值和拟合3- 1 丄agrange 插值1.编写Lagrange函数M文件function fi=lagra nge(x,y,xi)fi=zeros(size(xi);n pl=le ngth(y);for i=1: nplz=on es(size(xi);for j=1: nplif i=j,z=z.*(xi-x(j)/(x(i)-x(j);en d;en d;fi=fi+z*y(i);end注意： 插值节点坐标x、y为数组，待插值的xi可以是一个标量，亦可以是数组。Lagrange函数M文件必须保存在当前目录内。2调用Lagrange插值函数例6已知插值节点8111013Xk2yk57,

32、,9 ,求xi=2.5、5、7.3、9.1时的函数值。在命令窗口内输入插值节点坐标x=2,4,6,8,10 y=5,7,9,11,13 xi=2.5,5,7.3,9.1x、y以及待插值的数据xi如下:调用Lagrange函数输入：yi=lagrange(x,y,xi)显示：yi=5.50008.000010.3000 12.10003- 2.代数多项式插值1. 代数多项式的表示方法已知(n+1)个节点信息(Xi,yji=0,1,n可以构造一个n次代数插值多项式n *1 cPn = a X + an- 1X+?+QX+a)如，多项式tiab=!2,上述多项式用一个数P=表示),4其元素为多项式的

33、系数，并且从左至右按降幕排列。例2. 代数多项式的插值计算分成二步进行：通过polyfit 命令，对已知(n+1)个节点唯一确定一个 n次代数插值多项式；利用这个代数插值多项式计算待插点处的函数值。1) 构造n次代数插值多项式 调用格式polyfit(x,y, n)其中数组x,y是已知节点坐标，n是代数插值多项式的阶数。利用该命令可求得n次代数多项式p(x) = c x+Cn一 1x n-1+?+CX+C0 的系数。例如，输入节点数据：x=1.1,2.3,3.9,5.1y=3.887,4.276,4.651,2.117 调用命令：a=polyfit(x,y,le ngth(x)-1)显示：a=

34、-0.20151.4385-2.74775.4370其中length(x)-1=4-1=3,故所得的多项式是：y=-0.20153+1.4385x - 2.7477x+5.43702) 代数多项式的插值计算调用格式 polyval(a,xi)其中a为插值多项式的系数数值，xi为待插点数组。若要计算xi=1.8,1.95,2.93,4.8 时的函数值，调 用该命令即可。输入：xi=1.8,1.95,2.93,4.8yi=polyval(a,xi)显示：yi=3.97714.05524.66853.1127若调用Lagrange插值，可得到与上面相同的结果，说明了插值多项式存在的唯一性。3- 3.

35、插值误差以下用作图方法观察插值的误差。x例7已知曲线f(x = e.5 - 2sin x,采用3个节点x1.12.35.1yf(1.1)f(2.3)f(5.1)构造插值多项式P2(x)作为f(x)的近似，通过作图的方法观察在x 1.1,5.1的误差。运行以下脚本M文件clear,clf;hold on;Xx=1.1,2.3,5.1;y=exp(x/1.5)-2*s in (x);plot(x,y,o);n=len gth(x)-1; coeff=polyfit(x, y,n);xp=1.1:0.05:5.1; yp=polyval(coeff,xp);plot(xp,yp,r);yh=exp(

36、xp/1.5)-2*si n(xp); plot(xp,yh,k); xlabel(X);ylabel(f(x) P(x),data poi nts:o);x图中黑线为曲线f(M = e.5 - 2sin x,红线为由3个节点构成的2次插值多项式，圆圈表示插值节点。3-4.分段线性插值调用格式yi=i nterpl(x,y,xi,li near)其中数组x和y为已知节点坐标，xi是待插值的标量或数组，对应的 yi值通过线性插值算出。注意，调用格式中的数组全部用列向量表示。例8已知数据Xk135.281013yk46389.1-4求 xi=4,5.5,11.7 时线性3插值函数之值。 在命令窗口

37、内输入列数组x、y和xix=1,3,5.2,8,10,13y=4,6,3,8,9.1,-4xi=4,5.5,11.7 调用分段线性插值命令yi=i nterp1(x,y,xi,li near) 显示：yi=4.63643.53571.67673-5.数据的曲线拟合例9已知一组数据，作拟合曲线i012345Xi0.10.40.50.70.70.9yi0.610.920.991.521.472.031) 方法一：求解正规方程 将上述一组数据标在坐标纸上，观察到各点在一条直线附近，可设拟合曲线为y=C0+C1X其正规方程为：?(?0,?0 ) (?0,?1 )金？_ ?0, f)? ?1?) ? 1

38、?)? ?= ?(? , 1)? ?其中洛(x)=1,?1 (x)= x55(?0, ?0 )= 7?(X)? (x) = 1 ? = 6i=0=055(?0, ?1)= (?1, ?0)= 口(X)?1(X)=刀?X = 3.30i=05?，?)=?(x )?(x )=5X2 = 2.21(?0, f)=壬?(x)?y= !?y= 7.54=0i=055(?1, f)= 口(X)?y= D?y = 4.844i=0=0在Matlab命令窗口中输入： a=6,3.3;3.3,2.21b=7.54,4.844c=ab显示： c=0.28621.7646即所要求的拟合曲线为y=0.2862+1.7

39、646x2) 方法二：调用 polyfit(x,y, n)p(X = c X, +cn- 1X 1+?+c(x+(S 的系数，其中 x,y例9中的节点数据，求拟合直线方程，操作如下：调用polyfit(x,y,n)，可求得拟合代数多项式为节点坐标数组，n为拟合多项式的阶数。对于输入： x=0.1,0.4,0.5,0.7,0.7,0.9, y=0.61,0.92,0.99,1.52,1.47,2.03 c=polyfit(x,y,1)显示： c=1.76460.2862即所要求的拟合曲线为y=1.7646x+0.28624.数值积分4- 1.复合梯形求积公式n）。例101用复合梯度公式求1=訴亍

40、X（此题的解析解为复合梯形求积公式为f(x)dx/ah?f(a+2n-1f(x ) + 他其中h= b a, xn?k=1=a+lh k=1,2,?n -1。s2=0;201)输入复合梯形公式的函数M文件function t=trapez ia(a,b,n,f_name) h=(b-a)/n; fa=feval(f_ name,a);fb=feval(f_ name,b); t1=0;for k=1: n-1 x(k)=a+k*h; t1=t1+feval(f_ name,x(k); en d;t=h*(fa+2*t1+fb)/2;2) 输入待积分的函数M文件例10中的被积函数为f(M =4

41、，其M文件如下，并以1 + x2;_ts.m文件名存盘。fun ctio n y=fu nc_ts(x)y=4/(1+xA2);3) 调用复合梯形公式若取 n=8,则输入：t=trapezia(0,1,8,func_ts)显示：t=3.13904- 2.复合Simpson求积公式复合Simpson求积公式为h?1?_3?2nf 佻k)?f(a + f(b)+ 口2 f%-1)k=1其中 h= k = a+ k -rk= 1,2?,2n。n21)输入复合Simpson公式的函数 M文件 function s=simps on( a,b ,n,f_n ame) h=(b-a)/n;for k=1:

42、2*nx(k)=a+k*h/2;en d;fa=feval(f_ name,a);fb=feval(f_ name,b);s仁0;for k=1: n33s仁 s1+feval(f_ name,x(2*k-1); s2=s2+feval(f_ name,x(2*k);en d;s=h*(0.5*(fa-fb)+2*s1+s2)/3;2) 输入待积分的函数M文件此处仍以例10中的被积函数为例，且以fun_ts.m文件名存盘。3) 调用复合Simpson公式若取 n=4,则输入：s=simpson(0,1,4,func_ts)显示：s=3.1416通过对复合梯形公式和复合Simps on公式的比较，对于同一个被积函数，前者划分为8等分、计算结果为3.1390；后者划分为4等分、计算结果为3.141。两种算法计算被积函数值的次数相近(复合梯形公式计算被积函数值9次、复合Simpson公式10次)，但复合Simpson公式的精度(代数精度3)高于复合梯形公式(代数精度1 )。5.常微分方程的数值解法5- 1.Euler 方法对