与数学实验（第二版）第5章数值运算ppt课件.pptx

《与数学实验（第二版）第5章数值运算ppt课件.pptx》由会员分享，可在线阅读，更多相关《与数学实验（第二版）第5章数值运算ppt课件.pptx（84页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、数学软件与数学实验（第二版）第5章 数值运算 电子课件第五章 数值运算数值计算是指在利用计算机解决数学问题的过程中，所涉及的变量必须有一个确定的数值（如实数或整数），计算的结果会有一定的误差。我国古代数学家刘微利用“割圆术”求得圆周率=3.14，祖冲之在前人的基础上求出在3.1415926与3.1415927之间，都用的是数值逼近的方法，这一结果比外国数学家获得同样结果早1000多年。秦九韶在数书九章中提出了高次代数方程的数值解法也是数值逼近。5.1 多项式多项式5.1.1 多项式的创建多项式的创建多项式的创建方法主要有以下几种：1.直接输入法2.通过矩阵的特征多项式来创建多项式3.由根向量创

2、建多项式5.1 多项式多项式5.1 多项式多项式2.通过矩阵的特征多项式来创建多项式设A为n阶方阵，则称的n次多项式f()=|E-A|为方阵A的特征多项式。求矩阵的特征多项式，由函数poly实现。调用格式为：lp=poly(A)求矩阵A的特征多项式系数，要求输入参数A是是nn的的方阵，输出参数p是包含是包含n+1个个元素的行向量，是A的特征多项式系数向量。5.1 多项式多项式5.1 多项式多项式3.由根向量创建多项式如果我们已知一多项式方程的根，可以用poly函数反求出这个根的多项式方程的系数。调用格式为：lp=poly(r)返回一个行向量，该行向量是以r为根的多项式系数向量。例例5-3 求以

3、-5，-3，4为根的多项式方程。r=-5-3 4p=poly(r)用向量p表示多项式不直观，可以用poly2sym函数将向量p转换为符号表达式的形式。poly2sym(p)5.1 多项式多项式由给定根向量创建多项式时应注意：如果希望创建实系数多项式，则根向量中的复数根必须共轭成对共轭成对。例例5-4 根据根向量r=-1+2i，-1-2i，0.2创建多项式。r=-1+2i，-1-2i，0.2;p=poly(r)%求多项式系数向量p=1.0000 1.8000 4.6000 -1.00005.1 多项式多项式5.1.2 多项式运算多项式运算1、求多项式的值求多项式的值可以有两种形式，对应两种算法：

4、(1)以标量或矩阵中每个元素为计算单元的，对应函数为polyval;(2)以矩阵为计算单元的，进行矩阵式运算来求得多项式的值，对应函数为polyvalm。115.1 多项式多项式调用格式为：lpolyval(p,x)：求多项式p在在x点的值点的值，x可以是标量或矩阵，x是矩阵时，表示求多项式p在x中各元素中各元素的值。lpolyvalm(p,x)：求多项式p对于矩阵矩阵x的值，x可以是标量或矩阵。x如果是标量，求得的值与函数polyval相同，如果x是矩是矩阵则必须是阵则必须是方阵方阵。125.1 多项式多项式5.1 多项式多项式145.1 多项式多项式3.多项式多项式的的乘法和除法乘法和除法

5、多项式的乘法和除法实质就是多项式系数向量的卷积卷积（Convolution）和解卷（）和解卷（Deconvolution）运算。lc=conv(a,b)求多项式a和b的乘法，如果向量a的长度为m，b的长度为n，则c的长度为m+n-1。多项式的除法用函数deconv实现，此函数也是向量的卷积函数的逆函数。lb,r=deconv(c,a)求多项式c除以a的商b与余项r。155.1 多项式多项式165.1 多项式多项式4、多项式的微积分、多项式的微积分Matlab中多项式导数的函数为polyder，求多项式不定积分的函数为polyint。两个函数的调用格式为：lpolyder(a)求系数行向量为a的

