MATLAB数值运算课件

1、MATLAB基础及其应用教程第3章: MATLAB数值运算第3章 提纲3.1 多项式3.2 插值和拟合3.3 数值微积分3.4 线性方程组解3.5 稀疏矩阵3.6 常微分方程的数值解polyadd(x,y)conv(a,b)q,r=deconv(a,b)poly(r) polyval(p,x) poly2sym(p)roots(a) interp1(x,y,xi,method)polyfit(x,y,n) diffcumsum trapz quad quadl rrefjacobi gseidel sor 3.1 多项式 3.1.1 多项式的表示和创建 在高等数学中，多项式一般可表示为以下形式

2、：在 MATLAB 中多项式用一个行向量表示，向量中的元素为该多项式的系数，按照降序排列。 如多项式 可以表示为向量 p=9 7 4 3。若多项式某些系数为0，则在向量对应位置补0。1、多项式相加、减（polyadd） a=1 2 3 4; b=1 4 9 16; c=polyadd(a,b)例1：完成两个同阶次多项式相加：例2：完成两个同阶次多项式相减： a=1 2 3 4; b=1 4 9 16; c=polyadd(a,-b)2、多项式乘法（conv）例3：完成两个同阶次多项式相乘： a=1 2 3 4; b=1 4 9 16; e=conv(a, b)3、多项式除法（deconv）例4

