《数值计算与符号计算教材》由会员分享，可在线阅读，更多相关《数值计算与符号计算教材（83页珍藏版）》请在金锄头文库上搜索。
1、第4章 MATLAB数值计算与符号计算 4.1曲线拟合与插值运算 4.2 数值微积分 4.3 线性方程组求解 4.4 常微分方程的数值求解 4.5 MATLAB符号计算 4.6 级数 4.7 实验五 数值工具箱与符号工具箱的应用,4.1 曲线拟合与插值运算,1多项式的建立与表示方法 在MATLAB中，n次多项式用一个长度为n + 1的行向量表示，其元素为多项式的系数，按降幂排列，缺少的幂次项系数为0。 例如，多项式 在MATLAB中用向量p=1 -12 0 25 116表示。,2多项式的运算,(1) 多项式的加减运算 多项式的加减运算就是其所对应的系数向量的加减运算。相加减的多项式必须表示成相
2、同的次数，如果次数不同，应该把低次的多项式不足的高次项用0补足。,(2) 多项式的乘除运算 命令w=conv(u, v)表示多项式u和v相乘，例如在MATLAB中输入 u=1 2 3 4，v=10 20 30，c=conv(u, v) 返回 c= 10 40 100 160 170 120 conv指令可以嵌套使用，如conv(conv(a, b), c)。 命令q, r=deconv(v, u)表示u整除v。 向量q表示商，向量r表示余， 即有v=conv(u, q)+r。,(3) 多项式的导函数 对多项式求导数的函数有 k=polyder(p)，返回多项式p的导函数； k=polyder(
3、a, b)，返回多项式a与b的乘积的导函数； q, d=polyder(b, a)，返回多项式b整除a的导函数，其分子多项式返回给q，分母多项式返回给d。,(4) 多项式求值 MATLAB中提供了两种求多项式值的函数。 y=polyval(p, x)，代数多项式函数求值，若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。 Y=polyvalm(p, x)，矩阵多项式求值，要求x为方阵。设A为方阵，p代表多项式x3-5x2+8，那么polyvalm(p, A)的含义是 A*A*A-5*A*A+8*eye(size(A) 而polyval(p, A)的
4、含义是 A.*A.*A-5*A.*A+8*ones(size(A),例4.1 多项式P = x4-29x3 + 72x2 - 29x + 1，以4阶pascal矩阵为自变量分别用polyval和polyvalm计算该多项式的值。 在命令窗口输入如下命令： p= 1 -29 72 -29 1; X=pascal(4); A=polyval(p, X), B=polyvalm(p, X),返回 A= 16 16 16 16 16 15 -140 -563 16 -140 -2549 -12089 16 -563 -12089 -43779 B= 0 0 0 0 0 0 0 0 0 0 0 0 0
5、0 0 0,(5) 多项式的根 使用函数roots可以求出多项式等于0的根，根用列向量表示，其调用格式为 r=roots(p) 若已知多项式等于0的根，函数poly可以求出相应多项式，调用格式为 p=poly(r),例4.2 求多项式x4 + 8x3 - 10的根。 命令如下： A=1, 8, 0, 0, -10; x=roots(A) 返回 x= -8.0194 1.0344 -0.5075 + 0.9736i -0.5075 - 0.9736i 再输入 p=poly(x) 返回 p= 1.0000 8.0000 -0.0000 -0.0000 -10.0000,3曲线拟合,在MATLAB中
6、用polyfit函数来求得最小二乘拟合多项式的系数，再用polyval函数按所得的多项式计算所给出的点上的函数近似值。 polyfit函数的调用格式为 P, S=polyfit(X, Y, m) 函数根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样点的误差向量S。其中X, Y是两个等长的向量，P是一个长度为m + 1的向量，P的元素为多项式系数。,例4.3 用一个6次多项式在区间0, 2内逼近函数,MATLAB程序如下： x=linspace(0, 2*pi, 50); y=sin(x); P=polyfit(x, y, 6) %得到6次多项式的系数和误差 程序运行结果如下： P=
7、 0.0000 -0.0056 0.0874 -0.3946 0.2685 0.8797 0.0102,绘出sinx和多项式P(x)在给定区间的函数曲线，如图4.1所示。,图4.1 用6次多项式对正弦函数进行拟合,4数据插值,(1) 一维数据插值 在MATLAB中，实现一维数据插值的函数是interp1，被插值函数是一个单变量函数，其调用格式为 yi=interp1(x, Y, xi, method) 函数根据x, y的值，计算函数在xi处的值。x, y是两个等长的已知向量，分别描述采样点和样本值，xi是一个向量或标量，描述欲插值的点，yi是一个与xi等长的插值结果。method是插值方法，允
8、许的取值有linear、nearest、cubic、spline，分别表示线性插值、最近点插值、3次多项式插值、3次样条插值。 注意：xi的取值范围不能超出x的给定范围，否则，会给出“NaN”错误。,例4.4 向量t和p表示从19001990年的每隔10年的美国人口普查数据： t=1900:10:1990; p=75.995 91.972 105.711 123.203 131.669. 150.697 179.323 203.212 226.505 249.633; 根据人口普查数据估计1975年的人口，并利用插值估计从19002000年每年的人口数。,首先在命令窗口输入插值命令interp
9、1(t, p, 1975)，估计1975年的人口，返回结果如下： ans= 214.8585 再利用插值估计19002000年每年的人口数，并利用画图命令得出曲线分布。在MATLAB命令窗口输入如下命令语句： x=1900:1:2000; y=interp1(t, p, x, spline); %样条插值 plot(t, p, :o, x, y, -r) %带圆圈标记的虚线(:o)为普查数据，红实线 %(-r)为插值数据 所得曲线如图4.2所示。,图4.2 用一维数据插值得到的美国100年人口分布图,(2) 二维数据插值 在MATLAB中，提供了解决二维插值问题的函数interp2，其调用格式
10、为 Z1=interp2(X, Y, Z, X1, Y1, method) 其中，X, Y是两个向量，分别描述两个参数的采样点；Z是与参数采样点对应的函数值；X1, Y1是两个向量或标量，描述欲插值的点；Z1是根据相应的插值方法得到的插值结果。method的取值与一维插值函数相同。X, Y, Z也可以是矩阵形式。 同样，X1, Y1的取值范围不能超出X, Y的给定范围，否则，会给出“NaN”错误。,例4.5 利用插值运算对峰值函数peaks插入更多的栅格。 MATLAB程序如下： X, Y=meshgrid(-3:.25:3); Z=peaks(X, Y); XI, YI=meshgrid(-
11、3:.125:3); ZI=interp2(X, Y, Z, XI, YI); mesh(X, Y, Z), hold, mesh(XI, YI, ZI+15) hold off axis(-3 3 -3 3 -5 20) 程序运行结果如图4.3所示。,图4.3 用二维数据插值得到的peaks函数曲线图,例4.6 给定雇员数据如下： years=1950:10:1990; %工作年份 service=10:10:30; %服役时间，即在职时间 wage=150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767
12、 226.505 153.706 426.730 249.633 120.281 598.243; %工资 利用二维数据插值计算一个雇员在工作15年后在1975年所获得的工资。,MATLAB程序如下： w=interp2(service, years, wage, 15, 1975) 程序运行结果如下： w= 190.6287,4.2.1 数值微分 DX=diff(X)，计算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1, 2, , n-1。 DX=diff(X, n)，计算X的n阶向前差分，例如diff(X, 2)=diff(diff(X)。 DX=diff(A, n, dim
13、)，计算矩阵A的n阶差分，dim=1时(默认状态)，按列计算差分；dim=2时，按行计算差分。,4.2 数值微积分,例4.7 设x由0, 2间均匀分布的10个点组成，求sinx的13阶差分。 MATLAB程序如下： X=linspace(0, 2*pi, 10) Y=sin(X) DY=diff(Y) %计算Y的一阶差分 D2Y=diff(Y, 2) %计算Y的二阶差分，也可用命令diff(DY)计算 D3Y=diff(Y, 3) %计算Y的三阶差分，也可用diff(D2Y)或diff(DY, 2),程序运行结果如下： X= 0 0.6981 1.3963 2.0944 2.7925 3.49
14、07 4.1888 4.8869 5.5851 6.2832 Y= 0 0.6428 0.9848 0.8660 0.3420 -0.3420 -0.8660 -0.9848 -0.6428 -0.0000 DY= 0.6428 0.3420 -0.1188 -0.5240 -0.6840 -0.5240 -0.1188 0.3420 0.6428 D2Y= -0.3008 -0.4608 -0.4052 -0.1600 0.1600 0.4052 0.4608 0.3008 D3Y= -0.1600 0.0556 0.2452 0.3201 0.2452 0.0556 -0.1600,4.2
15、.2 数值积分,quad(filename, a, b, tol, trace) quadl(filename, a, b, tol, trace) 其中filename是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，默认时为10-6。trace控制是否展现积分过程，若取非0则展现积分过程，若取0则不展现，默认时取trace=0。,例4.9 用两种方法求 在MATLAB命令窗口输入 F=inline(1./(x.3-2*x-5); Q=quad(F, 0, 2) 或 Q=quadl(F, 0, 2),返回 Q= -0.4605 或者采用函数句柄 Q=quad(myfun, 0, 2) 也可以返回 Q= -0.4605,这里myfun.m为一个M文件： function y=myfun(x) y=1./(x.3-2*x-5);,I=dblquad(f, a, b, c, d, tol, trace) 该函数求f在a, b c, d区域上的二重定积分。参数tol, trace的用法与函数quad完全相同。,例4.11 计算二重定积分 首先建立一个函数文件fxy.m， function f=fxy(x, y) global ki; ki=ki+
