计算机应用基础积分方程及应用

计算机应用基础积分方程及应用第一页，共四十页，编辑于2023年，星期五一、插值基础

设函数y=f(x)在区间[a,b]上有定义，且已知在点a≤x0<x1<…<xn≤b上的值y0,y1,…,yn,若存在一简单函数P(x),使

P（xi）=yi(i=0,1,2,…,n)（4-1）

成立，就称P(x)为f(x)的插值函数，点x0，x1，…，xn称为插值节点，包含插值节点的区间[a,b]称为插值区间，求插值函数P(x)的方法称为插值法。

若P(x)为插值多项式，即

P(x)=a0+a1x+…+anxn

（4-2）

其中ai为实数，就称P(x)为插值多项式，相应的插值法称为多项式插值。若P(x)为分段多项式，就称分段插值。

4.1插值函数第二页，共四十页，编辑于2023年，星期五

从几何图形上看，插值法就是求曲线y=P(x),使其通过给定的n+1个点（xi,yi),i=0,1,…,n,并用它近似已知曲线y=f(x)。

插值多项式存在的唯一性

设P(x)是形如（4-2）的插值多项式，则满足条件（4-1）的插值多项式（4-2）存在唯一。4.1插值函数第三页，共四十页，编辑于2023年，星期五一线性插值与抛物插值

通过平面的两点可以确定一条直线经过两点，这就是拉格朗日线性插值；对于不在一条直线的3个点得到的插值多项式称为抛物线插值。二、拉格朗日插值多项式插值多项式将上式展开，就得到插值的基函数，插值多项式为4.1插值函数第四页，共四十页，编辑于2023年，星期五拉格朗日插值函数通用程序functionf=Lang(x,y,x0)symst;if(length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等！');return;end%检错f=0;for(i=1:n)l=y(i);for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end;%x为已知数据点的x坐标向量；%y为已知数据点的y坐标向量；%x0为插值点的x坐标；%f为求得的拉格朗日插值多项式；%f0为在x0处的插值4.1插值函数第五页，共四十页，编辑于2023年，星期五for(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));%计算拉格朗日基函数

end;f=f+l;%计算拉格朗日插值函数

simplify(f);%化简

if(i==n)if(nargin==3)f=subs(f,'t',x0);%计算插值点的函数值

elsef=collect(f);%将插值多项式展开

f=vpa(f,6);%将插值多项式的系数化成6位精度的小数

endendend4.1插值函数第六页，共四十页，编辑于2023年，星期五三Matlab自带函数插值

3.1一维插值

Matlab提供了interp1函数用于一维插值，其调用格式为yi=interp1(x,Y,xi,method)对节点（x,Y)进行插值，计算插值点xi的函数值；Method插值算法，默认的为线性插值：

nearest:线性最近插值

linear:线性插值（默认）

spline:三次样条插值

pchip:分段三次埃尔米特插值

cubic:双三次插值4.1插值函数第七页，共四十页，编辑于2023年，星期五例【4-1】已知函数f(x)=(x2-12x+9)e-4xcosx生成数据点，试根据生成的数据进行插值处理。比较不同算法，计算程序“ME_7_1.m”4.1插值函数例【4-2】：某观测站测得某日6:00时至18:00时之间每隔2小时的室内外温度(℃)，用3次样条插值分别求得该日室内外6:30至17:30时之间每隔2小时各点的近似温度(℃)。第八页，共四十页，编辑于2023年，星期五设时间变量h为一行向量，温度变量t为一个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。命令如下：h=6:2:18;t=[18,20,22,25,30,28,24;15,19,24,28,34,32,30]';XI=6.5:2:17.5YI=interp1(h,t,XI,‘spline’)%用3次样条插值计算4.1插值函数第九页，共四十页，编辑于2023年，星期五3.2二维插值

