Matlab的高级计算-2.ppt

（三）Matlab的高级数值计算 关系运算 逻辑运算 多项式计算 数值积分与微分 数据插值 曲线拟合 方程组求解 傅立叶分析 matlab语言把多项式表达成一个行向量， 该向量中的元素是按多项式降幂排列的。 f(x)=anxn+an-1xn-1+a0 可用行向量 p=an an-1 a1 +a0表示 (1) poly 产生特征多项式系数向量 特征多项式一定是n+1维的 3. 多项式运算 例:a=1 2 3;4 5 6;7 8 0; p=poly(a) p =1.00 -6.00 -72.00 -27.00 这是多项式p(x)=x3-6x2-72x-27的matlab描述方法， 可用： p1=poly2str(p,x) 函数文件，显示数学多项式的形式 p1 =x3 - 6 x2 - 72 x - 27 利用roots求多项式的根 r=roots(p) r = 12.1229 -5.7345 -0.3884 当然可用poly令其返回多项式形式p2=poly(r) p2 = 1.00 -6.00 -72.00 -27.00 vmatlab规定多项式系数向量用行向量表示，一组根用列向 量表示。 (2) 、多项式求根 求多项式零点: roots(p); 例1.求方程 x3 4x + 5 = 0 的解 P=1 0 -4 5; roots(P) ans = -2.4567 1.2283 + 0.7256i 1.2283 - 0.7256i 例2 求 x3 8x2 +6x 30=0 的解 P=1 -8 6 -30; r=roots(P) r = 7.7260 0.1370 + 1.9658i 0.1370 - 1.9658i (3) Polyval 计算系数为p的多项式在标量或向量x处的值 X = pascal(4) X = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 p = poly(X) p = 1 -29 72 -29 1 polyval(p,X) ans = 16 16 16 16 16 15 -140 -563 16 -140 -2549 -12089 16 -563 -12089 -43779 (4).conv多项式乘运算 例:a(x)=x2+2x+3; b(x)=4x2+5x+6; c = (x2+2x+3)(4x2+5x+6) a=1 2 3;b=4 5 6; c=conv(a,b) c = 4.00 13.00 28.00 27.00 18.00 p=poly2str(c,x) p = 4 x4 + 13 x3 + 28 x2 + 27 x + 18 (5).deconv多项式除运算 a=1 2 3; c = 4.00 13.00 28.00 27.00 18.00 d=deconv(c,a) d =4.00 5.00 6.00 d,r=deconv(c,a) 余数 c除以a后的商 conv(a,d) %因余数为零，可通过a,d的乘积返回 原多项式 ans = 4 13 28 27 18 (6).多项式微分 matlab提供了polyder函数多项式的微分。 命令格式： polyder(p): 求p的微分 polyder(a,b): 求多项式a,b乘积的微分 p,q=polyder(a,b): 求多项式a,b商的微分 例：a=1 2 3 4 5; poly2str(a,x) ans = x4 + 2 x3 + 3 x2 + 4 x + 5 b=polyder(a) b = 4 6 6 4 poly2str(b,x) ans =4 x3 + 6 x2 + 6 x + 4 k = conv(p,q) k,r = deconv(p,q) k = polyder(p) k = polyder(p,q) k,d = polyder(p,q) y = polyval(p,x) x = roots(p) 多项式运算小结 多项式运算中， 使用的是多项式 系数向量， 不涉及符号计算！ 4.数值微积分 4.1 数值积分 数值积分基本原理 求解定积分的数值方法多种多样，如简单的 梯形法、欧拉法、辛普生(Simpson)法、牛 顿柯特斯(Newton-Cotes)法等 基本思想都是将整个积分区间a,b分成n个子 区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b 。这样求定积分问题就分解为求和问题。 被积函数由一个表格定义 在MATLAB中，对由表格形式定义的函数关 系的求定积分问题用trapz(X,Y)函数。其中 向量X,Y定义函数关系Y=f(X)。 例 用trapz函数计算定积分。 命令如下： X=1:0.01:2.5; Y=exp(-X); %生成函数关系数据向量 trapz(X,Y) ans = 0.28579682416393 例：用数值积分法求 在x=0到 x=10之间所围面积，并讨论步长和积分 方法对精度的影响 建模： 将x的被积区间分为n段，各段长度为xi （i=1 ，2，n+1）, 欧拉法数值积分为 梯形法为 MATLAB程序 for dx=2,1,0.5,0.1 % 设不同步长 x=0:.1:10; y=-x.*x+115; plot(x,y,g),hold % 画出被积曲线 x1=0:dx:10; %根据不同步长取样点 y1=-x1.*x1+115; % 求取样点上的y1 n=length(x1); s=sum(y1(1:n-1)*dx; % 用欧拉法求积分 q=trapz(y1)*dx; % 用梯形法求积分 stairs(x1,y1), plot(x1,y1) % 画出欧拉法及梯形法的积分区域 dx,s,q,pause(1),hold off %给出不同步长下的结果 end 运行结果： 步长 欧拉法解 梯形法解 2 910 810 1 865 815 0.5 841.25 816.25 0.1 821.65 816.65 解析法得到的精确解为2450/3=816.6667 右图显示曲线的积分面积，在 曲线斜率为负的情况下，欧拉 法解偏大，梯形法偏小，精确 解位于两者之间 步长相同时，梯形法的精度比欧 拉法要高 变步长辛普生法 基于变步长辛普生法，MATLAB给出了quad函数来求定积分。 该函数的调用格式为： I,n=quad(fname,a,b,tol,trace) 其中fname是被积函数名。a和b分别是定积分的下限和上限。 tol 用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积 分过程，若取非0则展现积分过程，取0则不展现，缺省时取 trace=0。 需要事先建立函数文件fname 上例中，建立函数文件ex5f.m function y=ex5f(x) y=-x.*x+115; 应用quad函数 s=quad（ex5f,0,10） 运行结果：s = 816.6667 牛顿柯特斯法 基于牛顿柯特斯法，MATLAB给出 了quad8函数来求定积分。该函数的调 用格式为： I,n=quad8(fname,a,b,tol,trace) 其中参数的含义和quad函数相似，只是 tol的缺省值取10-6。该函数可以更精确 地求出定积分的值，且一般情况下函数 调用的步数明显小于quad函数，从而保 证能以更高的效率求出所需的定积分值 。 例： 分别用quad函数和quad8函数求定积分的 近似值，并在相同的积分精度下，比较函 数 的调用次数。 调用函数quad求定积分： format long; fx=inline(exp(-x); %定义函数 I,n=quad(fx,1,2.5,1e-10) 结果： I = 0.28579444254766 n = 65 调用函数quad8求定积分： format long; fx=inline(exp(-x); I,n=quad8(fx,1,2.5,1e-10) 运行结果 I = 0.28579444254754 n = 33 二重定积分的数值求解 使用MATLAB提供的dblquad函数就可 以直接求出上述二重定积分的数值解。 该函数的调用格式为： I=dblquad(f,a,b,c,d,tol,trace) 该函数求f(x,y)在a,bc,d区域上的二重定积 分。 参数tol，trace的用法与函数quad完全相同。 例 计算二重定积分 (1) 建立一个函数文件fxy.m： function f=fxy(x,y) f=exp(-x.2/2).*sin(x.2+y); (2) 调用dblquad函数求解。 I=dblquad(fxy,-2,2,-1,1) I = 1.57449318974494 4.2 数值微分 在MATLAB中，没有直接提供求数值 导数的函数，只有计算向前差分的函数 diff， 其调用格式为： 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)：计算矩阵A的n阶差分 ，dim=1时(缺省状态)，按列计算差分 ；dim=2，按行计算差分。 例 生成以向量V=1,2,3,4,5,6为基础的范得蒙矩 阵，按列进行差分运算。 命令如下：V=vander(1:6) DV=diff(V) %计算V的一阶差分 V = 1 1 1 1 1 1 32 16 8 4 2 1 243 81 27 9 3 1 1024 256 64 16 4 1 3125 625 125 25 5 1 7776 1296 216 36 6 1 DV = 31 15 7 3 1 0 211 65 19 5 1 0 781 175 37 7 1 0 2101 369 61 9 1 0 4651 671 91 11 1 0 5.数据插值 插值方法 一维插值的定义已知n个节点，求任 意点处的函数值。 分段线性插值 linear 多项式插值 cubic 样条插值 spline 最近邻插值方法 nearest y=interp1(x0,y0,x,method) 举例 多项式插值原理 已知 f(x) 在点 xi 上的函数值 yi=f(xi), (i=0,1,2,n) 求多项式 P(x)=a0 + a1x + anxn 满足: P(xk)= yk (k = 0,1,n) P10(t) f(t) f(x) 一维插值: yi = interp1(x, y, xi, method ) method nearest 最近邻邻点插值值 linear 线线性插值值 spline 样样条插值值 cubic 三次多项项式 插值值 例： x = 0:10; y = sin(x); xi = 0:.25:10; yi = interp1(x,y,xi);%默认线性插值 plot(x,y,o,xi,yi) %已知数据 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; %使用不同的方法进行插值运算 x = 1900:0.01:1990; yi_linear=interp1(t,p,x); yi_spline=interp1(t,p,x,spline); yi_cubic=interp1(t,p,x,cubic); yi_nearest=interp1(t,p,x,nearest); %绘制对比图形 subplot(2,1,1); plot(t,p,ko);hold on; plot(x,yi_linear,g,LineWidth,1.5);hold on plot(x,yi_spline,y,LineWidth,1.5); grid on; title(Linear Vs Spline); subplot(2,1,2); plot(t,p,ko);hold on; plot(x,yi_cubic,m,LineWidth,1.5);hold on plot(x,yi_nearest,k,LineWidth,1); grid on; title(Cubic Vs nearest); 对某城市人口数量统计统计插值 二维插值 高维插值的一种，主要应用于图像处理和数据的可视化 ，对两个变量的函数z=f(x,y)进行插值 zi=interp2(x, y, z, xi, yi, method) method 用法类似一维插值，只是参数linear 表示的是双线性插值方法 例：以二元函数 为例，首先获取基础 数据，而后使用二维插值得到更细致的数据，完成 绘图 %生成基础测量数据 x=-3*pi:3*pi; y=x; X,Y=meshgrid(x,y); R=sqrt(X.2+Y.2)+eps; Z=sin(R)./R; % eps正的极小值，防止出 现0/0 %进行二维插值运算 xi=linspace(-3*pi,3*pi,100); yi=linspace(-3*pi,3*pi,100); XI,YI=meshgrid(xi,yi); ZI=interp2(X,Y,Z,XI,YI,cubic) ; %meshgrid向量数据变为矩阵数据， 生成网格点 其他插值方法： 除了MATLAB内置的插值函数以外 ，还可以根据工程实践的需要，编写适 用于具体领域的插值函数文件来加以调 用， 如： 牛顿插值方法 拉格朗日（Lagrange）插值 Chebyshev多项式插值方法 6. 多项式拟和 在科学与工程领域，曲线拟合的主要功能是寻 求平滑的曲线来最好的表现带有噪声的测量数 据，从测量数据中找到变量函数的变化趋势， 最后得到曲线拟合的函数表达式 已知数据表 x x1 x2 xm f(x) y1 y2 ym 求拟合函数: (x) = a0 + a1x + + anxn 使得 离散数据的多项式拟合 达到最小 设 i=1,2,N 表示按拟合直线 求得的近似值，一般地说，它不同于 实测值，两者之差称 残差。显然，残差的大小是衡量 拟合好坏的重要标志，具体地说，可以采用下列三种准则： 使残差的最大绝对值为最小： 使残差的绝对值之和最小： 使残差的平方和为最小： 直线拟合 对于给顶的数据点 直线拟合问题可用数学语言描述如下： ,求作一次式 ,使总误差 最小。 有时候所给出数据点用直线拟合不合适，这 时可考虑用多项式拟合，用数学语言描述如 下： 对于给定的一组数据 ， ，寻 求作次多项式（ ） 使总误差 为最小。 曲线拟合 多项式拟合命令 P=polyfit(X,Y,n) 求出(最小二乘意义下)n次拟合多项式 P(x)=a0xn + a1xn-1 + + an-1x + an 计算结果为系数: P= a0, a1, , an-1, an 多项式求值命令: y1=polyval(P, x) 其中,P是n次多项式的系数,x是自变量的值,y1是多 项式在x处的值 例：如何预报人口的增长 人口的增长是当前世界上引起普遍关注的问题，并且我们会 发现在不同的刊物预报同一时间的人口数字不相同，这显然是 由于用了不同的人口模型计算的结果。 例如：1949年1994年我国人口数据资料如下： 年 份 xi 1949 1954 1959 1964 1969 1974 1979 19841989 1994 人口数 yi 5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8 建模分析我国人口增长的规律,预报1999年我国人口数。 介绍两个简单模型 模型一：假设：人口随时间线性地增加 拟合的精度: Q = ei 2 = (yi - a b xi)2, 误差平方和。 模型：y = 1.93 + 0.146 x 参数估计 观测值的模型： yi = a + b xi + ei ，i = 1,n 模型：y = a + b x 可以算出：a = 1.93， b = 0.146 模型二： 指数增长模型 （用简单的线性最小二乘法） 用Matlab软件计算得： a=2.33,b=0.0179 即: 程序如下： x=1949 1954 1959 1964 1969 1974 1979 1984 1989 1994; y=5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8 ; a=polyfit(x,y,1); x1=1949:10:1994; y1=a(2)+a(1)*x1; b=polyfit(x,log(y),1); y2=exp(b(2)*exp(b(1)*x1); plot(x,y,*) hold on plot(x1,y1,-r) hold on plot(x1,y2,-k) legend(原曲线,模型一曲线,模型二曲线) 程序执行后得到下面图形 结论的比较如下表： 年 份 xi 1949195419591964196919741979198419891994 人口 数 yi 5.466.778.19.19.810.311.311.8 模型 一值 5.245.976.77.438.168.99.6210.3611.0911.82 误 差 0.160.030-0.43-0.060.20.18-0.060.01-0.02 模型 二值值 5.556.066.627.237.98.649.4410.3111.2612.31 误误差-0.15-0.060.08-0.230.20.460.36-0.01-0.13-0.51 结果分析： 模型 I 2005年13.43亿， 2010年14.16亿 1. 误差Q1 = 0.2915 m时，此方程成为“超定”方程 当n A=1 2 1; 1 0 1; 1 3 0; b=2;3;8; x=linsolve(A,b) b是列向量！ 非线性方程的根 q Matlab 非线性方程的数值求解 fzero(f,x0)：求方程 f=0 在 x0 附近的根。 l 方程可能有多个根，但 fzero 只给出距离 x0 最近的一个 l fzer