MATLAB在数值分析中的应用

MATLAB在数值分析中的应用何仁斌11MATLAB在数值分析中的应用11.1为什么要学习数值分析11.2多项式与插值11.3数据拟合11.4〔非〕线性方程组11.5微分方程组11.6积分11.1为什么要学习数值分析现实世界的问题可以归结为各种各样的数学问题●方程求根问题●解线性方程组的问题●定积分问题●常微分方程初值问题…等等在科学计算中常要遇到求解各种方程对于高次代数方程，由代数根本定理知多项式根的个数和方程的阶相同但对超越方程就复杂的多，如果有解，其解可能是一个或几个，也可能是无穷多个。11.1.1方程求根问题例如：高次代数方程x5

–3x＋1=0超越方程e-x

–cosx=0看似简单，但难求其精确解。方程求根问题由线性代数知识可知：当线性方程组Ax=b的系数矩阵A非奇异(即detA≠0)时，方程组有唯一解，可用克莱默法那么求解但它只适合于n很小的情况，而完全不适合于高次方程组。如用克莱默法那么求解一个n阶方程组，要算n+1个n阶行列式的值，总共需要n!(n-1)(n+1)次乘法。当n充分大时，计算量是相当惊人的解线性方程组的问题如用克莱默法那么求解一个n阶方程组，要算n+1个n阶行列式的值，总共需要n!(n-1)(n+1)次乘法。当n充分大时，计算量是相当惊人的一个20阶不算太大的方程组，大约要做1021次乘法，这项计算即使每秒1万亿次浮点数乘法计算的计算机去做，也要连续工作2000万亿年才能完成。当然这是完全没有实际意义的，故需要寻找有效算法11.1.2解线性方程组的问题由微积分知识知，定积分的计算可以使用牛顿——莱布尼兹公式：其中F(x)为被积函数f(x)的原函数。为何要进行数值积分？定积分问题原因之一：许多形式上很简单的函数，例如等，它们的原函数不能用初等函数表示成有限形式。定积分问题原因之二：有些被积函数的原函数过于复杂，计算不便。例如的一个原函数是定积分问题原因之三：f(x)以离散数据点形式给出：xix0x1…xnyi=f(xi)y0y1…yn定积分问题对一些典型的微分方程，如可别离变量方程、一阶线性方程等，有可能找出它们的一般解表达式，然后用初始条件确定表达式中的任意常数，这样即能确定解但是对于常微分方程初值问题：那么无法求出一般解常微分方程初值问题11.2多项式与插值来源于实际、又广泛用于实际。多项式插值的主要目的是用一个多项式拟合离散点上的函数值，使得可以用该多项式估计数据点之间的函数值。可导出数值积分方法，有限差分近似关注插值多项式的表达式、精度、选点效果。11.2.1关于多项式MATLAB命令一个多项式的幂级数形式可表示为：也可表为嵌套形式或因子形式

N阶多项式n个根，其中包含重根和复根。假设多项式所有系数均为实数，那么全部复根都将以共轭对的形式出现幂系数：在MATLAB里，多项式用行向量表示，其元素为多项式的系数，并从左至右按降幂排列。

例：

被表示为>>p=[2145]>>poly2sym(p)ans=2*x^3+x^2+4*x+5Roots:多项式的零点可用命令roots求的。

例：>>r=roots(p)得到

r=0.2500+1.5612i0.2500-1.5612i-1.0000所有零点由一个列向量给出。polyval:可用命令polyval计算多项式的值。例：计算y(2.5)

>>c=[3,-7,2,1,1]；xi=2.5;yi=polyval(c,xi)yi=23.8125如果xi是含有多个横坐标值的数组，那么yi也为与xi长度相同的向量。>>c=[3,-7,2,1,1];xi=[2.5,3];>>yi=polyval(c,xi)yi=23.812576.0000polyfit:给定n+1个点将可以唯一确定一个n阶多项式。利用命令polyfit可容易确定多项式的系数。

