《CHAPT4数值计算》由会员分享，可在线阅读，更多相关《CHAPT4数值计算（27页珍藏版）》请在金锄头文库上搜索。
1、Chapt4 Chapt4 数值计算数值计算4.1 插值与数据拟合插值与数据拟合4.2 多项式多项式4.3 微分与积分微分与积分4.4 线性方程与非线性方程线性方程与非线性方程4.5 常微分方程常微分方程14.1 插值与数据拟合插值与数据拟合u在分析试验数据中，常常要面临将试验数据作解析描述的任在分析试验数据中，常常要面临将试验数据作解析描述的任务，这个问题有插值和曲线拟合两种方法务，这个问题有插值和曲线拟合两种方法。u当数据量不够，需要补充，且认定已有数据可信时当数据量不够，需要补充，且认定已有数据可信时, , 通常利通常利用函数插用函数插值，插值是一个非常有价值的工具。值，插值是一个非常有
2、价值的工具。u在曲线拟合中，假定已知曲线的规律，作曲线的最佳逼近，在曲线拟合中，假定已知曲线的规律，作曲线的最佳逼近，但不需要经过所有的数据点但不需要经过所有的数据点；4.1.1 插值函数插值函数多项式插值原理多项式插值原理已已知知函函数数 f(x) 在在区区间间a,b上上有有定定义义且且已已知知，则则可可以以求求出出点点 xi （a xi b）上上的的函函数数值值 yi=f(xi), (i=0,1,2,n)。如果如果xi出界，则对应值为出界，则对应值为NaN21、曲线插值函数、曲线插值函数interp1一维数值插值一维数值插值t=interp1(x,y,x0,method) x、y：原始数据
3、点原始数据点，是两个等长的已知向量。是两个等长的已知向量。 x0为进行插值的数组为进行插值的数组，计算函数在计算函数在x0处的值。处的值。 method为插值算法为插值算法:该缺省的使用假定为线性插值该缺省的使用假定为线性插值 3当当x0为等距时可以用快速插值法，使用快速插值法的格式为为等距时可以用快速插值法，使用快速插值法的格式为*nearest、*linear、*spline、*cubic以上五种方法得出的数据值一个比一个精确以上五种方法得出的数据值一个比一个精确, ,而所需内存及计算时间也一个而所需内存及计算时间也一个比一个要大要长比一个要大要长. .qnearestnearest方法方
4、法: :将插值点的数据设为与之最接近的邻点的数值将插值点的数据设为与之最接近的邻点的数值. .执执行速度最快，输出结果为直角转折行速度最快，输出结果为直角转折qlinearlinear方法方法: :利用每两个相邻插值节点作线性插值利用每两个相邻插值节点作线性插值. .默认值，在默认值，在样本点上斜率变化很大样本点上斜率变化很大q样条插值方法：分段多项式的在连接处光滑连接。最花时间，样条插值方法：分段多项式的在连接处光滑连接。最花时间，但输出结果也最平滑但输出结果也最平滑qHermiteHermite方法：对插值函数，不仅要求它在节点处与函数同值，方法：对插值函数，不仅要求它在节点处与函数同值，
5、而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值。而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值。qCubicCubic：三次：三次插值，要求插值，要求x x0 0必须为等距。最占内存，输出结果与必须为等距。最占内存，输出结果与splinespline差不多差不多4例例1：已知：已知1900年到年到2010年每隔十年的数据如下：年每隔十年的数据如下：75.995 91.972 105.711 123.203 131.669 150.697179.323 203.212 226.505 249.633 256.344 267.893用插值法求用插值法求1995年的数据。年的数据。命令如
6、下：命令如下：T =1900:10:2010;Y=75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 256.344 267.893;Y1=interp1(T,Y,1995)设时间变量设时间变量T为一行向量，数据用为一行向量，数据用Y为一个行为一个行向量向量,x0=1995.Y1 = 252.9885plot(T,Y,1995,Y1,Or)t=interp1(x,y,x0,method)5例例2：已知：已知1900年到年到2010年每隔十年的数据如下：年每隔十年的数据如下：75.995 91
7、.972 105.711 123.203 131.669 150.697179.323 203.212 226.505 249.633 256.344 267.893用三种不同的插值方法进行比较。用三种不同的插值方法进行比较。T =1900:10:2010;Y=75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 256.344 267.893;T1=1900:5:2010Y1=interp1(T,Y,T1,nearest);Y2=interp1(T,Y,T1,linear);Y3=inter
8、p1(T,Y,T1,spline);plot(T,Y,-b,T1,Y1,-r,T1,Y2,xk,T1,Y3,-y);legend(原曲线原曲线,最近点插值最近点插值,线性插值线性插值,样条插样条插值值)；6注意：注意： 1.首先首先x0的取值范围不能超出的取值范围不能超出x的给定范围，否则，会给出的给定范围，否则，会给出“NaN”错误。错误。 2.其次，独立变量必须是单调的。即独立变量在值上必须总是增加的其次，独立变量必须是单调的。即独立变量在值上必须总是增加的或总是减小的。或总是减小的。 3.最后，不同的插值方法得到的结果是不相同的，所以，在确定方法最后，不同的插值方法得到的结果是不相同的，
9、所以，在确定方法之前，要考虑该插值的应用领域以及插值的曲线形状。之前，要考虑该插值的应用领域以及插值的曲线形状。7例例3 某观测站测得某日某观测站测得某日6:00时至时至18:00时之间每隔时之间每隔2小时的室内外温小时的室内外温度度()，用，用3次样条插值分别求得该日室内外次样条插值分别求得该日室内外6:30至至17:30时之间每时之间每隔隔2小时各点的近似温度小时各点的近似温度()。设时间变量设时间变量h为一行向量，温度变量为一行向量，温度变量t为一个两列矩阵，其中第一列为一个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。存放室内温度，第二列储存室外温度。 h =6:2:18;t=
10、18,20,22,25,30,28,24;15,19,24,28,34,32,30;X1=6.5:2:17.5Y1=interp1(h,t,XI,spline) plot(h,t,-b,X1,Y1,-+r);legend(原曲线原曲线,样条插值样条插值)；YI = 18.5020 15.6553 20.4986 20.3355 22.5193 24.9089 26.3775 29.6383 30.2051 34.2568 26.8178 30.959482、曲面插值函数、曲面插值函数interp2二维数值插值二维数值插值zi=interp2(x,y,z,xi,yi,method) 其中其中X，
11、Y是两个向量，分别描述两个参数的采样点，是两个向量，分别描述两个参数的采样点，Z是与参数采是与参数采样点对应的函数值。样点对应的函数值。 xi，yi是两个向量或标量，描述欲插值的点，是两个向量或标量，描述欲插值的点，Zi是根据相应的插值是根据相应的插值方法得到的插值结果。方法得到的插值结果。 method的取值与一维插值函数相同。的取值与一维插值函数相同。X,Y,Z也可以是矩阵形式。也可以是矩阵形式。9clc;clear;x,y=meshgrid(-3:3); %创建创建x和和y矩阵矩阵z=peaks(x,y);surf(x,y,z);xi,yi=meshgrid(-3:0.5:3);对数据插
12、值加密对数据插值加密z1=interp2(x,y,z,xi,yi,nearest);z2=interp2(x,y,z,xi,yi,linear);z3=interp2(x,y,z,xi,yi,bicubic);figure,subplot(121); surf(xi,yi,z1);subplot(122);contour(xi,yi,z1);等值线等值线title(nearest);figure,subplot(121); surf(xi,yi,z2);subplot(122);contour(xi,yi,z2);title(linear);figure,subplot(121); surf(
13、xi,yi,z3);subplot(122);contour(xi,yi,z3);title(bicubic););10114.1.2 4.1.2 数据拟合数据拟合所谓数据拟合是求一个简单的函数所谓数据拟合是求一个简单的函数,例如是一个低次多项式例如是一个低次多项式,不要不要求通过已知的这些点求通过已知的这些点,而是要求在整体上而是要求在整体上“尽量好尽量好”的逼近原函的逼近原函数。这时数。这时,在每个已知点上就会有误差在每个已知点上就会有误差,数据拟合就是从整体上使数据拟合就是从整体上使误差误差,尽量的小一些。尽量的小一些。已知数据表已知数据表 x x1 x2 xmf(x) y1 y2 ym
14、最小12p=polyfit(x,y,n)p,s=polyfit(x,y,n)p,s,mu=polyfit(x,y,n)多项式拟合命令多项式拟合命令x=0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1;y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; p=polyfit(x, y, 2) 其中其中X，Y是两个向量，是已知的数据点。是两个向量，是已知的数据点。 n是阶次，如果我们选择是阶次，如果我们选择n=1，得到最简单的线性近似。通常称得到最简单的线性近似。通常称为线性回归。相反，如果我们选择为线性回归。相反，如果
15、我们选择n=2，得到一个得到一个2阶多项式。阶多项式。 p = -9.8108 20.1293 -0.0317polyfit 的输出是一个多项式系数的行向量。的输出是一个多项式系数的行向量。其解是其解是y = 9.8108x2 20.1293x0.0317多项式阶次的选择是有点任多项式阶次的选择是有点任意的意的，但是，但是选择的过选择的过小，拟小，拟合效果不好；阶次选择的过合效果不好；阶次选择的过高，虽然数据点上看到效果高，虽然数据点上看到效果好，数据点之间会好，数据点之间会出现有数出现有数据振荡的问题，阶数不宜过据振荡的问题，阶数不宜过高，小于高，小于5阶。阶。13xi=linspace(0
16、, 1, 100); z=polyval(p, xi);函数值函数值plot(x,y,o,x,y,xi,z,:) 为了将曲线拟合解与数据点比较，让我们把二者都绘成图。为了将曲线拟合解与数据点比较，让我们把二者都绘成图。 144.2 多项式多项式1、多项式表示、多项式表示多项式表达成一个行向量，该向量中的元素是按多项式降幂排多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。列的。 （表示中需要包含零系数的项表示中需要包含零系数的项）f(x)=anxn+an-1xn-1+a0可用行向量可用行向量 p=an an-1 a0表示表示例如例如:x412x325x11 p=1 -12 0 25
17、1115规定规定：多项式用行向量，根用列向量。多项式用行向量，根用列向量。给出多项式的根，使用给出多项式的根，使用polypoly函数也可以构造出相应的多项式函数也可以构造出相应的多项式2、roots 求多项式的根求多项式的根p=1 -12 0 25 11r=roots(p) r = 11.8142 1.7471 -1.0573 -0.50411611给出一个多项式的根（向量），也可以构造相应的多项式。给出一个多项式的根（向量），也可以构造相应的多项式。 3、poly 产生特征多项式系数向量产生特征多项式系数向量p=1 -12 0 25 11r=roots(p)p2=poly(r)p2= 1.
18、0000 -12.0000 -0.0000 25.0000 11.0000注意：注意：特征多项式一定是特征多项式一定是n+1n+1维的。维的。特征多项式第一个元素一定是特征多项式第一个元素一定是1 1。17a=1 2 3;4 5 6;7 8 0;p=poly(a)22给出一个矩阵，返回的是特征多项式的系数给出一个矩阵，返回的是特征多项式的系数 p = 1.00 -6.00 -72.00 -27.00p是矩阵是矩阵a的特征多项式的系数的特征多项式的系数r=roots(p)r = 12.12 -5.73 r是矩阵是矩阵a的特征值的特征值 -0.39p2=poly2sym(p,x)表达成多项式形式表
19、达成多项式形式p2=x3-6*x2-72*x-27p1=poly2str(p,x) 将多项式系数向量表示成将多项式系数向量表示成符号多项式符号多项式p1 = x3 - 6 x2 - 72 x - 27184、多项式运算、多项式运算对多项式加法，对多项式加法，MATLABMATLAB不提供一个直接的函数。如果两个多不提供一个直接的函数。如果两个多项式向量大小相同，标准的数组加法有效。项式向量大小相同，标准的数组加法有效。 多项式运算多项式运算命令命令乘法乘法（卷积）卷积）conv，convs除法除法 （反卷积）反卷积）deconv微分微分polyder(p)求求值值polyvald=a+b 19
20、1conv多项式乘法运算多项式乘法运算a=1 2 3;b=4 5 6;c=conv(a,b）p=poly2str(c,y)例例:a(x)=x2+2x+3; b(x)=4x2+5x+6;c = (x2+2x+3)(4x2+5x+6)c = 4 13 28 27 18p = 4 y4 + 13 y3 + 28 y2 + 27 y + 18202deconv多项式除法运算多项式除法运算a=1 2 3; c = 4 10 28 26 18 b,d=deconv(c,a)Q,r=deconv(P1,P2)返回多项式返回多项式P1除以除以P2的商式的商式r返回返回P1除以除以P2的余式的余式c=conv(
21、a,b)+db = 4 2 12d = 0 0 0 -4 -18c = 4 10 28 26 18213多项式微分多项式微分命令格式：命令格式：polyder(p): 求求p的微分的微分polyder(a,b): 求多项式求多项式a,b乘积的微分乘积的微分p,q=polyder(a,b): 求多项式求多项式a,b商的微分商的微分a=1 2 3 4 5; b=polyder(a)b = 4 6 6 4224多项式求值多项式求值可以求出多项式在某一点的值。可以求出多项式在某一点的值。 MATLAB提供了两种求多项式值的函数：提供了两种求多项式值的函数：polyval与与polyvalm，它们的输入
22、参数均为多项式系数向量它们的输入参数均为多项式系数向量P和自变量和自变量x。两者区别：前者是代数多项式求值，而后者是矩阵多项式求值。两者区别：前者是代数多项式求值，而后者是矩阵多项式求值。23已知多项式已知多项式x4+8x3-10，分别取，分别取x=1.2和一个和一个23矩阵为自变矩阵为自变量计算该多项式的值。量计算该多项式的值。a=1 8 0 0 -10; y1=polyval(a,1.2)x=2 3;4 5;y2=polyval(a,x)y1 = 5.8976y2 = 70 287 758 16151代数多项式求值代数多项式求值polyval函数用来求代数多项式的值，其调用格式为：函数用来
23、求代数多项式的值，其调用格式为：p 若若x为一数值，则求多项式在该点的值；为一数值，则求多项式在该点的值；p 若若x为向量或矩阵，则对向量或矩阵中的每个元素求其多为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。项式的值。Y=polyval(P,x)242矩阵多项式求值矩阵多项式求值polyvalm函数用来求矩阵多项式的值，其调用格式与函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同。相同，但含义不同。 polyvalm函数要求函数要求x为方阵，它以方阵为自变量求多项式的值。为方阵，它以方阵为自变量求多项式的值。Y=polyval(P,x)a=1 2 3;4 5 6;
24、7 8 0;p=poly(a)；x=pascal(3);创建一个正定矩阵创建一个正定矩阵polyvalm(p,x)；ans = -98.0000 -63.0000 -51.0000 -63.0000 -146.0000 -168.0000 -51.0000 -168.0000 -374.000025第三次上机作业第三次上机作业1.1.复习本次的课件，熟悉本次课程所提到的全部内容。复习本次的课件，熟悉本次课程所提到的全部内容。2.2.根据所学的内容，与本专业相联系，动手操作书上的根据所学的内容，与本专业相联系，动手操作书上的例题。例题。3.3.给出下表数据，求二次和三次拟合多项式给出下表数据，求
25、二次和三次拟合多项式。x x-1-1-0.75-0.75-0.5-0.5-0.25-0.250 025250.50.50.750.751 1y y-0.2209-0.22090.32950.32950.88260.88261.43921.43922.00032.00032.56452.56453.13343.13343.70613.70614.28364.28364 4、以函数、以函数 , , 为例，用为例，用 “nearest,linear,spline”“nearest,linear,spline”三插值方法进行插值，并对插值三插值方法进行插值，并对插值结果绘图进行比较。结果绘图进行比较。265.5.给出概论积分给出概论积分 的数据表的数据表x x0.460.460.470.470.480.480.490.49y y0.4846560.4846560.4937450.4937450.502750.502750.5116680.511668求（求（1）x=0.472时，该积分值等于多少？时，该积分值等于多少？ （2）当）当x为何值时积分值等于为何值时积分值等于0.5？6.（1）求方程）求方程 x3 4x + 5 = 0 的解的解 （2）求）求 x3 8x2 +6x 30=0的解的解 （3）求）求 的的“商商”及及“余余”多项多项式式。27
