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

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

2、.图形的交互编辑 11第二部分数值计算 12§1 .方程求根121-1.牛顿迭代法121-2.图解法确定迭代的初始点 13§2 性方程组132-1.迭代法的收敛性 132-2.线性方程组白病态问题142-3.求解线性方程组15§3 插值和拟合163-1.Lagrange 插值163-2.代数多项式插值173-3.插值误差 173-4.分段线性插值183-5.数据的曲线拟合 18§4 .数值积分204-1.复合梯形求积公式204-2.复合 Simpson 求积公式20§5 .常微分方程的数值解法 215-1.Euler 方法215-2.改进的Eu

3、ler方法235-3.四阶龙格-库塔方法24习题27一、方程求根 27二、线性方程组 27三、插值与拟合 28四、数值积分29五、常微分方程30计算方法实验报告 31一、方程求根 31二、线性方程组31三、插值与拟合 32四、数值积分32五、常微分方程33III第一部分基本操作§1.Matlab的使用Matlab的使用方法有两种：（1）在 Matlab的命令窗口（ Matlab Command Windows）中直接输入命 令，即可得到结果；（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文件，最后点击该窗口上方菜单条中的File-Save保存文件。3 .例示一一 221）计算函数 f （x,为）=2x+7& - 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)在命令窗口中调用该函数键入

6、y=demof_1(2,3)显示y=40键入y=demof_1(3,-5)显示 y=196Bia2)计算一组数据x i(i = 1,2,?, ne S =1- 和S2 = 产nn 在编辑窗口中键入该函数的 M文件function s1,s2=demof_2(x)n=length(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)显示t1=3t2=3.3166或者，通过主程序调

7、用该函数i) 首先，在编辑窗口输入脚本M文件如下：x=input('x= =');s1,s2=demof_2(x);fprintf('s1=%12.5f s2=%12.5f ' ,s1,s2);ii) 保存该文件，备下次使用。iii)运行该脚本M文件，即在命令窗口点击File一Open,调用该文件(该文件在编辑窗口)；点击Debug一Run,在命令窗口显示待输入的数值变量提示x=输入一组数据，如 2,-3,0,9,13,55则显示结果s1=12.66667s2=23.409402 M文件程序的主要语句和主要函数2-1.Matlab的数字特征1 .数字类型Matl

8、ab中所有变量均为双精度，整型和实型之间没有差别。几个特殊的数字，如pi表示兀的值、inf表示8、eps表示计算机精度 2.2204 X10-16等。2 .算术运算符加减乘除乘方+-*/A3 .逻辑比较符大于小于大于等于小于等于andorif语句中的不等于if语句中的等于><>=<=&|=4 .数组变量（1）数组的表示 一维数组（等同于向量），例如：x=1,3,-4,7,21y=3：7（等同于 y=3,4,5,6,7）z=3:0.5:6（等同于 z=3,3.5,4,4.5,5,5.5,6）v='g','k','m'

9、二维数组（等同于矩阵），例如一个3X2数组：m=0.1,0.2,0.3;0.7,0.8,0.9（2 ）数组的元素对应于上述数组，x（2）=3, y（3）=5,m（2,1）=0.7等等。二维数组的整行或整列可以一个冒号表示，例如：m（1,:）=0.1 0.2 0.3,m（:,2）=0.20.8（3 ）数组的运算两数组的加（+）、减（-）运算符与向量的运算相同，而乘（.*）、除（./）、乘方（人）运算符不同。例如,对于数组a和b相加Z=a+ba和b的对应兀素相加相减Z=a-ba和b的对应兀素相减相乘Z=a.*ba和b的对应兀素相乘相除Z=a./ba和b的对应兀素相除乘方Z=a.Al.3a的所白兀素

10、的1.3次方2-2.主要语句1.if-end 语句每个if语句必须以一个end结束，二者一一对应。例如: if n= =2,price=17; end if n %price=15;else price=12; end if a>5,tap=10;elseif a<5,tap=-10;（如有必要，可多次使用elseif）else tap=1; end2. for-end 语句 for x=1:0.5:9 y=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:

