项目二 数值计算.docx
- 文档编号:9913806
- 上传时间:2023-02-07
- 格式:DOCX
- 页数:14
- 大小:73.01KB
项目二 数值计算.docx
《项目二 数值计算.docx》由会员分享,可在线阅读,更多相关《项目二 数值计算.docx(14页珍藏版)》请在冰豆网上搜索。
项目二数值计算
项目二数值计算
实验类型验证型实验学时8
实验目的学习Matlab软件的数值计算功能及操作
实验原理Matlab软件的数值计算基本功能及操作
实验内容
1.数据插值
插值法是实用的数值方法,是函数逼近的重要方法。
在生产和科学实验中,自变量x与因变量y的函数y=f(x)的关系式有时不能直接写出表达式,而只能得到函数在若干个点的函数值或导数值。
当要求知道观测点之外的函数值时,需要估计函数值在该点的值。
如何根据观测点的值,构造一个比较简单的函数y=φ(x),使函数在观测点的值等于已知的数值或导数值。
用简单函数y=φ(x)在点x处的值来估计未知函数y=f(x)在x点的值。
寻找这样的函数φ(x),办法是很多的。
φ(x)可以是一个代数多项式,或是三角多项式,也可以是有理分式;φ(x)可以是任意光滑(任意阶导数连续)的函数或是分段函数。
函数类的不同,自然地有不同的逼近效果。
在许多应用中,通常要用一个解析函数(一、二元函数)来描述观测数据。
根据测量数据的类型:
1.测量值是准确的,没有误差。
2.测量值与真实值有误差。
这时对应地有两种处理观测数据方法:
1.插值或曲线拟合。
2.回归分析(假定数据测量是精确时,一般用插值法,否则用曲线拟合)。
MATLAB中提供了众多的数据处理命令。
有插值命令,有拟合命令,有查表命令。
1.1插值命令
命令1interp1
功能一维数据插值(表格查找)。
该命令对数据点之间计算内插值。
它找出一元函数f(x)在中间点的数值。
格式yi=interp1(x,Y,xi)%返回插值向量yi,每一元素对应于参量xi,同时由向量x与Y的内插值决定。
参量x指定数据Y的点。
若Y为一矩阵,则按Y的每列计算。
yi是阶数为length(xi)*size(Y,2)的输出矩阵。
yi=interp1(Y,xi)%假定x=1:
N,其中N为向量Y的长度,或者为矩阵Y的行数。
yi=interp1(x,Y,xi,method)%用指定的算法计算插值:
’nearest’:
最近邻点插值,直接完成计算;
’linear’:
线性插值(缺省方式),直接完成计算;
’spline’:
三次样条函数插值。
对于该方法,命令interp1调用函数spline、ppval、mkpp、umkpp。
这些命令生成一系列用于分段多项式操作的函数。
命令spline用它们执行三次样条函数插值;
’pchip’:
分段三次Hermite插值。
对于该方法,命令interp1调用函数pchip,用于对向量x与y执行分段三次内插值。
该方法保留单调性与数据的外形;
’cubic’:
与’pchip’操作相同;
’v5cubic’:
在MATLAB5.0中的三次插值。
例1.1已知某产品从1900至2010每间隔10年的产量,试用不同的插值方法计算其在1995年的产量,并对不同的插值方法得到的结果进行比较。
year=1900:
10:
2010;
product=[75.99591.972105.711123.203131.669150.697179.323203.212226.505249.633256.344267.893];
p19951=interp1(year,product,1995,'pchip');%分段三次插值
p19952=interp1(year,product,1995,'nearest');%最邻近插值
p19953=interp1(year,product,1995,'linear');%线形插值
p19954=interp1(year,product,1995,'spline');%三次样条插值
x=1900:
1:
2010;
y1=interp1(year,product,x,'pchip');
y2=interp1(year,product,x,'nearest');
y3=interp1(year,product,x,'linear');
y4=interp1(year,product,x,'spline');
subplot(2,2,1);plot(year,product,'o',x,y1)
subplot(2,2,2);plot(year,product,'o',x,y2)
subplot(2,2,3);plot(year,product,'o',x,y3)
subplot(2,2,4);plot(year,product,'o',x,y4)
结果:
p19951=75.9950;p19952=75.9950;p19953=75.9950;p19954=75.9950
命令2interp2
功能二维数据内插值(表格查找)
格式ZI=interp2(X,Y,Z,XI,YI)%返回矩阵ZI,其元素包含对应于参量XI与YI(可以是向量、或同型矩阵)的元素,即Zi(i,j)←[Xi(i,j),yi(i,j)]。
用户可以输入行向量和列向量Xi与Yi,此时,输出向量Zi与矩阵meshgrid(xi,yi)是同型的。
同时取决于由输入矩阵X、Y与Z确定的二维函数Z=f(X,Y)。
参量X与Y必须是单调的,且相同的划分格式,就像由命令meshgrid生成的一样。
若Xi与Yi中有在X与Y范围之外的点,则相应地返回nan(NotaNumber)。
ZI=interp2(Z,XI,YI)%缺省地,X=1:
n、Y=1:
m,其中[m,n]=size(Z)。
再按第一种情形进行计算。
ZI=interp2(Z,n)%作n次递归计算,在Z的每两个元素之间插入它们的二维插值,这样,Z的阶数将不断增加。
interp2(Z)等价于interp2(z,1)。
ZI=interp2(X,Y,Z,XI,YI,method)%用指定的算法method计算二维插值:
’linear’:
双线性插值算法(缺省算法);
’nearest’:
最临近插值;
’spline’:
三次样条插值;
’cubic’:
双三次插值。
例2已知1950到1990每间隔10年,工龄10年,20年,30年的劳动报酬,试用不同的二维插值方法计算1975年,工龄15年的平均工资,并对结果进行比较。
>>years=1950:
10:
1990;
>>service=10:
10:
30;
>>wage=[150.697199.592187.625
179.323195.072250.287
203.212179.092322.767
226.505153.706426.730
249.633120.281598.243];
>>w=interp2(service,years,wage,15,1975)
2.数据拟合
问题:
给定一批数据点,需确定满足特定要求的曲线或曲面
若要求所求曲线(面)通过所给所有数据点,就是插值问题;
若不要求曲线(面)通过所有数据点,而是要求它反映对象整体的变化趋势,这就是数据拟合,又称曲线拟合或曲面拟合.
函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似,由于近似的要求不同,二者在数学方法上是完全不同的.
2.1多项式拟合:
命令:
p=polyfit(X,Y,n)
给出已知数据(X,Y)的n次多项式拟合,输出为多项式系数(按降幂排列)
命令:
polyval(p,x0)
利用上面拟合出的多项式p,计算在x0处的值,x0可以是单个点,也可以是数组。
例2.1求如下数据的拟合曲线:
X=0.5:
0.5:
3;
Y=[1.752.453.814.80,7.008.60];
P=polyfit(X,Y,2)
Y1=polyval(P,X);
f=inline('c
(1)*exp(c
(2)*x)','c','x');
c=lsqcurvefit(f,[0,0],X,Y);
sum((Y-Y1).^2)
sum((Y-f(c,X)).^2)
P=
0.56140.82871.1560
c=
1.46320.6003
0.17810.4043
2.2plot(X,Y,'+',X,Y1,X,f(c,X),'*')非线性拟合
命令:
c=lsqcurvefit(fun,c0,x,y)
说明:
fun是预先定义的一个函数文件(M文件),c0是初值。
例xdata=[-0.200.10.20.150.3];
ydata=[1.501.060.950.840.860.72];
f=inline('c
(1)*exp(c
(2)*x)','c','x');
c=lsqcurvefit(f,[0,0],xdata,ydata);
plot(xdata,ydata,'*',xdata,f(c,xdata))
sum((f(c,xdata)-ydata).^2)
例3根据美国人口从1790年到1990年间的人口数据(如下表),确定人口指数增长模型(Logistic模型)中的待定参数,估计出美国2010年的人口,同时画出拟合效果的图形。
表1美国人口统计数据
年份
1790
1800
1810
1820
1830
1840
1850
人口(×106)
3.9
5.3
7.2
9.6
12.9
17.1
23.2
年份
1860
1870
1880
1890
1900
1910
1920
人口(×106)
31.4
38.6
50.2
62.9
76.0
92.0
106.5
年份
1930
1940
1950
1960
1970
1980
人口(×106)
123.2
131.7
150.7
179.3
204.0
226.5
提示:
可采用以下两种模型
指数增长模型:
T=[1790:
10:
1980];
X=[3.95.37.29.612.917.123.231.438.650.262.976.092.0106.5123.2131.7150.7179.3204.0226.5];
f=inline('r
(1)*exp(r
(2)*T)','r','T');
r=lsqcurvefit(f,[0,0],T,X)
plot(T,X,'*',T,f(r,T))
norm(feval(f,r,T)-X)^2
Logistic模型:
T=[1790:
10:
1980];
X=[3.95.37.29.612.917.123.231.438.650.262.976.092.0106.5123.2131.7150.7179.3204.0226.5];
f=inline('r
(1)/(1+(r
(1)/r
(2)-1)*exp(-r(3)*T)','r','T');
r=lsqcurvefit(f,[227,3.9,0],T,X)
plot(T,X,'*',T,f(r,T))
norm(feval(f,r,T)-X)^2
可参考拟合函数:
a=lsqcurvefit('example_curvefit_fun',a0,x,y);
3.方程及方程组的求解
命令:
solve,fsolve
格式:
solve的基本功能是求解符合方程(组)SOLVE('eqn1','eqn2',...,'eqnN')
SOLVE('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')
X=FSOLVE(FUN,X0);
X=FSOLVE(FUN,X0,OPTIONS)
还可以使用null命令求解线性方程组AX=0的基础解系(解空间),格式为null(A);
还可以使用矩阵左除运算‘\’求解方程组AX=B,即使用命令:
X=A\B.
例1求方程组
的解。
解:
>>A=[56000
15600
01560
00156
00015];
B=[10001]';
R_A=rank(A)%求秩
X=A\B%求解
或用函数rref求解:
>>C=[A,B]%由系数矩阵和常数列构成增广矩阵C
>>R=rref(C)%将C化成行最简行
R=
1.000000002.2662
01.0000000-1.7218
001.0000001.0571
0001.00000-0.5940
00001.00000.3188
则R的最后一列元素就是所求之解。
例2求方程组
的一个特解。
解:
>>A=[11-3-1;3-1-34;15-9-8];
>>B=[140]';
>>X=A\B%由于系数矩阵不满秩,该解法可能存在误差。
X=[00-0.53330.6000]’(一个特解近似值)。
若用rref求解,则比较精确:
>>A=[11-3-1;3-1-34;15-9-8];
B=[140]';
>>C=[A,B];%构成增广矩阵
>>R=rref(C)
例3求解方程组的通解:
解:
>>A=[1221;21-2-2;1-1-4-3];
>>formatrat%指定有理式格式输出
>>B=null(A,'r')%求解空间的有理基
注意:
非齐次线性方程组需要先判断方程组是否有解,若有解,再去求通解。
因此,步骤为:
第一步:
判断AX=b是否有解,若有解则进行第二步
第二步:
求AX=b的一个特解
第三步:
求AX=0的通解
第四步:
AX=b的通解=AX=0的通解+AX=b的一个特解。
例4求解方程组
解:
在Matlab中建立M文件如下:
A=[1-23-1;3-15-3;212-2];
b=[123]';
B=[Ab];
n=4;
R_A=rank(A)
R_B=rank(B)
formatrat
ifR_A==R_B&R_A==n%判断有唯一解
X=A\b
elseifR_A==R_B&R_A X=A\b%求特解 C=null(A,'r')%求AX=0的基础解系 elseX='equitionnosolve'%判断无解 end 例5求解方程组的通解: 解法一: 在Matlab编辑器中建立M文件如下: A=[11-3-1;3-1-34;15-9-8]; b=[140]'; B=[Ab]; n=4; R_A=rank(A) R_B=rank(B) formatrat ifR_A==R_B&R_A==n X=A\b elseifR_A==R_B&R_A X=A\b C=null(A,'r') elseX='Equationhasnosolves' end 解法二: 用rref求解 A=[11-3-1;3-1-34;15-9-8]; b=[140]'; B=[Ab]; C=rref(B)%求增广矩阵的行最简形,可得最简同解方程组。 以上例子,请根据程序运行结果,自行写出其(通,特)解。 4微分方程及方程组的求解 要求解微分方程(组)dy/dt=f(t,y),可如下调用: [T,Y]=ode45(fun,[t0,tn],y0) 1)函数在求解区间[t0,tn]内,自动设立采样点向量T,并求出解函数y在采样点T处的样本值Y。 2)fun是一个函数文件,可以是自定义的M文件,也可以是符号表达式,如fun=inline(‘……’)要有两个参数,第一个参数是自变量t,第二个参数是因变量y。 3)y0=y(t0)给定方程的初值。 例: 求微分方程初值问题dy/dx=-2y/x+4x,y (1)=2在[1,3]区间内的数值解,并将结果与解析解进行比较。 先建立一个该函数的m文件fxy1.m: functionf=fxy1(x,y) f=-2.*y./x+4*x%注意使用点运算符 再输入命令: [X,Y]=ode45('fxy1',[1,3],2); X'%显示自变量的一组采样点 Y'%显示求解函数与采样点对应的一组数值解 (X.^2+1./X.^2)'%显示求解函数与采样点对应的一组解析解 Plot(X,Y) 例: 求解常微分方程组初值问题在区间[0,2]中的解。 建立一个函数文件fxy2.m: functionf=fxy2(x,y) f (1)=y (2); f (2)=-x.*y (2)+x.^2-5; f=f'; 在MATLAB命令窗口,输入命令: [X,Y]=ode45('fxy2',[0,2],[5,6])
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 项目二 数值计算 项目 数值 计算