现代数值计算方法(MATLAB版)第3章(2)

第三章数值分析一．多项式例ex3_1.m 多项式的定义、求根、求导%多项式的显示、系数矩阵及多项式的值disp('y(x)=x^3+5*x^2-9*x+3') %显示多项式表达式p=[1 5 -9 3]; %多项式系数矩阵x=1; %x赋初值y1=polyval(p,x); %计算x点处多项式的值%多项式求根、求导r1=roots(p); %求多项式的根p1=poly(r1); %用根构造多项式dy1=polyder(p); %对多项式求导数%显示结果disp(p); %显示多项式系数矩阵disp(y1); %显示x点处多项式的值disp(r1); %显示多项式的根disp(p1); %显示用根构造多项式的系数矩阵disp(dy1); %显示多项式一阶导数例ex3_2.m 多项式的乘、除disp('a(x)=x^3+2*x^2+3*x+4');disp('b(x)=x^3+4*x^2+9*x+16');a=[1 2 3 4]; %多项式a系数矩阵b=[1 4 9 16]; %多项式b系数矩阵c=conv(a,b) %两个多项式相乘，实际是求两个向量的卷积(Convolution):%∑--+=kiikbiakc1)1()()([d,r]=deconv(c,b) %多项式的除法，实际是卷积的逆运算%(Deconvolution), d,r为商多项式与余多项式:%∑=-+= -kiikqjakrkc1)1()()()(例ex3_3.m 多项式拟合x=0:0.1:1;y=[2.1,2.3,2.5,2.9,3.2,3.3,3.8,4.1,4.9,5.4,5.8];n=5; %拟合多项式的阶数取5p=polyfit(x,y,n); %用5阶多项式拟合x、y向量给定的数据y1=polyval(p,x); %计算x点处多项式的值plot(x,y,'o',x,y1,'-'); %绘x、y点图和拟合后的x、y1点图legend('y(x)','y1(x)'); %图例标注例ex3_4.m阶数对拟合效果的影响x=0:0.5:10;y=sqrt(x)+3*sin(x);n=2;p=polyfit(x,y,n);p2=polyval(p,x);n=4;p=polyfit(x,y,n);p4=polyval(p,x);n=6;p=polyfit(x,y,n);p6=polyval(p,x);n=8;p=polyfit(x,y,n);p8=polyval(p,x);plot(x,y,'o',x,p2,'-',x,p4,'-',x,p6,'-',x,p8,'-'); legend('y(x)','n=2','n=4','n=6','n=8');二、插值例ex3_5.m一维插值x=0:1:2*pi;y=sin(x);xi=0:0.1:6.5;yi1=interp1(x,y,xi,'nearst'); %最临近插值yi2=interp1(x,y,xi,'linear'); %线形插值yi3=interp1(x,y,xi,'cubic'); %三次多项插值yi4=interp1(x,y,xi,'spline'); %三次样条插值plot(x,y,'o',xi,yi1,xi,yi2,xi,yi3,xi,yi4);legend('y=sin(x)','nearst','linear','cubic','spline'); 例ex3_6.m高维函数插值（不讲！）[x,y]=meshgrid(-3:0.3:3);z=peaks(x,y);[xi,yi]=meshgrid(-3:0.1:3);zi=interp2(x,y,z,xi,yi,'spline');%zi=interp2(x,y,z,xi,yi,'cubic');%zi=interp2(x,y,z,xi,yi,'linear');%zi=interp2(x,y,z,xi,yi,'nearest');%surf(x,y,z);mesh(x,y,z);hold on;%surf(xi,yi,zi+15);mesh(xi,yi,zi+15);hold off;axis tight;三、快速富叶变换与逆变换（不讲！）例ex3_7.m快速富叶变换t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+2*randn(size(t));subplot(2,1,1);plot(y(1:50));%---------------y=fft(y); %快速富叶变换f=0:500;subplot(2,1,2);plot(f,y(f+1));例ex3_8.m滤波clear;x=linspace(0,2*pi,64);s=5*sin(x)+2*sin(5*x)+randn(size(x));f=fft(s);f1(1:9)=f(1:9);f1(56:64)=f(56:64);s1=ifft(f1);subplot(2,2,1);plot(x,s);%---------------subplot(2,2,3);fx=0:63;plot(fx,f(fx+1));%---------------subplot(2,2,2);plot(fx,f1(fx+1));%---------------subplot(2,2,4);plot(x,s1);四、稀疏矩阵（不讲！）例ex3_9.ma=[0 0 0 50 2 0 01 3 0 00 0 4 0];s=sparse(a) %转为稀疏矩阵形式b=full(s) %转为全元素矩阵形式c=sparse([3 2 3 4 1],[1 2 2 3 4],[1 2 3 4 5],4,4)%创建稀疏矩阵，c=[I,J,S,m,n,] I、J--行下标、列下标,S—按列排列的所有非零元素构成%的向量,m、n—待生成的稀疏矩阵行、列维数例ex3_10.mc=sparse([3 2 3 4 1],[1 2 2 3 4],[1 2 3 4 5],4,4)nnz(c) %稀疏矩阵的非零元素总数nonzeros(c) %稀疏矩阵的非零元素数值[i,j,s]=find(c) %找出稀疏矩阵的所有非零元素，按列排列[n,m]=size(c) %稀疏矩阵行、列维数d=c+ones(4,4) %稀疏矩阵加1矩阵五、数值积分例ex3_11.m(1)建立函数fn1function y=fn1(x)y=exp(-x.*x)(2)对被积函数fn1进行积分s1=quad('fn1',0,1,0.001); %采用Simpson法计算积分%0,1: 下限、上限%0.001: 精度s2=quad8('fn1',0,1,0.001,1);%采用8样条Newton-Cutes%公式求数值积分%最后的1:显示积分过程，0：不显示积分过程六、微分方程的数值解先将高阶微分方程降阶处理，转换为一阶微分方程组，再利用ode45函求解。

第7章MATLAB科学计算¾方程求解¾概率统计¾插值、拟合¾数值微积分¾最优化求解其它常用的matlab 数值计算命令¾max,min¾mean, median¾sum 求和, prod 求积¾cumsum 求和, cumprod 求积¾std 标准方差¾corrcoef 相关系数¾sort 元素排序¾离散傅里叶变换fft,fft2,fftn__ifft第7章MATLAB 数值计算作业¾1.编写傅立叶变换的matlab 程序与matlab 自带的fft 进行比较，并分析冲击信号的傅立叶变换。

（若不了解冲击信号，可计算方波的傅里叶变换，方波幅度为1，周期为10，方波个数为10，占空比为0.5）。

∑=−−−=Nm Nk m j em f k F 1/)1)(1(2)()(π编写的DFT 函数：function X=mydft(x )N = length(x );W=exp(-2*i*pi/N);X=zeros(1,N);for k=1:NX(k )=sum(x .*W.^((0:N -1)*(k -1)));end∑=−−−=N m N k m j em f k F 1/)1)(1(2)()(πx = [0 0 0 0 0 1 1 1 1 1]; X = [x x x x x x x x x x]; y = mydft(X);plot(abs(y))y1=fft(X);plot(abs(y1));¾y = fftshift(mydft(X));¾>> plot(abs(y))第3章MATLAB符号计算¾Maple优势在于符号运算，¾Mathematic符号运算和数值计算均不差，图像处理或者数据可视化较差¾Matlab强项是数值计算和数据可视化，¾MathCAD各方面均弱一些，但易学。

实验3 MATLAB 数值计算（2）
目的和要求：
（1）了解多项式的运算。

（2）熟练掌握MATLAB二维曲线的绘制。

内容和步骤：
1．多项式的运算式的运算
（1）多项式的运算。

已知表达式G(x)=(x-4)(x+5)(x2-6x+9), 展开多项式形式；求导；并计算当x在[0，20]范围变化时G(x)的值；计算出G(x)=0的根。

多项式相乘conv;
求导polyder;
计算零点，即求根roots
解：
展开为多项式
求导：
求零点：
（2）多项式的拟合和插值。

x在[0，20]范围内，计算多项式y=x4-5x3-17x2+129x-180 的值y；并根据x和y进行二阶、三阶和四阶拟合。

并绘出拟合曲线。

对多项式y进行插值，计算在5.5处的值。

多项式拟合p=ployfit(x,y,n)
插值yi=interp1(x,y,xi,’method’)
2．绘制二维曲线
绘制的图形窗口分割为一行两列，窗口左上角画一正弦曲线，y=sin(2t),t:[0.2π];窗口右上角画3条衰减的单边指数曲线y=e-t, y=e-2t,和y=e-3t, t:[0,2]。

在图上添加标题，将3条曲线用不同的线型，并添加图例。

第3章 MATLAB 数值运算教学提示：每当难以对一个函数进行积分或者微分以确定一些特殊的值时，可以借助计算机在数值上近似所需的结果，从而生成其他方法无法求解的问题的近似解。

这在计算机科学和数学领域，称为数值分析。

本章涉及的数值分析的主要内容有插值与多项式拟合、数值微积分、线性方程组的数值求解、微分方程的求解等，掌握这些主要内容及相应的基本算法有助于分析、理解、改进甚至构造新的数值算法。

教学要求：本章主要是让学生掌握数值分析中多项式插值和拟合、牛顿-科茨系列数值求积公式、3种迭代方法求解线性方程组、解常微分方程的欧拉法和龙格-库塔法等具体的数值算法，并要求这些数值算法能在MATLAB 中实现。

3.1 多 项 式在工程及科学分析上，多项式常被用来模拟一个物理现象的解析函数。

之所以采用多项式，是因为它很容易计算，多项式运算是数学中最基本的运算之一。

在高等数学中，多项式一般可表示为以下形式：120121()n n n n n f x a x a x a x a x a −−−=+++++…。

当x 是矩阵形式时，代表矩阵多项式，矩阵多项式是矩阵分析的一个重要组成部分,也是控制论和系统工程的一个重要工具。

3.1.1 多项式的表达和创建在MATLAB 中，多项式表示成向量的形式，它的系数是按降序排列的。

只需将按降幂次序的多项式的每个系数填入向量中，就可以在MATLAB 中建立一个多项式。

例如，多项式43231529s s s s +−−+在MATLAB 中，按下面方式组成一个向量x = [1 3 -15 -2 9]MATLAB 会将长度为n +1的向量解释成一个n 阶多项式。

因此，若多项式某些项系数为零，则必须在向量中相应位置补零。

例如多项式41s +在MATLAB 环境下表示为y = [1 0 0 0 1]3.1.2 多项式的四则运算多项式的四则运算包括多项式的加、减、乘、除运算。

下面以对两个同阶次多项式MATLAB 基础及其应用教程·66··66·32()234a x x x x =+++，32()4916b x x x x =+++做加减乘除运算为例，说明多项式的四则运算过程。

数值计算MATLAB 数值计算第四章MATLAB1主要内容基本数据运算数据统计与分析数据插值与曲线拟合多项式计算数值微积分线性方程组求解非线方与问求解非线性方程与最优化问题求解2常微分方程的数值求解多项式计算N次多项式表示为–P(x)=a 0x n +a 1x n-1+a 2x n-2…a n-1x+a nMatlab中n次多项式用一个长度为n+1的行向量（系数向量）表示–[a a …a n-1a [01n 1n ]多项式的四则运算–多项式的加减运算»系数向量的加减运算要求次数相同不足时用“»要求次数相同，不足时用“0”补起——向量化处理»例54322()352756f x x x x x x =−+−++3()353g x x x =+−多项式乘法运算–函数conv(P1,P2)用于求多项式P1和P2的乘积。

这里，P1、P2是两个多项式系数向量 多项式除法运算–函数[Q,r]=deconv(P1,P2)用于对多项式P1和P2作除法运算。

其中Q 返回多项式P1除的商式的余式这里仍是多项式系数向量以P2的商式，r 返回P1除以P2的余式。

这里，Q 和r 仍是多项式系数向量。

–deconv 是conv 的逆函数，即有P1=conv(P2,Q)+r 。

5432−− 例–求f(x)+g(x)、f(x)-g(x)。

2()352756()353f x x x x x xg x x x =+++=+−–求f(x)×g(x)、f(x)/g(x)。

–f=[3,-5,2,-7,5,6];g=[3,5,-3];g1=[0,0,0,g];–f+g1%求f(x)+g(x)f+g1 %求f(x)+g(x)–f-g1 %求f(x)-g(x)–conv(f,g) %求f(x)*g(x)[]()求()/()商式送余式送4–[Q,r]=deconv(f,g) %求f(x)/g(x)，商式送Q，余式送r。

第一章习题1. 序列满足递推关系，取及试分别计算，从而说明递推公式对于计算是不稳定的。

n1 1 0.01 0.00012 0.01 0.0001 0.0000013 0.0001 0.000001 0.000000014 0.000001 0.0000000110-105 0.00000001 10-10n1 1.000001 0.01 0.0000992 0.01 0.000099 -0.000099013 0.000099 -0.00009901-0.010000994 -0.00009901 -0.01000099-1.00015 -0.01000099-1.0001初始相差不大，而却相差那么远，计算是不稳定的。

2. 取y0=28，按递推公式，去计算y100，若取（五位有效数字），试问计算y100将有多大误差？y100中尚留有几位有效数字？解：每递推一次有误差因此，尚留有二位有效数字。

3．函数，求f(30)的值。

若开方用六位函数表，问求对数时误差有多大？若改用另一等价公式计算，求对数时误差有多大？设z=ln(30-y)，，y*， |E(y)| 10-4z*=ln(30-y*)=ln(0.0167)=-4.09235若改用等价公式设z=-ln(30+y)，，y*， |E(y)|⨯10-4z*=-ln(30+y*)=-ln(59.9833)=-4.094074．下列各数都按有效数字给出，试估计f的绝对误差限和相对误差限。

1)f=sin[(3.14)(2.685)]设f=sin xyx*=3.14, E(x)⨯10-2, y*=2.685, E(y)⨯10-3,sin(x*y*)=0.838147484, cos(x*y*)=-0.545443667⨯(-0.5454) ⨯⨯10-2+3.14(-0.5454) ⨯⨯10-3|⨯10-2⨯10-2|E r(f)| ⨯10-2⨯10-2<10-22)f=(1.56)设f = x y ,x*=1.56, E(x)⨯10-2, y*=3.414, E(y)⨯10-3,⨯⨯⨯10-2⨯⨯⨯10-3|⨯⨯⨯10-2⨯⨯⨯10-3|=0.051|E r(f)| =0.01125．计算，利用下列等式计算，哪一个得到的结果最好，为什么？6．下列各式怎样计算才能减少误差？7. 求方程x2-56x+1=0的二个根，问要使它们具有四位有效数字，至少要取几位有效数字？如果利用伟达定理， 又该取几位有效数字呢？解一:若要取到四位有效数字,如果利用伟达定理,解二:由定理二,欲使x1,x2有四位有效数字,必须使由定理一知，∆至少要取7位有效数字。