例：>>x=[1.1,2.3,3.9,5.1];>>y=[3.887,4.276,4.651,2.117];>>a=polyfit(x,y,length(x)-1)a=-0.20151.4385-2.74775.4370>>poly2sym(a)ans=-403/2000*x^3+2877/2000*x^2-27477/10000*x+5437/1000多项式为Polyfit的第三个参数是多项式的阶数。m阶多项式与n阶多项式的乘积是d=m+n阶的多项式:计算系数的MATLAB命令是:c=conv(a,b)多项式除多项式的除法满足:其中是商,是除法的余数。多项式和可由命令deconv算出。例：[q,r]=deconv(a,b)例

>>a=[2,-5,6,-1,9];b=[3,-90,-18];>>c=conv(a,b)c=6-195432-4539-792-162>>[q,r]=deconv(c,b)q=2-56-19r=0000000>>poly2sym(c)ans=6*x^6-195*x^5+432*x^4-453*x^3+9*x^2-792*x-162

分段插值：通过插值点用折线或低次曲线连接起来逼近原曲线。MATLAB实现可调用内部函数。命令1interp1功能：一维数据插值〔表格查找〕。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。格式1yi=interp1(x,Y,xi)%返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。参量x指定数据Y的点。假设Y为一矩阵，那么按Y的每列计算。11.2.2一维插值2024/3/14Method可取：‘nearest’

：最邻近插值；‘spline’

：三次样条插值；‘cubic’

：立方插值;缺省时：分段线性插值。注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。1一维插值一维插值函数：yi=interp1(x，y，xi，'method')插值方法被插值点插值节点xi处的插值结果2024/3/14

例：在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。用MATLAB作分段线性插值计算2024/3/14hours=1:12;temps=[589152529313022252724];h=1:.1:12;t=interp1(hours,temps,h)plot(hours,temps,'+',h,t)title('线性插值下的温度曲线'),xlabel('Hour'),ylabel('DegreesCelsius')用MATLAB作分段线性插值计算2024/3/14程序运行结果：用MATLAB作分段线性插值计算2024/3/14返回

2024/3/14Method可取：‘nearest’

：最邻近插值；‘spline’

：三次样条插值；‘cubic’

：立方插值;缺省时：分段线性插值。注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。用MATLAB作插值计算一维插值函数：yi=interp1(x，y，xi，'method')插值方法被插值点插值节点xi处的插值结果2024/3/14例:在上取11个点，作三次样条插值,观察三次样条插值曲线与g(x)的误差.用MATLAB作三次样条插值计算2024/3/14x0=linspace(-5,5,11);y0=1./(1+x0.^2);x=linspace(-5,5,100);y=interp1(x0,y0,x,'spline');x1=linspace(-5,5,100);y1=1./(1+x1.^2);plot(x1,y1,'k',x0,y0,'+',x,y,'r');M文件yangcha.m用MATLAB作三次样条插值计算2024/3/14返回

2024/3/14二维插值的提法mn个节点(xi,xj,zij)(i=1,2,…,m;j=1,2,…,n)其中xi,yj互不相同，不妨设a=x1<x2<…<xm=bc=y1<y2<…<yn=d求任一插值点(x*,y*)(≠(xi,yj))处的插值Z*第一种〔网格节点〕11.2.3二维插值2024/3/14二维插值的提法

xyx1x2x3x4x5y1y2y3

(x*,y*)2024/3/14第二种〔散乱节点〕n个节点其中互不相同，求任一插值点处的插值二维插值的提法2024/3/14二维插值方法的根本思路

构造一个二元函数通过全部已知节点,即再用计算插值，即或返回2024/3/14Method可取：

‘nearest’

最邻近插值；‘linear’

双线性插值；

‘cubic’

双三次插值；缺省时,双线性插值。二维插值:已有程序z=interp2(x0,y0,z0,x,y,’method’)被插值点插值方法插值节点被插值点的函数值网格节点数据的插值

注意：x0,y0为向量，但z0是矩阵，其列数等于x0的长度，行数等于y0的长度。例某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点0:2.5:10(米)，用h表示测量时间0:30:60(秒)，用T表示测试所得各点的温度(℃)。试用线性插值求出在一分钟内每隔20秒、钢轨每隔1米处的温度TI。命令如下：x=0:2.5:10;h=0:30:60;T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];xi=[0:10];hi=[0:20:60]';TI=interp2(x,h,T,xi,hi)2024/3/14用MATLAB作散点数据的插值计算