Matlab提供interp2用于二维插值，适用于二维网格数据插值，其调用格式如下：ZI=interp2(X,Y,Z,XI,YI,method)X,Y均为矩阵在2-D区域的点，这些点的数值Z已知，因此构造插值函数Z=F(X,Y)4.1插值函数第十页，共四十页，编辑于2023年，星期五3.3二维随机点的插值对于任意多组数据（xi,yi,zi)采用griddata函数，其调用格式：[XI,YI,ZI]=griddata(x,y,z,XI,YI,method)x,y,z为已知样本点，[XI,YI]插值点的位置【例4-3】在区间【-2,2】上创建100个随机点，用griddata函数对这些数据点进行插值计算程序“ME_7_2.m”4.1插值函数第十一页，共四十页，编辑于2023年，星期五4.1插值函数3.4三次样条插值对于插值基点很多时，采用多项式插值效果并不是很好，采用三次样条插值效果更好一些，其调用格式：yy=spline(x,Y,xx),计算三次样条插值的分段多项式，采用ppval(pp,x)计算多项式在x点的值。pp=spline(x,Y)例4-4利用spline函数对y=tan(pi*x/24)进行样条插值。ME_7_6.m第十二页，共四十页，编辑于2023年，星期五functionME_7_6clear;x=0:10;y=tan(pi*x/24);xi=linspace(0,10);yi=spline(x,y,xi);plot(x,y,'rp',xi,yi);title('插值效果图');4.1插值函数第十三页，共四十页，编辑于2023年，星期五3.5B样条插值是一类常用的样条函数，其调用格式为：Spline=spapi(knots,x,y)spapi(k,x,y),k为B样条阶数4.1插值函数第十四页，共四十页，编辑于2023年，星期五一最小二乘曲线拟合4.2拟合函数有一组离散数据xi,yi,i=1,2…N,满足某一函数的原型yˇ(x)=f(a,x),其中a为待定系数向量，则最小二乘曲线拟合的目标：求出这一组待定系数的值，使得目标函数为最小。Matlab调用格式如下：第十五页，共四十页，编辑于2023年，星期五4.2拟合函数其中Fun为原型函数的Matlab表示，可以使M函数inline()函数；a0

为最优化的初值；

x,y为原始输入和输出数据向量；

a为返回的待定系数向量；

Jm为在此待定系数下目标函数的值。第十六页，共四十页，编辑于2023年，星期五4.2拟合函数

二最小二乘多项式拟合对于离散型函数，若数据点较多，若将每个数据点都当做插值节点，运算显得非常复杂。在工程试验中，常测得一组离散数据点（xi,yi),(i=1,2…N),要求y=(x），这种应变量只有一个自变量的数据拟合方法称之为直线拟合。(仍然采用最小二乘方法)p=polyfit(x,y,n)其中x和y为原始的样本点构成的向量n为选定的多项式阶数p为多项式系数按降幂排列得出的行向量第十七页，共四十页，编辑于2023年，星期五4.2拟合函数函数polyval()

根据由polyfit（）计算得到的多项式系数向量，计算在x处多项式函数的值y，它通常与polyfit（）联合使用，调用格式：y=polyval(p,x)p:多项式函数；x:需要计算点的x值。第十八页，共四十页，编辑于2023年，星期五4.2拟合函数三两种拟合方法的比较最小二乘曲线拟合最小二乘多项式拟合p=polyfit(x,y,n)拟合的函数是固定的,且是降阶的多项式,因此拟合调用不需要注明拟合函数,不需要初值。相同之处:都是采用最小二乘法优化获得拟合曲线系数拟合的函数形式可任意,因此拟合调用需要注明拟合函数,即需要建立一个Fun的函数，需要初值，而且结果与初值密切相关。第十九页，共四十页，编辑于2023年，星期五【例4-5】进行多项式拟合x0=0:0.1:1;y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);p3=polyfit(x0,y0,3)vpa(poly2sym(p3),10)x=0:0.01:1;%绘图ya=(x.^2-3*x+5).*exp(-5*x).*sin(x);y1=polyval(p3,x);plot(x,y1,x,ya,x0,y0,‘o’)拟合程序“ME_7_3”结果如下：4.2拟合函数

已知的数据点来自f(x)=(x2-3x+5)e-5xsinx第二十页，共四十页，编辑于2023年，星期五p3=2.8400-4.78981.94320.0598f(x)=2.84x3-4.79x2+1.94x+0.05984.2拟合函数第二十一页，共四十页，编辑于2023年，星期五例【4-6】4.2拟合函数已知数据可能满足y(x)=ax+bx2e-cx+dXi0.10.20.30.40.5Yi2.322.652.973.293.6Xi0.60.70.80.91Yi3.914.214.524.825.13求满足最小二乘解a,b,c,d的值令a1=a,a2=b,a3=c,a4=d,则原型函数可写成：