6、多项式的导数。lpolyint(a)求系数行向量为a的多项式的不定积分。175.1 多项式多项式185.1 多项式多项式5.1 多项式多项式 在MALAB中，两个多项式之比用部分分式展开的函数为residue，有两种调用方法：lr,p,k=residue(b,a)求多项式之比b(x)/a(x)的部分分式展开，输出参数r为留数，p为极点和k为直项。lb,a=residue(r,p,k)从部分分式得出多项式表达式b(x)和和a(x)的系数向量的系数向量，结果为对于表达式分母的归一形式。205.1 多项式多项式5.1 多项式多项式5.2 线性方程组求解5.2.1 齐次线性方程组的解法齐次线性方程组的

7、解法对于齐次线性方程组AX=0而言，可以通过求系数矩阵系数矩阵A的秩的秩来判断解的情况：(1)如果系数矩阵的秩=n（方程组中未知数的个数），则方程组只有零解只有零解。(2)如果系数矩阵的秩n，则方程组有无穷无穷多多个个解解，可以利用函数null(A)求出它的基本解系求出它的基本解系。22235.2 线性方程组求解5.2.2 非齐次线性方程组的非齐次线性方程组的解法解法对于非齐次线性方程组AX=b而言，则要根据系数矩阵系数矩阵A的秩和的秩和增广矩阵增广矩阵B=A b的秩和未知数个数的秩和未知数个数n的关系的关系，才能判断方程组AX=b的解的情况。（1）如果系数矩阵的秩=增广矩阵的秩=n，则方程组

8、有唯一解。（2）如果系数矩阵的秩=增广矩阵的秩n，则方程组有无穷多解。（3）如果系数矩阵的秩增广矩阵的秩，则方程组无解。245.2 线性方程组求解5.2 线性方程组求解5.2 线性方程组求解R=rank(A);B=A b;Rr=rank(B);if R=Rr&R=n%n为未知数的个数，判断是否有唯一解x=Ab;elseif R=Rr&R A=1，2，5，10，25;diff(A)ans=1 3 5 155.3 数值微积分数值微积分2.数值微分(numerical differentiation)根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。通常用差商代替微商，或者用

9、一个能够近似代替该函数的较简单的可微函数（如多项式或样条函数等）的相应导数作为所求导数的近似值。给定一个函数y=f(x),求其数值微分可以利用差分来近似计算。5.3 数值微积分数值微积分5.3 数值微积分数值微积分5.3 数值微积分数值微积分5.3 数值微积分数值微积分5.3 数值微积分数值微积分5.3.2 数值积分数值积分积分是微积分学与数学分析里的一个核心概念。通常分为定积分和不定积分两种。直观地说，对于一个给定的正实值函数，在一个实数区间上的定积分可以理解为在坐标平面上，由曲线、直线以及轴围成的曲边梯形的面积值曲边梯形的面积值（一种确定的实数值）。在多数情况下，求某函数的定积分时，被积函

10、数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的。另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解5.3 数值微积分数值微积分5.3 数值微积分数值微积分lq=quad2d(fun,a,b,c,d)计算z=f(x,y)在平面区域axb和c(x)yd(x)上的二重积分，其中c,d可以是一个数，也可以是一个函数句柄。lQ=trapz(X,Y)用梯形法计算Y的数值积分，X为一个分割。说明:fun是被积函数，指定为函数句柄。5.3 数值微积分数值微积分方法二：方法二：用trapz函数。在在

11、命命令令行行窗窗口口依依次次输输入入下下面面的的命令：命令：x=0:0.1:1;y=log(1+x)./(1+x.2);J=trapz(x,y)J=0.27135.3 数值微积分数值微积分方法二：方法二：用quad2d函数 fun=(x,y)x.2.*exp(-y.2);xmin=(x)x;I=quad2d(fun,0,1,xmin,1)I=0.04405.4 插值和拟合问题提出：5.4 插值和拟合5.4 插值和拟合42MATLAB中的一维插值函数是interp1，其调用格式为：lyi=interp1(x,y,xi,method)，输入参数x为样本数据点的横坐标向量，y为纵坐标向量或矩阵，me