注意：x0,y0,z0均为向量，长度相等。Method可取‘nearest’,’linear’,’cubic’;‘linear’是缺省值。z=griddata(x0,y0,z0,x,y,’method’)被插值点插值节点被插值点的函数值x=rand(100,1)*4-2;y=rand(100,1)*4-2;z=x.*exp(-x.^2-y.^2);ti=-2:.25:2;[XI,YI]=meshgrid(ti,ti);ZI=griddata(x,y,z,XI,YI);mesh(XI,YI,ZI),holdplot3(x,y,z,'o'),holdoff11.3曲线拟合一组〔二维〕数据，即平面上n个点〔xi,yi)i=1,…n,寻求一个函数〔曲线〕y=f(x),使f(x)在某种准那么下与所有数据点最为接近，即曲线拟合得最好。+++++++++xyy=f(x)(xi,yi)

ii为点〔xi,yi)与曲线y=f(x)的距离问题的数学模型步骤：1)选定一类函数f(x，a1，a2，…，am)〔1〕其中a1,a2,…am为待定常数。2)确定参数a1,a2,…am准那么(最小二乘准那么〕：使n个点〔xi,yi)与曲线y=f(x，a1，a2，…，am)的距离i的平方和最小。曲线拟合的根本原理记

问题归结为:求a1,a2,…am使J(a1,a2,…am)最小。这样的拟合称为最小二乘拟合。问题的数学模型曲线拟合的根本原理最小二乘拟合函数的选取+++++++++++++++f=a1+a2x将数据(xi,yi)i=1,…n作图，通过直观判断确定f：f=a1+a2x+a3x2f=a1+a2x+a3x22.通过机理分析建立数学模型来确定f。+++++++++++++++f=a1+a2/xf=a1exp(a2x)f=a1exp(a2x)最小二乘拟合函数的选取多项式拟合：作多项式f(x)=a1xm+…+amx+am+1函数拟合,可利用已有程序polyfit,其调用格式为:a=polyfit(x,y,m)用MATLAB作最小二乘拟合数据点系数拟合多项式次数例.热敏电阻数据：温度t(0C)20.532.751.073.095.7电阻R()7658268739421032求电阻R随温度t的变化规律。1.多项式拟合用MATLAB作最小二乘拟合2)用命令polyfit(x,y,m)作最小二乘拟合3)编写MATLAB程序dianzu.m，并运行得到：a1=3.3940,a2=702.49181)选取拟合函数R=a1t+a21.多项式拟合用MATLAB作最小二乘拟合t=[20.532.5517395.7];r=[7658268739421032];aa=polyfit(t,r,1);a=aa(1)b=aa(2)y=polyval(aa,t);plot(t,r,'k+',t,y,'r')M文件dianzu.m1.多项式拟合用MATLAB作最小二乘拟合2.曲线拟合：

作一般的最小二乘曲线拟合，可利用已有程序lsqcurvefit,其调用格式为：

[a,resnorm,residual]=lsqcurvefit(‘f’,a0,x,y)

待定常数函数M文件拟合函数y=f(a,x)的函数M—文件各数据点处的误差误差平方和a的初值数据点用MATLAB作最小二乘拟合xdata=[0:0.1:2];ydata=[5.89553.56392.51731.97901.89901.39381.13591.00961.03430.84350.68560.61000.53920.39460.39030.54740.34590.13700.22110.17040.2636]例：用函数y(x)=a1*exp(-a2*x)+a3*exp(-a4*x)拟合以下数据点：1.用命令lsqcurvefit(fun,a0,x,y)2.曲线拟合用MATLAB作最小二乘拟合f=@(a,x)a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x);a0=[1,1,1,0];[a,resnorm,residual,flag,output]=lsqcurvefit(f,a0,xdata,ydata)xi=linspace(0,2,200);yi=f(a,xi);plot(xdata,ydata,'ro',xi,yi)xlabel('x'),ylabel('y=f(x)'),title('nonlinearcurvefitting')2.曲线拟合用MATLAB作最小二乘拟合1、方程(组)，f1(x)=0,…,fn(x)=0,x=(x1,…,xn)求解方程(组)的MATLAB语句2、方程(组)，f1(x)=0,…,fn(x)=0,x=(x1,…,xn)fun.mfunctionf=fun(x)f(1)=f1(x);……f(n)=fn(x)初值1）可以省略。2）options=1,表示输出中间结果。solve('f1(x)’,'f2(x)’,…,'fn(x)

’)X=fsolve(‘fun’,X0,options)11.4解方程组3、多项式方程：amxm+am-1xm-1+…+a0=0p=[am,am-1,…,a0]；roots(p)4、线性方程组：AX=b

