Matlab 科学计算

1Matlab与科学计算计算机学院刘咏梅Email:[e-mail]2第二章MATLAB数值计算基础•IntroductiontoMATLAB•MATLABBasicsofNumericalComputing•MATLABProgramming(branchingandLoops)•MATLABProgramming(PredefinedFunctions)•MATLABGraphicandImageProcessing•ScientificComputingusingMATLAB3本章学习的目标：•利用MATLAB进行插值和拟合•利用MATLAB求积分和微分•利用MATLAB求解线性方程组4插值和拟合•插值：利用函数f(x)在某区间中若干个点的函数值，做适当特定的函数f’(x)，使区间上其它点的值用f’(x)作为f(x)的近似。•拟合：已知某函数f(x)若干离散值，通过调整f(x)中若干待定系数，使f(x)与已知点集的差别最小。5插值和拟合示意图插值示意图拟合示意图6插值与拟合的应用背景•在生产实际和科学研究中，会遇到大量函数，其中相当一部分是通过测量和实验观测得到的，虽然其函数关系y=f(x)在某区间上客观存在，但却不知道具体的解析表达式。只能得到函数在区间[a,b]上的一些离散点的函数值、导数值等。•我们希望对这样的函数用一个比较简单的函数的表达式来近似地给出其整体描述。有时是虽然有明确表达式，但函数本身过于复杂而不便于进行数值计算，同样希望构造一个既能反映函数性质又便于计算的简单函数，来近似代替原来的函数，插值就是寻求近似函数的方法之一。7计算方法中的插值法•多项式插值：设函数y=f(x)在区间[a,b]上有定义，已知[a,b]区间上n+1互异点x0,x1,…,xn及其函数值y0,y1,…,yn,若存在一个简单函数y=p(x),使其经过y=f(x)上的这n+1个已知点(x0,y0),(x1,y1),…,(xn,yn),那么函数p(x)称插值函数，点x0,x1,…,xn称为插值节点，点(x0,y0),(x1,y1),…,(xn,yn)称为插值点，包含插值点的区间称为插值区间，求p(x)的方法叫做插值法，f(x)称为被插函数。若p(x)表示为pn(x)=a0+a1x+a2x2+…+anxn则称pn(x)为n次插值多项式，相应的插值法称为多项式插值。8计算方法中的插值法•Lagrange插值法•分段线性插值法9Lagrange插值法•因为pn(x)=a0+a1x+a2x2+…+anxn且pn(xi)=yi，i=0,1,2,…,n，插值问题实际上就是根据给定的n+1个节点上的函数值，确定n次多项式pn(x)的n+1个系数，即由n+1个条件确定n+1个待定系数。也就是求解线性方程组的解。•但是系数矩阵为Vandermonde矩阵，条件数常常很大，是病态问题，小的扰动就会产生大的偏移量，求精确解失去意义。10Lagrange插值法•Lagrange插值公式如下：•对任意k(k=0,…,n),Lagrange插值基函数：nknkjjjkjkxxxxyxy00)()())(())(()())(())(()(11101110nkkkkkkknkkkxxxxxxxxxxxxxxxxxxxxxl11分段线性插值法•高次插值多项式可能产生极大的误差（这一现象称为荣格（Runge）现象）因此不宜用太多的点来做插值多项式。•为了解决荣格现象，引入了分段线性插值，即通过n+1个插值点用折线段连接起来逼近原曲线，这也是计算机绘制图形的基本原理。12荣格（Runge）现象x=[-5:1:5];y=1./(1+x.^2);x0=[-5:0.1:5];y0=lagrange(x,y,x0);y1=1./(1+x0.^2);plot(x0,y0,'--r');holdonplot(x0,y1,'-b');13一维插值函数interp1•YI=interp1(X,Y,XI,’method’)interpolatestofindYI.•ThevaluesoftheunderlyingfunctionYatthepointsinthearrayXI.•XmustbeavectoroflengthN.IfYisavector,thenitmustalsohavelengthN,andYIisthesamesizeasXI.•Method=nearest/linear(default)/spline/pchip/cubic14插值实例x=0:10;y=sin(x);xi=0:.25:10;yi=interp1(x,y,xi);plot(x,y,'o',xi,yi)15插值实例x=[-5:1:5];y=1./(1+x.^2);x0=[-5:0.1:5];y0=interp1(x,y,x0);plot(x0,y0,'*m')holdonplot(x,y,‘bo')16二维插值函数interp2•ZI=interp2(X,Y,Z,XI,YI,’method’)interpolatestofindZI.•Thevaluesoftheunderlying2-DfunctionZatthepointsinmatricesXIandYI.•MatricesXandYspecifythepointsatwhichthedataZisgiven.•Method=nearest/linear(default)/spline/cubic17拟合•在实际中，往往需要从一组实验数据(xi,yi),中找到自变量和因变量的函数关系y=f(x)。由于观测数据往往不够准确，因此不要求y=f(x)经过所有的点(xi,yi),而只要求在给定点xi上的误差det(i)=f(xi)-yi按某种要求达到最小。•最小二乘拟合采用的策略是使误差的平方和达到最小。18Matlab拟合函数•利用polyfit函数进行多项式拟合•格式：p=polyfit(X,Y,N)p(1)*X^N+p(2)*X^(N-1)+…+p(N)*X+p(N+1)•实例：x=[0.51.01.52.02.53.0];y=[1.752.453.814.808.008.60];a=polyfit(x,y,2)a=0.56140.82871.156019x1=[0.5:0.05:3.0];y1=a(1)*x1.^2+a(2)*x1+a(3);plot(x,y,’*’)holdonplot(x1,y1,’-r’)20计算方法中的求积公式•微积分是理工科大学生踏入大学校门后首先要接触到的一门重要基础课程。•实践证明，微积分运算也是实际工作中广泛应用的基本工具，它是科学计算最重要的组成部分21计算方法中的求积公式•Newton-Cotes系列数值求积公式梯形公式复合梯形公式22Newton-Cotes系列数值求积公式•一种等间隔闭式求积公式•将积分区间[a,b]划为n等分，步长h=(b-a)/n，选取等距离节点xk=a+kh，构造下面的插值型求积公式：nkknkbanxfCabdxxfI0)()(Cotes系数23Cotes系数当n=1时，c0=1/2,c1=1/2当n=2时，c0=1/6,c1=2/3,c2=1/6当n=3时，c0=1/8,c1=3/8,c2=3/8,c3=1/8………梯形公式Simpson公式24梯形公式))()((2afbfabdxxfIban25复合梯形求积公式•复合梯形求积公式的由来：为了提高精度，常把区间[a,b]分成若干个小区间，在每个小区间上使用次数不高的Newton-Cotes求积公式，如梯形公式和Simpson公式，再将所有的子区间上的积分值累加起来。•复合梯形求积公式：f(x)在区间[a,b]上的积分值为])()(2)([211nkkbfxfafnab26MATLAB求积函数trapz•trapz:trapezoidalnumericalintegration(梯形数值积分)•Form：Z=TRAPZ(Y)computesanapproximationoftheintegralofYviathetrapezoidalmethod(withunitspacing).Tocomputetheintegralforspacingdifferentfromone,multiplyZbythespacingincrement.•Comments:Forvectors,TRAPZ(Y)istheintegralofY.Formatrices,TRAPZ(Y)isarowvectorwiththeintegralovereachcolumn.27trapz实例•Z=TRAPZ(X,Y,DIM)orTRAPZ(Y,DIM)integratesacrossdimensionDIMofY.•Example:IfY=[012•345]•thentrapz(Y,1)is[1.52.53.5]andtrapz(Y,2)is[28];28trapz实例x=[0:2*pi/100:2*pi];y=sin(x);z=trapz(y)z=0x=[0:pi/100:pi];y=sin(x);z=trapz(y)*pi/100z=1.999829思考。。。•如何利用trapz求下列积分的值305.0)6/sin(dttet30解线性方程组的应用背景•求解线性方程组在科学计算中最常遇到，数值方法如三次样条、有限元法等都要涉及到求解线性方程组。•求解线性方程组的实际应用也是非常广泛：天气预报、石油勘探、地质测绘、电路分析、分子结构预测等。31解线性方程组的理论方法•克莱姆法则（Cramer’sRule）：若线性方程组⑴的系数行列式D≠0，则线性方程组⑴有唯一解，其解为xj=Dj/D，其中Dj是把D中第j列元素对应地换成常数项，j=1,…,n,而其余各列保持不变所得到的行列式。32解线性方程组的理论方法•如果按克莱姆法则求解，假设计算机每秒可计算1亿次乘法，共n+1个行列式，每个行列式计算需要n!(n+1)次乘法，故共需(n2-1)n!次乘法。当问题规模n=10时的计算量耗时3.59秒，而n=100时已经非常大了，不能接受。•因此说，理论上比较完善的方法在实际中却用处不大。33利用计算机来解线性方程组的方法•解线性方程组的数值方法有两类：直接法和迭代法。•直接法是指通过有限次的四则运算求解方程组的解的方法。但由于舍入误差的存在，此类方法只能得到近似解，只能适用于低阶方程。•迭代法是先给出一个初始近似值，然后按一定的法则，逐步求得更加精确的近似解的方法。高阶方程组的数值解常用迭代法。•上述两种方法均包含许多种算法，有的解法是针对一般性方程组，有的则是针对一些有特点的方程组。34MATLAB求解线性方程组•AX=B的解可以用X＝A\B求解，XA=B的解可以用X=A/B求解。如果A是m×n的矩阵，当m＝n时可以得到唯一解。•在MATLAB中，只需用左除或右除就可以解得，虽然表面上是一个简单的符号，而它内部却包含着许多自适应算法，对不同类型的方程组用不同的方法。35高斯消去法•高斯消去法是通过一系列初等行变换来实现将增广矩阵简化为一个上三角阵，向量b在这一过程中也将被修改，再由回代过程解得向量x，类似于我们中学学过的消元法。36高斯消去法实例•用Gauss消去法解方程组解消元得37再由回代过程，得解x1=x2=x3=138第六章作业•第六章需要提交实验报告•实验目的掌握利用Matlab进行插值和拟合的方法掌握利用Matlab求积分和微分的方法掌握利用Matlab求解线性方程组的方法39实验（5）利用Matlab进行科学计算本次实验课，完成下面各项实验的其中之一：（1）编制函数文件，实现Lagrange插值法，然后f(x)=lnx