12、thod为插值方法选项，xi为插值点的横坐标，yi是在xi指定位置计算出的插值结果。method的值常用的有：nearest、next、previous、linear、spline、pchip、makima。若省略method则默认为 linear。5.4 插值和拟合5.4 插值和拟合(1)nearest，最近邻近点插值。插值点的值为最近样本数据点的值，速度快但精确度差，不平滑。(2)next，下一邻近点插值。插值点的值为后面一个样本数据点的值，速度快但精确度差，不平滑。(3)previous，前一邻近点插值。插值点的值为前面一个样本数据点的值，速度快但精确度差，不平滑。(4)linear，线

13、性插值。在两个数据点之间连接直线，根据给定的插值点计算出它们在直线上的值，作为插值结果，速度较快，精度一般，默认形式。(5)spline，三次样条插值。通过数据点拟合出三次样条曲线，根据给定的插值点计算出它们在曲线上的值，将此作为插值结果。速度最慢，但精度高，平滑性最好。5.4 插值和拟合(6)pchip，保形分段三次插值，通过三次多项式计算插值结果。速度较慢，精度高，平滑性好。(7)makima，修正Akima三次Hermite插值。在插入点插入的值基于次数最大为 3 的多项式的分段函数。产生的波动比 spline 小，但不像 pchip 那样急剧变平。lyi=spline(x,y,xi)，

14、三次样条插值，其中输入、输出参数含义同上。lyi=pchip(x,y,xi)，分段三次Hermite插值，其中输入、输出参数含义同上例例5-18 已知正弦函数值如下表，试用四种插值方法估算x=55时正弦值。x0/6/3/22/35/6y01/23/213/21/20clearx=0,pi/6,pi/3,pi/2,2*pi/3,5*pi/6,pi;y=0 1/2 sqrt(3)/2 1 sqrt(3)/2 1/2 0;method=nearest,linear,spline,pchip;xi=55*pi/180;for i=1:4 yq(i)=interp1(x,y,xi,methodi);en

15、dy siny=sin(xi)yq 5.4 插值和拟合5.4 插值和拟合2、二维插值二维插值与一维插值的基本思想相同，它是对两个自变量两个自变量的函数z=f(x,y)进行插值。lzi=interp2(x,y,z,xi,yi,method)在已知(x,y,z)三维栅格点基础上，在(xi,yi)这些点上用这些点上用method插值方法估计函数值插值方法估计函数值。二维插值有五种插值方法：(1)nearest，最邻近点插值。(2)linear，双线性插值，该方法是interp2的默认插值方法。(3)spline，三次样条插值。(4)cubic，二重立方插值。(5)makima，修正Akima三次He

16、rmite 插值。475.4 插值和拟合8583808487788285868880818886848588908783例例5-19 测得一铁板表面45个网格点的温度如下表，试用四种插值方法绘制铁板表面温度的分布曲面5.4 插值和拟合在在MATLAB中建立命令文件如下中建立命令文件如下：%文件名为ex5_19clearx=1:5;y=1:4;X,Y=meshgrid(x,y);Z=85,83,80,84,87;78,82,85,86,88;80,81,88,86,84;85,88,90,87,83;xi=1:0.1:5;yi=1:0.1:4;XI,YI=meshgrid(xi,yi);5.4

17、插值和拟合Z1=interp2(X,Y,Z,XI,YI);Z2=interp2(X,Y,Z,XI,YI,nearest);Z3=interp2(X,Y,Z,XI,YI,cubic);Z4=interp2(X,Y,Z,XI,YI,spline);stem3(X,Y,Z),title(样本数据)axis(1 5 1 4 75 90)view(-24 17)figure,surf(XI,YI,Z1),title(linear)；figure,surf(XI,YI,Z2),title(nearest)figure,surf(XI,YI,Z3),title(cubic)；figure,surf(XI,Y

18、I,Z4),title(spline)5.4 插值和拟合5.4 插值和拟合5.4.2 拟合拟合拟合与插值方法不同，在构造函数时，不要求构造的函数一定通过样本点，而只要求构造的函数与样本点间的误差的平方和最小。数据拟合的实现方法有两种：(1)多项式拟合用最小二乘法求出拟合多项式的系数，在MATLAB中调用函数polyfit实现。lp=polyfit(x,y,n)，应用最小二乘法求出n阶拟合多项式p(x)，即用p(x)逼近f(x)。x、y为数据的横、纵坐标向量，n为拟合多项式的阶数，输出参数p为多项式p(x)的系数向量。525.4 插值和拟合例例5-20 测得铜导线在不同温度T()时的电阻如下表，