其中A是m×n阶矩阵，b是m维向量。x=A\borx=inv(A)*b特点：可以找出全部根。特点：只能求出一个特解。输出:[1/2/a*(-b+(b^2-4*a*c)^(1/2))][1/2/a*(-b-(b^2-4*a*c)^(1/2))]①单变量方程solve()语句的用法例1：求解方程ax2+bx+c=0输入: x=solve('a*x^2+b*x+c')或solve('a*x^2+b*x+c=0')1〕符号解例2：解方程：x3-2x2=x-1解：

s=solve('x^3-2*x^2=x-1')double(s)2〕数值解solve()语句的用法s=7/(9*(108^(1/2)*((7*i)/108)+7/54)^(1/3))+((108^(1/2)*7*i)/108+7/54)^(1/3)+2/32/3-((108^(1/2)*7*i)/108+7/54)^(1/3)/2+(3^(1/2)*(7/(9*(108^(1/2)*((7*i)/108)+7/54)^(1/3))-((108^(1/2)*7*i)/108+7/54)^(1/3))*i)/2-7/(18*(108^(1/2)*((7*i)/108)+7/54)^(1/3))2/3-((108^(1/2)*7*i)/108+7/54)^(1/3)/2-(3^(1/2)*(7/(9*(108^(1/2)*((7*i)/108)+7/54)^(1/3))-((108^(1/2)*7*i)/108+7/54)^(1/3))*i)/2-7/(18*(108^(1/2)*((7*i)/108)+7/54)^(1/3))ans=2.2470+0.0000i0.5550-0.0000i-0.8019-0.0000i输入：[x,y]=solve('x^2*y^2,x-(y/2)-b')输出：x=[0]，y=[-2*b][0]，[-2*b]〔符号解〕[b]，[0][b]，[0]v=[x,y]②多变量方程组例4solve()语句的用法如果有10个方程的系统，输入[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10]=solve(……)不但费时，而且非常笨拙。solve可以给出结构输出形式。S=solve('x^2*y^2-2*x-1=0','x^2-y^2-1=0')S=x:[8x1sym]y：[8x1sym]〔给出结构解〕输出具体的解：S.x,S.y,s1=[S.x(2),S.y(2)]；〔取出解空间的第二组解〕M=[S.x,S.y]；〔创立解矩阵〕例5solve()语句的用法例6：求解方程组fsolve()语句的用法解①输入：[x,y,z]=solve('sin(x)+y^2+log(z)-7=0','3*x+2^y-z^3+1=0','x+y+z-5=0','x','y','z')x=0.59905375664056731520568183824539y=2.3959314023778168490940003756591z=2.0050148409816158357003177860955输出：解②：1〕建立方程组的M-函数文件(nxxf.m)functioneq=nxxf(x)eq(1)=sin(x(1))+x(2)^2+log(x(3))-7;eq(2)=3*x(1)+2^x(2)-x(3)^3+1;eq(3)=x(1)+x(2)+x(3)-5;运行程序(test4.m)y=fsolve('nxxf',[1,1,1],1)3〕运行结果：OptimizationTerminatedSuccessfullyy=0.59902.39592.0050fsolve()语句的用法输出:

-1.2131-0.9017+0.5753i-0.9017-0.5753i-0.2694+0.9406i-0.2694-0.9406i0.4168+0.8419i0.4168-0.8419i0.8608+0.3344i0.8608-0.3344i例7：求解多项式方程x9+x8+1=0roots()语句的用法输入:p=[1,1,0,0,0,0,0,0,0,1];roots(p)例8:AX=b,A\b和inv(A)*b语句的用法解:输入: A=[123;456;789];b=[6;14;-3];x1=A\b,x2=inv(A)*b输出: 警告:系统的秩缺乏.解不唯一.特殊情形：1、微分方程的一般形式：F(x,y,y’,…,y(n))=0