11、0.1:10 y=sin（x）; if y<0,y=0;（循环中可以插入if-end语句）endend3. while-end 语句 i=1;c=0;x=-8,0,12,33,-6,5,-7; while i<=length(x)if x(i)<0,c=c+1; end i=i+1;endfprintf('C=%d',c);上述语句为统计数组x中小于0的元素个数。 while 1 ?if x>xlimit,break; end ?endwhile 1-end将可以无限循环下去，当条件 x>xlimit得到满足时，通过 break语句中止循环。4.

12、输入、输出语句输入(input)语句举例如下：z=input('input z=');输入一个数zp=input('Your name=','s');输入一个字符或字符串输出(fprintf)与C语言的输出语句类似，如fprintf('Volume=%12.5fn',vol);除120 格式，还可以输出 12.5e 格式、%d 格式、%s格式等。2-3.常用函数1 .常用教学函数 三角 函数： sin(x)、cos(x)、tan(x)、asin(x)、acos(x)、atan(x) 初等函数：abs(x)、sqrt(x)、roun

13、d(x)(四舍五入取整)、fix(x)(去尾数取整)、mod(x,y)(取余数)、 exp(x)(以e为底的指数)、log(x)(以e为底的对数)、log10(x)(以10为底的对数)2 .常用功能函数sum(x)求向量/数组元素之和 max(x)求向量/数组的最大值 min(x)求向量/数组的最小值rand(n)生成随机数，当n=1时返回一个随机数；n>1时，返回nXn随机矩阵 feval(f_name,x)计算以x为参数，名为f_name的函数，例如,s_name='sin'贝Uy=feval(s_name,x)等价于 y=sin(x)eval(s)执行字符串s所代表

14、的命令例如s='x=cos(pi/3)'则eval(s)等价于执行x=cos(pi/3)2-4.几个常用的命令1. cd 显示或改变当前目录cd 显示当前目录 cd c:matlabrllwork将此目录设置为当前目录2. dir 列出当前目录或列出指定目录中的文件dir显示当前目录中的内容 dir c:matlabrllbin显示此目录中的内容3. disp显示变量内容4. type列出指定文件的全部内容5. clear清除内容中的变量和函数6. home 清屏并将光标移至左上角7. exit 或 quit 退出 Matlab做矩阵的有关计算3-1.矩阵的输入1 .在命令窗口