19、求电阻R与温度 T的近似函数关系。由于测得的数据接近一条直线，故我们可以用一条直线逼近函数。i0123456T()19.125.030.136.040.045.150.0R76.3077.8079.2580.8082.3583.9085.105.4 插值和拟合在在命令行窗口依次输入下面命令：命令行窗口依次输入下面命令：t=19.1 25.0 30.1 36.0 40.0 45.1 50.0;r=76.30 77.80 79.25 80.80 82.35 83.90 85.10;p=polyfit(t,r,1)%n=1根据拟合函数，我们可以推测铜导线电阻为0时的温度：T=roots(p)T=-2

20、42.1373%大约在零下242度左右。5.4 插值和拟合(2)矩阵除法拟合若多项式函数并未得出适合您数据的满意模型，您可尝试使用带非多项式项的线性模型。首先通过样本点的判断，构造一包含待定系数的函数表达式y=f(x)，然后将样本点数据代入y=f(x)，构造一线性方程组，最后解线性方程组，求出待定系数。5.4 插值和拟合在在MATLAB中建立命令文件如下：中建立命令文件如下：%文件名为ex5_21x=1,2,3,4,5,6,7,8,9;y=1.3,0.5,2.6,5.6,7.2,8.4,8.5,9.1,9.2;e=ones(size(x),exp(-x),x.*exp(-x);%方程组的系数矩

21、阵a=ey%用左除法解方程组x2=1:0.1:10;%在1,10区间上进行分割y2=ones(size(x2),exp(-x2),x2.*exp(-x2)*a;%用拟合函数计算分割点上的函数值figureplot(x2,y2,k,x,y,r*)%用红色*表示样本数据点，黑色实线为拟合曲线5.4 插值和拟合5.5 微分方程求解微分方程求解微分方程是指描述未知函数的导数与自变量之间关系的方程。微分方程的应用十分广泛，物理中许多涉及变力的运动学、动力学问题，如空气的阻力为速度函数的落体运动等问题，很多可以用微分方程求解。在微分方程中，自变量的个数只有一个，称为常微分方程，自变量的个数为两个或两个以上

22、的微分方程，称为偏微分方程，其中导数实际出现的最高阶数，称为微分方程的阶。微分方程的解是一个符合方程的函数。只有少数简单的微分方程可以求得解析解。在无法求得解析解时，可以利用数值分析的方法找到其数值解。5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解lt,y=ode15s(odefun,tspan,y0,options)采用可变阶的数值微分公式算法（NDFs）,属于多步算法，精度在低到中之间。适用于刚性方程，若 ode45 失败或效率低下且怀疑面临刚性问题，可尝试 ode15s。此外，该函数用于求解微分代数方程(DAE)。lt,y=ode23s(

23、odefun,tspan,y0,options)采用基于经过改进的2阶 Rosenbrock 公式。属于单步算法，精度低。适用于刚性方程，在求解允许宽松容差的问题或者解变化很快的问题时，效率可能高于 ode15s。lt,y=ode23t(odefun,tspan,y0,options)采用自由插值的梯形规则，精度低，适用于中等刚性问题。5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解建立函数文件如下：建立函数文件如下：f

24、unction dy=myode(x,y)dy(1,1)=y(2);dy(2,1)=y(2)+x;在命令行窗口输入下面命令：在命令行窗口输入下面命令：x,y=ode45(myode,0,1,0,1)plot(x,y(:,1),-o,x,y(:,2),-*)5.5 微分方程求解微分方程求解4常微分方程求解函数参数的设置常微分方程求解函数参数的设置利用MATLAB所提供的一系列求解函数求解常微分方程时，可以通过options控制方程求解的相对误差、绝对误差等参数，它是一个结构体，通过odeset函数来设定。loptions=odeset(Name,Value,.)，用参数名和参数值来设定求解函数的

25、参数。loptions=odeset(oldopts,Name,Value,.)，用参数名和参数值修改求解函数原来的参数oldopts。loptions=odeset(oldopts,newopts)，用newopts修改求解函数原来的参数。lodeset，显示求解函数所有的参数名及参数值。常用的参数名有：RelTol（相对误差），AsbTol（绝对误差），OutputFcn（输出函数）等，更多的参数可以通过odeset命令了解。5.5 微分方程求解微分方程求解例例5-25 求解例5-23的微分方程组的解，设置输出函数为二维相位平面图。首先建立函数文件如下：首先建立函数文件如下：functio

26、n dy=myfun(t,y)dy(1,1)=-2*y(1)+y(2)+2*sin(t);%dy要写成列向量的形式dy(2,1)=10*y(1)-9*y(2)+9*(cos(t)-sin(t);在命令行窗口输入下面命令：在命令行窗口输入下面命令：option=odeset(RelTol,1e-6,OutputFcn,odephas2)t,y=ode45(myfun,0,9,2,3,option);xlabel(y)ylabel(y的一阶导数)5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解MATLAB提供了bvp4c和bvp5c函数求解常微分方程的边值问题，其调用格式如下：ls

27、ol=bvp4c(odefun,bcfun,solinit)lsol=bvp4c(odefun,bcfun,solinit,options)lsol=bvp5c(odefun,bcfun,solinit)lsol=bvp5c(odefun,bcfun,solinit,options)说明说明：(1)bvp4c和bvp5c在大多数情况下，可以互换使用。它们之间的主要区别在于 bvp4c 实现四阶公式，而 bvp5c 实现五阶公式。(2)odefun是描述常微分方程组的函数句柄，其格式为dydx=odefun(x,y)，也可以包含未知参数，dydx=odefun(x,y,parameters)。5

28、.5 微分方程求解微分方程求解(3)bcfun是描述边界条件的函数句柄，其格式为res=bcfun(ya,yb)或包含未知参数res=bcfun(ya,yb,parameters)(4)solinit是对方程的猜测解，solinit可由bcpinit函数得到。(5)options是一结构体，用来设置bcp4c函数的参数，options结构体可由函数bcpset函数获得。lsolinit=bvpinit(x,yinit,params)计算常微分方程的猜测解猜测解，其中x是初始网格的有序节点，在区间a,b上求解方程，x取linspace(a,b,10)就足够了，yinit是猜测解因变量y的取值，p

29、arams是未知参数的猜测解。5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解(3)编写函数求解常微分方程满足边界条件的解。在MATLAB中建立命令文件如下：clear,clcx=linspace(0,4,5);%x取值yinit=1;0;%y取值solinit=bvpinit(x,yinit);%求猜测解sol=bvp4c(odefun,bcfun,solinit);%求数值解plot(sol.x,sol.y(1,:)5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解lsol=pdepe(m,pdefun,icfun,bcfun,

30、xmesh,tspan)(1)m 是对称常量。(2)pdefun，定义要求解的方程函数句柄。(3)icfun，定义初始条件的函数句柄。(4)bcfun，定义边界条件的函数句柄。(5)xmesh，是 x 的空间值向量x0 x1.xn。(6)tspan，是 t 的时间值向量t0 t1.tf。xmesh 和 tspan 向量共同构成一个二维网格，pdepe 在该网格上计算解。要使用 pdepe 求解 PDE，必须定义 c、f 和 s 的方程系数、初始条件、解在边界处的行为以及在其上计算解的点网格。5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解创建函数如下：function c,f,s=pdex1pde(x,t,u,dudx)c=pi2;f=dudx;s=0;end(2)编写初始条件代码编写初始条件代码function u0=pdex1ic(x)u0=sin(pi*x);end5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解5.5 微分方程求解微分方程求解(6)绘制解函数的图像绘制解函数的图像 u=sol(:,:,1);surf(x,t,u)title(偏微分方程的数值解)xlabel(距离 x)ylabel(时间 t)5.5 微分方程求解微分方程求解