3、：求 a=1 6 20 50 75 84 64; b=1 4 9 16; f,r=deconv(a, b)例5：求 的“商”及“余”多项式。p1=conv(1,0,2, conv(1,4,1,1);p2=1 0 1 1;q,r=deconv(p1,p2);3.1.3 多项式求值、求根运算和多项式构造 p=1 2 -12 -1 7; z=polyval(p,3)例6：求函数 在 处的值 例7：估计矩阵多项式 在已知矩阵 处的值 ，其中X=1 2 1; -1 0 2; 4 1 2 X=1 2 1; -1 0 2; 4 1 2; P=1,0,-2,-1;Y=polyvalm(P,X)1、多项式求值（

4、polyval、polyvalm）例10：用多项式的根构造多项式 r=roots(p)； poly(r)3.2 插值和拟合插值：根据给定有限个已知数据（样本），构造一简单函数，使得函数通过全部的数据点。拟合：根据给定的数据，选取适当阶数的多项式逼近给定的数据，不要求多项式通过全部的给定数据。区别：拟合要找出一个曲线函数，插值仅求出插值数值即可。根据实验数据来确定函数的方法 x=0:10;y=cos(x); xi=0:.25:10; y0=cos(xi); %精确值 y1=interp1(x,y,xi); %线性插值结果 y2=interp1(x,y,xi,cubic); %三次方程式插值结果

5、y3=interp1(x,y,xi,spline); %三次样条式插值结果 plot(xi,y0,o,xi,y1,xi,y2,-.,xi,y3) %画出函数计算值和3种插值比较图例11：取余弦曲线上11个点的自变量和函数值点作为已知数据，再选取41个自变量点，分别用3种插值方法计算确定插值函数的值。3种插值方法比较图2、多项式拟合（polyfit） x=-2.8 -1 0.2 2.1 5.2 6.8; y=3.1 4.6 2.3 1.2 2.3 -1.1; p3=polyfit(x,y,3); p4=polyfit(x,y,4); p5=polyfit(x,y,5); xcurve=-3.5:

6、0.1:7.2; p3curve=polyval(p3,xcurve); p4curve=polyval(p4,xcurve); p5curve=polyval(p5,xcurve); plot(xcurve,p3curve,-,xcurve,p4curve,-.,xcurve,p5curve,-,x,y,*) grid on例12：对向量X=-2.8 -1 0.2 2.1 5.2 6.8和Y=3.1 4.6 2.3 1.2 2.3 -1.1分别进行3、4、5阶的多项式拟合。35阶多项式拟合曲线图3.3.1 微分和差分（diff）返回x对预设独立变量的一次微分值；返回x对独立变量t的一次微分值

7、；返回x对预设独立变量的n次微分值；返回x对独立变量t的n次微分值；3.3 数值微积分 S1=6*x3-4*x2+b*x-5; S2=sin(a); S3=(1-t3)/(1+t4); diff(S1); diff(S1,2); diff(S1,b); diff(S2); diff(S3);例13：计算方程式的微分： ， ，梯形法数值积分 x=linspace(0,pi,100); y=sin(x); t=trapz(x,y)辛普森数值积分 q=quad(sin,0,pi) quad(1./(x.3-2.*x-5),0,2) F= 1./(x.3-2.*x-5) ; quad(F,0,2):采

8、用低阶的自适应递归 Simpson 方法求积分q = quad(fun,a,b)，计算函数 fun 在区间a,b上的积分，其精确度为 1e-6。q = quad(fun,a,b,tol)，指定允许误差，指定的误差 tol 需大于 1e-6。q = quad(fun,a,b,tol,trace)，跟踪迭代过程，分别为计算函数值的次数、当前积分区间的左边界、步长和该区间内的积分值。MATLAB 中二重积分和三重积分分别由函数 dblquad() 和函数 triplequad() 来实现。首先介绍函数 dblquad()，该函数的基本格式如下：q = dblquad(fun,xmin,xmax,ym

9、in,ymax)，函数的参数分别为函数句柄、两个自变量的积分限，返回积分结果。q = dblquad(fun,xmin,xmax,ymin,ymax,tol)，指定积分结果的精度。triplequad() 函数的调用格式和 dblquad() 基本相同，在调用 triplequad() 函数时，需要六个参数指定积分限。 二重积分和三重积分3.4 线性方程组的数值解1、直接法矩阵相除法例14：求解线性方程组 a=1/2 1/3 1;1 5/3 3;2 4/3 5;b=1;3;2;c=ab例15：求解线性方程组 a=1 -1 1 -1;1 -1 -1 1;1 -1 -2 2;b=1;0;-0.5;

10、c=ab例16：求解线性方程组 a=1/2 1/3 1;1 5/3 3;2 4/3 5;1 2/3 1;b=1;3;2;2;c=ab2、迭代法Jacobi迭代法function tx=jacobi(A,b,imax,x0,tol)del=10-10;tx=x0;n=length(x0);for i=1:n dg=A(i,i); if abs(dg)del disp(diagonal element is too small); return endendfor k=1:imax for i=1:n sm=b(i); for j=1:n if j=i sm=sm-A(i,j)*x0(j); en

11、d end x(i)=sm/A(i,i); end tx=tx;x; if norm(x-x0) A=10 -1 2 0;-1 11 -1 3;2 -1 10 -1;0 3 -1 8; b=6 25 -11 15; tol=1.0*10-6; imax=10; x0=zeros(1,4); tx=jacobi(A,b,imax,x0,tol); for j=1:size(tx,1) fprintf(%4d %f %f %f %fn,j,tx(j,1),tx(j,2),tx(j,3),tx(j,4) endGauss-Seidel迭代法function tx=gseidel(A,b,imax,x

12、0,tol)del=10-10;tx=x0;n=length(x0);for i=1:n dg=A(i,i); if abs(dg)del disp(diagonal element is too small); return endendfor k=1:imax x=x0; for i=1:n sm=b(i); for j=1:n if j=i sm=sm-A(i,j)*x (j); end end x(i)=sm/A(i,i); end tx=tx;x; if norm(x-x0) A=10 -1 2 0;-1 11 -1 3;2 -1 10 -1;0 3 -1 8; b=6 25 -11

13、 15; tol=1.0*10-6; imax=10; x0=zeros(1,4); tx=gseidel(A,b,imax,x0,tol); for j=1:size(tx,1) fprintf(%4d %f %f %f %fn,j,tx(j,1),tx(j,2),tx(j,3),tx(j,4) endSOR迭代法function tx=gseidel(A,b,imax,x0,tol)del=10-10;tx=x0;n=length(x0);for i=1:n dg=A(i,i); if abs(dg)del disp(diagonal element is too small); retu

14、rn endendfor k=1:imax x=x0; for i=1:n sm=b(i); for j=1:n if j=i sm=sm-A(i,j)*x (j); end end x(i)=sm/A(i,i); x(i)=w*x(i)+(1-w)*x0(i); end tx=tx;x; if norm(x-x0) A=10 -1 2 0;-1 11 -1 3;2 -1 10 -1;0 3 -1 8; b=6 25 -11 15; tol=1.0*10-6; imax=10; x0=zeros(1,4);w=1.02; tx=sor(A,b,imax,x0,tol); for j=1:siz

15、e(tx,1) fprintf(%4d %f %f %f %fn,j,tx(j,1),tx(j,2),tx(j,3),tx(j,4) end3.5 稀疏矩阵1、稀疏矩阵的创建例20：创建45阶矩阵的稀疏矩阵 i=1 2 4; j=1 3 5; s=6 7 8; A=sparse(i,j,s)例21：创建矩阵A的对角稀疏矩阵 A=0 5 0 10 0 0 0 0 6 0 11 0 3 0 0 7 0 12 1 4 0 0 8 0 0 2 5 0 0 9;B,d=spdiags(A) s=spdiags(B,d,A)2、稀疏矩阵的存储例22：比较一个1010单位矩阵A及对应的稀疏矩阵B=spars

16、e(A)所需的存储量。 A=eye(10); B=sparse(A); whos A B稀疏矩阵的存储格式：（1）第一个数组以浮点格式存放数组中所有的非零元素；（2）第二个数组存放非零元素对应的行号，行号为整数；（3）第三个数组存放n+1个整形指针。3.6 常微分方程的数值解常微分方程求解思路：对求解区间进行剖分，然后把微分方程离散成在节点上的近似公式或近似方程，最后结合定解条件求出近似解。常微分方程求解方法：泰勒展开法、数值积分法（Runge-Kutta）、数值微分法（Adams）1、泰勒（Taylor）公式2、欧拉（Euler）法前向欧拉公式后向欧拉公式改进欧拉公式改进欧拉公式（预报）（校

17、正）例23：用前向欧拉法和改进的欧拉法解初值问题取步长h=0.1.（方程的解析解为 ）function f=exam1fun(x,y)f=y-2*x./y;f=f(:); xspan=0 1; y0=1; h=0.1;x1,y1=euler1(exam1fun,xspan,y0,h)x2,y2=euler2(exam1fun,xspan,y0,h)x=0:0.1:1; y=sqrt(1+2.*x); plot(x1,y1,-,x2,y2,*,x,y)3、Runge-Kutta法（1）二阶Runge-Kutta公式（2）三阶Runge-Kutta公式（3）四阶Runge-Kutta公式function f=exam2fun(x,y)f=-y-x*y.2;f=f(:); x1,y1=ode23(exam2fun,0:0.1:1,1)； x2,y2=ode45(