15、内直接输入?1 2 3?对于a=?4 5 6?,可采用如下输入歹8 9?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阵例如， a=zeros(3)生成 3X3 的 0 阵 a=zeros(3,2)生成 3 >2 的 0 阵 b=zeros(size(a)生成与矩阵a维数相同的0阵格式同z

16、eros格式同zeros ones生成1阵eye生成单位阵4 .向量的输入行向量x=1,2,3,4,5列向量x=1,2,3,4,5'3-2.矩阵/向量的运算1 .加减c=a+bc=a-bc=a+3（矩阵亦可与数量运算）2 .乘c=a*bc=3*a3 .求逆c=inv（a） （a须为非奇异方阵）4 .除 矩阵的左除c=ab（等效于a-1*b） 矩阵的右除c=a/b（等效于b*a-1）5 .行列式的值b=det（a） （a须为方阵）6 .转置c=a'对于矢邱车A=?3O?6即3-3.矩阵的范数111.无穷范数即矩阵元素的绝对值按行相加的最大值7H 义： A = max1110010

17、可j=1输入命令p=norm(A,inf)显示p=192.1- 范数定义:IAn=max Ta11勺叼台F即矩阵元素绝对值按列相加的最大值输入命令p=norm(A,1)显示p=223.2- 范数定义：IA| 2二(ATA),入max 6A)表示为ATA的最大特征值输入命令p=norm(A,2)或 p=norm(A)显示p=15.86873-4.向量的范数对于向量V=1 3 6 -71 .无穷范数定义：VII= max v1输入命令 p=norm(V,inf)即 max(abs(V)显示p=7输入命令p=norm(V,-inf)即min(abs(V)显示p=11.1- 范数n定义：V 1 = 1

18、> i=1输入命令p=norm(V,1)显示p=173.2- 范数定义：V2输入命令p=norm(V,2)或 p=norm(V)显示p=9.74683-5.矩阵的条件数定义：A为非奇异矩阵，cond(A)r = A1 口 A r为矩阵A的条件数。其中r= 8,1,2 ；分别称为无穷 范数条彳数、1-范数条件数、2-范数条件数。对于矢邱车A =? ? ? 09. 4-171 .无穷范数条件数输入命令p=cond(A,inf)显示p=29.52632.1-范数条件数输入命令p=cond(A,1)显示p=30.31583.2-范数条件数输入命令p=cond(A) 或 p=cond(A,2)显示

19、p=19.05753-6.矩阵的特征值和特征向量定义：设A是nxn的矩阵，满足关系式AXi = AX，则(i=1,2,n)为A的特征值，向量 X为人对应的特征向量。?3 2?对于矢I阵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= 4001?0.8944?以上结果即为1=4、其对应的特征向量为 ?；1W4472? ,?- 0.7071?2=1、其对应的特征向量为 ?。? 0.7071 ?§4.Matlab 绘图4-1.绘图的基本

20、命令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,'linewidth',3)设置线的粗细 plot(x,y,'r')设置线的颜色3时，线变细，反之亦然，此数的默认值是0.5。'linewidth'和3控制所绘线的粗细，当用较小的数代替r表示画出红色的线，其余颜色设置见下表:红黄紫青绿兰白里 八、rymcgbwk plot(x,y,'+')设置绘出线的标记上述命令表示用标记

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

22、sin(x)5 .title 给图形加上标题例如 title('pressure Rario')6 .clf和cla清除图形窗口中的全部内容7 .subplot绘制图形阵列调用格式为 subplot(m,n,k)在图形窗口内绘制 mxn个图形阵列，k是图形窗口的序号。见下述例2。8.示例例1 输入下面的M文件并运行之clf; x=0:0.2:10;10.6y1=sin(x);y2=exp(-0.4*x);y3=y1.*y2;plot(x,y1,'r');hold onplot(x,y2,'g');plot(x,y3,'linewidth&

23、#39;,3);grid onxlabel('X');ylabel('Y1=sin(X) Y2=exp(-0.4*X)Y3=Y1*Y2');例2输入下面的M文件并运行之-oa6 4 2 D 2 4 6 0.OI.Q 。-O.-OI1*>=£>-(各 FUKMm UkiV 二clf; t=0:0.1:30;subplot(2,2,1);plot(t,sin(t);title('SubPlot No:1');SubPlot Nd:1xlabel('t');ylabel('Y=sin(t)');su

24、bplot(2,2,2);plot(t,t.*sin(t);title('SubPlot No:2');xlabel('t');ylabel('Y=t*sin(t)');subplot(2,2,3);plot(t,t.*sin(t)A2);title('SubPlot No:3');xlabel('t');ylabel('Y=t*sin(t)A2');subplot(2,2,4);plot(t,t.A2.*sin(t).A2);title('SubPlot No:4');xlabel

25、('t');ylabel('Y=tA2*sin(t)A2');垂直叠放的两个图形可通过如下方式绘制:subplot(2,1,1);plot(;subplot(2,1,2);plot(；20口-i J410102030tSuhPhrt Nd:330r-4C 01000Figure 14-2.图形的交互编辑图4-2-1所示为图形窗口左上角局部。1.改变曲线的属性用鼠标点击菜单栏中按钮移动鼠标使其指向所要编辑曲线上的任意一点。 按下鼠 标右键，曲线变成连续的黑点，并弹出下拉菜File Edit View (ns：art Tools Desktop口方0昌|恃I默RO&

26、#174;SubPlot No.1Window He 制400.5-0.6) |单。点击菜单中各盘命令,可以改变曲线的 颜色、线宽、标记和曲线类型等属性。图.'42T.i30o o o 2 4 a三产192.其他的编辑命令通过编辑窗口内的Tools菜单项亦可完成若干编辑功能。读者可以自行摸索。ii第二部分数值计算§1.方程求根1-1.牛顿迭代法x- 0.01 一例 1求解万程(0.01x+1)sin x-2+1- - 0.0096 = 0在 xo=4 附近的一个根。1 .编写牛顿迭代法的函数M文件function x=Newt_n(f_name,x0)x=x0;xb=x+1;

27、k=0;n=50;del_x=0.01;while abs(x-xb)>0.000001k=k+1;xb=x;if k>=n break;end; y=feval(f_name,x);y_driv=(feval(f_name,x+del_x)-y)/del_x;x=xb-y/y_driv;fprintf('k=%3.0f,x=%12.5e,y=%12.5e,yd=%12.5en',k,x,y,y_driv);end;fprintf('n Final answer=%12.6en',x);用默认文件名 Newt_n.m存盘。2 .编写待求方程的函数M文

28、件function y=eqn_1(x)y=(0.01*x+1)*sin(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.82170e+000,y=-9.55313e-009,yd=-8.88819e-001Final answer=2.

29、821702e+000 注意：Newt_n.m和eqn_1.m应保存在当前目录。1-2.图解法确定迭代的初始点x- 0.01万程(0.0仅+1)sin x- 0.0096 = 0具有多重根，若要求出某一区间(如xC|0,10|)的所x +1有根，可用图解法确定其初始点，然后再调用Newt_n函数解之，具体步骤为：1 .绘方程曲线该方程的解也就是? y产(0.0仅+1) sin xx- 0.010.0096的交点，故应绘出曲线y1和y2,观测两条曲线交点的坐16标值。输入M文件如下:clf;hold on; x=0:0.01:10;y1=(0.01*x+1).*sin(x);y2=(x-0.01

30、)./(x.A2+1)-0.0096;plot(x,y1,'r');plot(x,y2,'g');观测图形，两条曲线在 xe0,10有3个交点，交点x坐标值约为3、6.5和9.5,可将这些点作为牛 顿迭代法的初始点。2 .求解键入：Newt_n('eqn_1',3)显示：ans=2.8217键入：Newt_n('eqn_1',6.5)显示：ans=6.4351 键入：Newt_n('eqn_1',9.5)显示：ans=9.31892线性方程组2-1.迭代法的收敛性?10x -为-2x = 7.2 ?例 2对于? -

31、x+-10x+; 523x=34=28.3建立 Jocobi 迭代格式为：1)?x1k+1 ? ?0 ? k+1 ?为?=夕.19 k+1解法一:.2k -0.1 0.2?x1 ?0.72?00.2?>2 k?+ ?0.83?,判断其收敛性。0.20 ?xk ? ?0.42?0 0.1 0.2? ? 输入迭代矩阵 B= ?0.100.2?如 0.2 0 ?调用求矩阵特征值命令(eig命令)输入：eig(B)显示：ans= - 0.10000.3372-0.2372矩阵B的谱半径均为该矩阵的最大特征值，故B)= 0.3372<1,该迭代格式收敛。2) 解法二：调用求矩阵范数的命令(n

32、orm命令)若矩阵B的任意一种算子的范数小于1,则迭代收敛。故可计算B |、|闿1或B|2中的任意一种，只要小于1,即可判断出迭代收敛。输入：norm(B,1)显示：ans=0.4000即|B | =0.4000<1 ,该迭代格式收敛。亦可调用norm(B,2)或norm(B,inf)来判断。例3对于例2的线性方程组，若建立迭代格式为?x1k+1? ?010- 2?x1 ? ?- 8.3? 9o? ,? 9 o?2k+1 ?= ?- 105 ?为 k?+ ? 4.2?,其收敛性如何？?x：+1 ? ?5 -0.50 ?x3? ?-3.6?2-2.线性方程组的病态问题?11 ?x ? ?2

33、?例4分析方程组? ?= ? ?的病态性。.夕 1.000?%? ?2?输入该方程的系数矩阵A=?11 ?夕 1.0001?cond命令) 调用求矩阵条件数命令(输入：cond(A,1)显示： ans=4.0004e+004即cond(A)= A-1 |1? A |1 =4.0004 X104>>1 ,该方程组是病态方程组。亦可调用 cond(A)或cond(A,inf)命令来判断。2-3.求解线性方程组1 .调用x=Ab命令求解方程组对于线性方程组 AX=b （其中A为n3的矩阵、x和b均为n维列向量），有X=A -1b,该式与Matlab 的矩阵左除运算等效，即 X=Ab 。?

34、3x + 2为=-1 例5求解? X - & =1?3输入系数矩阵A= ?12 ?- 1?和吊数歹U向重b= ? ?1 ?-1?输入：X=Ab 显示：X=0.2000-0.80002 .利用LU分解求解方程组1) LU分解? 2x + x - 3%=2一一 ?对于-x + 3x + 2x = 0?1 -3?2? 3x +&-叙=1?2? ? 输入系数矩阵A=? 1 32 ?和常数项列向量b= ?0? ?3 1 -3?1?调用LU分解命令（lu命令）输入：l,u,p=lu(A)显示：1=1.000000-0.33331.000000.66670.10001.0000u=3.000

35、01.0000-3.000003.33331.000000-1.1000?123I2p=0 01其中l为单位下三角阵100（对角线元素均为1）, u为上三角；p为置换矩阵，它记录了选主元高斯消去法中行交换的信息。它们之间的关系为p*A=l*u2）利用LU分解求解线性方程组 对于AX=b ,用p左乘方程两边，得pAX=pb21uX=yly=pb通过对原方程组系数矩阵 操作步骤为：在命令窗口内输入插值节点坐标x=2,4,6,8,10y=5,7,9,11,13xi=2.5,5,7.3,9.1x、y以及待插值的数据xi如下:(a)(b)A的LU分解，将求解AX=b的问题转化为求解方程组(a)和方程组(

36、b)。其D求方程组(b)中的未知变量y输入：y=l(p*b) 显示：y=1.00000.3333 1.3000D求方程组(a)中的未知变量X输入：X=uy显示：X=-1.00000.4545-1.1818亦可将上述 、的命令合并为 X=出(l(p*b)§3.插值和拟合3-1.Lagrange 插值1 .编写Lagrange函数M文件function fi=lagrange(x,y,xi)fi=zeros(size(xi);npl=length(y);for i=1:nplz=ones(size(xi);for j=1:nplif i=j,z=z.*(xi-x(j)/(x(i)-x(j

37、);end;end;fi=fi+z*y(i);end注意： 插值节点坐标x、y为数组，待插值的xi可以是一个标量，亦可以是数组。Lagrange函数M文件必须保存在当前目录内。2 .调用Lagrange插值函数例6已知插值节点xk246810yk .、._ _ 5 _一 7,一、9 1113调用Lagrange函数输入：yi=lagrange(x,y,xi)显示：yi=5.5000 8.000010.3000 12.10003-2.代数多项式插值1 .代数多项式的表示方法已知(n+1)个节点信息 保/N=0,1,n可以构造一个n次代数插值多项式n ,n-1Pn = a X + an- 1X+?

38、+&x+a)如，watiabx5上峪用桃示的数舞丽,4乔素为多项式的系数，并且从左至右按降哥排列。例2 .代数多项式的插值计算分成二步进行：通过polyfit 命令，又已知(n+1)个节点唯一确定一个n次代数插值多项式；利用这个代数插值多项式计算待插点处的函数值。1)构造n次代数插值多项式 调用格式polyfit(x,y,n)其中数组x,y是已知节点坐标，n是代数插值多项式的阶数。利用该命令可求得n次代数多项式p(x) = c x+cn-1x n-1+?+cx+G 的系数。例如，输入节点数据：x=1.1,2.3,3.9,5.1y=3.887,4.276,4.651,2.117 调用命令

39、：a=polyfit(x,y,length(x)-1)显示：a=-0.20151.4385-2.74775.4370其中length(x)-1=4-1=3,故所得的多项式是：y=-0.2015x3+1.438或-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插值，可得到

40、与上面相同的结果，说明了插值多项式存在的唯一性。3 - 3.插值误差以下用作图方法观察插值的误差。x例7已知曲线f(>) = e.5- 2sin x,采用3个节点x1.12.35.1yf(i.i)f(2.3)f(5.1)构造插值多项式P2(x)作为f(x)的近似，通过作图的方法观察在x41.1,5.1的误差。运行以下脚本M文件clear,clf;hold on;x=1.1,2.3,5.1;y=exp(x/1.5)-2*sin(x);plot(x,y,'o');n=length(x)-1;coeff=polyfit(x,y,n);xp=1.1:0.05:5.1;yp=pol

41、yval(coeff,xp);plot(xp,yp,'r');yh=exp(xp/1.5)-2*sin(xp);plot(xp,yh,'k');xlabel('X');ylabel('f(x) P(x),data points:。');x图中黑线为曲线f(R = e，5 - 2sin x,红线为由3个节点成的2次插值多项式，圆圈表示插值节点。3-4.分段线性插值调用格式 yi=interpl(x,y,xi,'linear')其中数组x和y为已知节点坐标，xi是待插值的标量或数组，对应的 yi值通过线性插值算出。注意

42、，调用格式中的数组全部用列向量表示。例8已知数据xk135.281013yk46389.1-4求xi=4,5.5；11.7时线性插值函数之值。在命令窗口内输入列数组x、y和xix=1,3,5.2,8,10,13' y=4,6,3,8,9.1,-4' xi=4,5.5,11.7'调用分段线性插值命令yi=interp1(x,y,xi,'linear') 显示：yi=4.63643.53571.67673-5.数据的曲线拟合例9已知一组数据，作拟合曲线i012345xi0.10.40.50.70.70.9yi0.610.920.991.521.472.031

43、）方法一：求解正规方程 将上述一组数据标在坐标纸上，观察到各点在一条直线 附近，可设拟合曲线为y=co+cix其正规方程为：?G，?0 ) (%?)夜?_ ?0, f)?1?)?,?)? ? ?(? , f)?，其中?0 (X)=1,?1 (x)= x55(?0, ?o)=刀(x)?0(x)= B?1 = 6 i=0=055(?0, ?i)= (?1, ?0)= X? (x)?i(x)=汇1?x = 3.3 ""i=0i=055?, ?)= ?(x )?(x)=x2 = 2.21112 1 i 1 i 2 i% f)= Z?：(x)?y=归&=7.54 =0i=05

44、5(?1, f)= X? (x)?y= »? = 4.844 i=0=0在Matlab命令窗口中输入： a=6,3.3;3.3,2.21b=7.54,4.844'c=ab显示： c=0.28621.7646即所要求的拟合曲线为y=0.2862+1.7646x2)方法二：调用 polyfit(x,y,n)调用polyfit(x,y,n),可求得拟合代数多项式4为=c义+cn.1x *+?+Qx+G的系数，其中x,y为节点坐标数组，n为拟合多项式的阶数。对于例9中的节点数据，求拟合直线方程，操作如下：输入： x=0.1,0.4,0.5,0.7,0.7,0.9,y=0.61,0.9

45、2,0.99,1.52,1.47,2.03c=polyfit(x,y,1)显示： c=1.76460.2862即所要求的拟合曲线为y=1.7646x+0.28624数值积分4-1.复合梯形求积公式复合梯形求积公式为24I = f(x)dx-h?f(a+2 n-1 f(x ) + f(b)? Z2 k1)2)3)其中h=b-?k=1=a+k k=1,2,?n-1。n输入复合梯形公式的函数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