隐式或y(n)=f(x,y,y’,…,y(n-1))

显式2、一阶微分方程组的一般形式：n阶1阶初始条件：y(x0)=y011.5解微分方程微分方程解的形式①解析解y=f(x)②数值解(xi,yi)③图形解xyo①简单的微分方程。②复杂、大型的微分方程。一阶微分方程:获取解析解的方法归类:①别离变量法；如dy/dx=x*y；②齐次方程的变换法；如dy/dx=f(y/x)③线性方程的常数变易法或公式法.……解析解MATLAB软件实现解析解dsolve('eqn1','eqn2',…,'c1',…,'var1',…)微分方程组初值条件变量组注意：①y'Dy，y''D2y②自变量名可以省略，默认变量名‘t’。例①输入：y=dsolve('Dy=1+y^2')y1=dsolve('Dy=1+y^2','y(0)=1','x')输出：y=tan(t-C1)〔通解，一簇曲线〕y1=tan(x+1/4*pi)〔特解，一条曲线〕例②常系数的二阶微分方程y=dsolve('D2y-2*Dy-3*y=0','x')y=dsolve('D2y-2*Dy-3*y=0','y(0)=1,Dy(0)=0','x')输入：x=dsolve('D2x-(1-x^2)*Dx+x=0','x(0)=3,Dx(0)=0')上述两例的计算结果怎样？由此得出什么结论？例③非常系数的二阶微分方程例③无解析表达式！x=dsolve('(Dx)^2+x^2=1','x(0)=0')例④非线性微分方程x=[sin(t)][-sin(t)]假设欲求解的某个数值解，如何求解？t=pi/2;eval(x)输入：[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y')[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y','x(0)=0,y(0)=1')例④输出：〔li3.m〕数值解1、欧拉法2、龙格—库塔法数值求解思想：(变量离散化)

研究常微分方程的数值解法是十分必要的。[t，x]=solver〔’f’,ts,x0,options〕ode23

ode45ode113ode15sode23s由待解方程写成的m-函数文件ts=[t0，tf]，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔算法ode45：运用组合的4/5阶龙格-库塔算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分别为设定的相对误差和绝对误差.Matlab软件计算数值解f=@(x,y)-y+x+1;[x,y]=ode23(f,[0,1],1)plot(x,y,'r');holdonezplot('x+exp(-x)',[0,1])%精确解例1

y’=-y+x+1，y(0)=1标准形式：y’=f(x,y)1、在解n个未知函数的方程组时，x0和x均为n维向量，m-函数文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意:注意1：1、建立M文件函数functionxdot=fun(t,x,y)xdot=[f1(t,x(t),y(t))；f2(t,x(t),y(t))]；2、数值计算〔执行以下命令〕[t,x,y]=ode23(‘fun',[t0,tf],[x0,y0])注意：执行命令不能写在M函数文件中。xd(1)=f1(t,x(t),y(t));xd(2)=f2(t,x(t),y(t));xdot=xd’;%列向量例如：令注意2：functionxdot=fun1(t,x,y)〔fun1.m)xdot=[f(t,x(t),y(t))；x(t)]；[t,x,y]=ode23(‘fun1',[t0,tf],[x0,y0])M-文件函数如何写呢？注意：y(t)是原方程的解。x(t)只是中间变量。如果方程形式是：z’’’=f(t,z,z’’)?例2Vanderpol方程：令y1=x(t),y2=x’(t);

该方程是否有解析解？编写程序如下:

f=@(t,y)[y(2);(1-y(1)^2)*y(2)-y(1)];[t,y]=ode23(f,[0,20],[3,0]);y1=y(:,1);%原方程的解y2=y(:,2);plot(t,y1,t,y2,'--')%y1(t),y2(t)曲线图

pause,plot(y1,y2),%相轨迹图，即y2(y1)曲线grid,蓝色曲线——y〔1〕；〔原方程解〕

红色曲线——y〔2〕；计算结果例3考虑Lorenz模型：其中参数β=8/3，σ=10，ρ=28解：1〕编写匿名函数；2)数值求解并画三维空间的相平面轨线。f=@(t,x)[-8/3,0,x(2);0,-10,10;-x(2),28,-1]*x;x0=[000.1]';[t,x]=ode45(f,[0,10],x0);pl