y(x)=a1x+a2x2e-a3x+a4第二十二页，共四十页，编辑于2023年，星期五functionME_7_4clear;x=0.1:0.1:1;y=[2.32,2.65,2.97,3.29,3.6,...3.91,4.21,4.52,4.82,5.13];a=lsqcurvefit(@c8f3,[1;2;2;3],x,y);a'y1=c8f3(a,x);plot(x,y,x,y1)计算程序4.2拟合函数functiony=c8f3(a,x)y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);编写原型函数程序第二十三页，共四十页，编辑于2023年，星期五4.2拟合函数四最小二乘拟合生成样条曲线函数csaps()的用法:

平滑生成三次样条曲线,调用格式:样条拟合曲线不要求曲线通过全部实验数据点，更适用实验观察得到的离散数据点sp=csaps(x,y,p)ys=casps(x,y,p,xx,w)函数span2()的用法:

用最小二乘法生成B样条曲线,调用格式:sp=span2(knots,k,x,y)sp=span2(knots,k,x,y,w)第二十四页，共四十页，编辑于2023年，星期五4.2拟合函数函数spaps()的用法:

平滑生成B样条曲线,调用格式:sp=spaps(x,y,tol)[sp,ys]=spaps(x,y,tol,m,w)函数fnval()的用法:计算函数f在给定点x处的函数值，调用格式:values=finval(f,x)第二十五页，共四十页，编辑于2023年，星期五t/h青霉素浓度t/h青霉素浓度001209430201061401095040160016010280603000180962080581020094001008600例4-6青霉素发酵实验数据如表所示，请分别用插值方法和最小二乘法估算t=10,30,50,70,90,110,150,190h时青霉素的浓度。并分别画出插值函数和拟合函数曲线，并标注插值点。（ME_7_5.m)4.2拟合函数第二十六页，共四十页，编辑于2023年，星期五4.3数值积分及应用一复合辛普森法积分二梯形求积三单变量数值积分问题求解四双重积分问题的数值解第二十七页，共四十页，编辑于2023年，星期五4.3数值积分问题对于一个数学积分关系式：

常用的解法有：矩形法、梯形法和辛普森法。

一复合辛普森法

基本思想：将函数f(x)的积分区间n等分（n为偶数），以两个小区间作为一个计算单元，用一个简单的二次函数近似被积函数，则函数f(x)在区间[a,b]积分可表述为：

S=h/3{f(a)-f(b)+2∑[2f(a+(2i-1)h)+f(a+2ih)]}

第二十八页，共四十页，编辑于2023年，星期五离散误差：4.3数值积分问题调用格式：q=quad(@fun,a,b,tol)fun:被积函数，必须是向量表达式，表达式必

须用点运算符。a,b：积分限tol：绝对误差限函数quad1()精度更高，使用方法同前。第二十九页，共四十页，编辑于2023年，星期五例4-4用复合simposn法计算以下积分，要求精度在1e-5内，并计算误差。第三十页，共四十页，编辑于2023年，星期五要使则M>=5,取m=5,h=0.2程序ME_3_2.m,S=3.8875,err=2.72864.3数值积分问题第三十一页，共四十页，编辑于2023年，星期五二由给定数据进行梯形求积4.3数值积分问题第三十二页，共四十页，编辑于2023年，星期五4.3数值积分问题利用trapz()函数进行积分求解,格式如下:S=trapz(x,y)例4-7用定步长方法求解积分画图：>>x=[0:0.01:3*pi/2,3*pi/2];y=cos(15*x);plot(x,y)计算程序：ME_4_9第三十三页，共四十页，编辑于2023年，星期五求理论值：不同步距：>>symsx,A=int(cos(15*x),0,3*pi/2)h0=[0.1,0.01,0.001,0.0001,0.00001,0.000001];v=[];Forh=h0x=[0:h:3*pi/2,3*pi/2];y=cos(15*x);I=trapz(x,y);v=[v;h,I,1/15-I];end4.3数值积分问题

>>symsx,A=int(cos(15*x),0,3*pi/2)第三十四页，共四十页，编辑于2023年，星期五二双重积分问题的数值解4.3数值积分问题求解格式:(1)y=dblquad(Fun,xm,xM,ym,yM)%矩形区域的双重积分(2)y=dblquad(Fun,xm,xM,ym,yM,)%限定精度的双重积分