46、+k*h;t1=t1+feval(f_name,x(k);end;t=h*(fa+2*t1+fb)输入待积分的函数M文件例10中的被积函数为可勾=4y，其M文件如下，并以funjts.m文件名存盘。1 + x2function y=func_ts(x)y=4/(1+xA2);调用复合梯形公式若取 n=8,则输入：t=trapezia(0,1,8,'func_ts')显示：t=3.13904-2.复合Simpson求积公式复合Simpson求积公式为b h?1I = /(Rdx ?2nf(a + f(切+。2 f%1) + f %)?k=1其中 h= b- a,x n=a+ kh

47、rk= 1,2?,2n。21)输入复合Simpson公式的函数 M文件 function s=simpson(a,b,n,f_name) h=(b-a)/n;for k=1:2*nx(k)=a+k*h/2;end;fa=feval(f_name,a); fb=feval(f_name,b); s1=0;s2=0;for k=1:ns1=s1+feval(f_name,x(2*k-1);s2=s2+feval(f_name,x(2*k);end;s=h*(0.5*(fa-fb)+2*s1+s2)/3;2) 输入待积分的函数M文件此处仍以例10中的被积函数为例，且以 fun_ts.m文件名存盘。3

48、) 调用复合Simpson公式若取 n=4,则输入：s=simpson(0,1,4,'func_ts') 显示：s=3.1416通过对复合梯形公式和复合Simpson公式的比较，对于同一个被积函数，前者划分为8等分、计算结果为3.1390;后者划分为 4等分、计算结果为3.141。两种算法计算被积函数值的次数相近(复合梯形公式计算被积函数值9次、复合Simpson公式10次)，但复合Simpson公式的精度(代数精度3)高于复合梯形公式(代数精度1)。§5,常微分方程的数值解法5-1.Euler 方法1)计算公式,、r一？y = f(x y1,对于初值问题?, Eul

49、er公式y+1 = y + h?f(x, y)。以下通过仞题，说明Euler方法的?y(x) = %使用。例11质量m=70kg的人在t=0时刻从飞机上跳出，假设跳伞者垂直下降，且初始时刻垂直速度速度v(t1)=0,空气阻力F=cv2 (N), c=0.27kg/m,取步长h=0.1。试求速度，并绘出 00t 020s的解的图 形。2)分析由牛顿第二定律，mFv= - F + ng,式中重力加速度g=9.8m/s2。那么，建立初值问题为： dt22?V= f (t, V = - + g(5-1)?m? v(t) = 0根据例题要求，tea,t=0,20,将求解区间离散如下：11切131用I|1

50、I0=Q1=20其中等分数n= -b-ah对于所求的(5-1)式，采用Euler公式，可得下述的计算格式:? = M+=(0 1)hi = 2,3,?,n+1i=2,3?,n+1? v = 0?V= v-i + h* f(tiNi)3) 编写f(t,v)函数程序function f=descent(t,v) c=0.27,m=70,g=9.8;f=-c*vA2/m+g;4) 编写该问题的Euler法求解程序clear,clc;a=0,b=20,h=0.1;n=(b-a)/h;t(1)=0；v(1)=0;for i=2:n+1t(i)=t(1)+(i-1)*h;v(i)=v(i-1)+h*des

51、cent(t(i-1),v(i-1);end;for i=1:n+1fprintf('Time=%f,t(i);fprintf(' Velocity=%fn',v(i);end; plot(t,v);xlabel('Time (s)');ylabel('Velocity (m/s)');5) 运行求解例11的结果如下：Time=0.000000Time=2.000000Time=4.000000Time=6.000000Time=8.000000Time=10.000000Time=12.000000Time=14.000000Time=16.000000Time=18.000000Time=20.000000Velocity=0.000000Velocity=18.730327Velocity=32.989678Velocity=41.671830Velocity=46.250031Velocity=48.480593Velocity=49.525261Velocity=50.005441Velocity=50.224250Velocity=50.323563Velocity=50.3685585-2.改进的Euler