第7章计算方法的MATLAB实现讲稿文档格式.docx
《第7章计算方法的MATLAB实现讲稿文档格式.docx》由会员分享，可在线阅读，更多相关《第7章计算方法的MATLAB实现讲稿文档格式.docx（21页珍藏版）》请在冰豆网上搜索。
7.1.2roots函数
用roots函数计算多项式的根，该函数的调用格式为：
●r=roots（c）：
返回一个列向量，其元素为多项式c的解。
例7－2 求多项式方程
P=[1-51-5];
%表示多项式
r=roots（P）%求多项式方程的解
P1=poly（r）%由多项式的根返回多项式的系数
r=
5.0000
-0.0000+1.0000i
-0.0000-1.0000i
P1=
1.0000-5.00001.0000-5.0000
注：
poly和roots互为逆函数，函数poly返回多项式的系数。
7.2非线性方程组的数值解法
用MATLAB求解线性方程组的解，在6.3中已介绍了基于矩阵变换的直接解法，除此方法之外，还有很多方法，比如，Jocabi迭代法、Gauss-Seidel迭代法和SOR（超松驰）迭代法等，这里不作详细介绍。
可参阅有关文献。
实际工程中得到的数学模型往往具有非线性的特点，得到其解析解比较困难。
一般采用迭代法求解非线性方程组。
比较常见的迭代方法有不动点迭代法、Newton迭代法和拟Newton迭代法等几种。
本节仅介绍不动点迭代法。
设含有
个未知数和
个方程的非线性方程组记为：
其中，
为由
个未知数构成的向量，
个函数构成的向量。
将方程组
改写为便于迭代的等价形式：
，由此得出不动点迭代法的迭代公式：
。
如果得到的序列
满足
，则
就是
的不动点。
这样就可以求出线性方程组的解（近似解）。
据此，编写不动点迭代法的M文件（文件名为example7_2）：
functions=example7_2（x,eps）
%用不动点迭代法求非线性方程组的解
%x为迭代初值，eps为允许误差值。
ifnargin==1
eps=1e-6;
elseifnargin<
1
error
return
end
x1=example7_2a（x）;
%example7_2a是函数
的文件名。
whilenorm（x1-x）>
=eps
x=x1;
x1=example7_2a（x）;
%循环迭代
s=x1;
return
例7－3 用不动点迭代法求方程组
解：
变形为
首先创建函数代码（文件名为example7_2a.m）：
functiony=example7_2a（x）
y
（1）=（x
（1）*x
（1）+x
（2）*x
（2）+9）/（-11）;
y
（2）=（x
（1）*x
（2）*x
（2）+x
（1）+8）/6;
再作图形，程序代码：
ezplot（'
x1^2+11*x1+x2^2+9'
[-12,6]）;
holdon
x1*x2^2+x1-6*x2+8'
holdoff
title（'
x1^2+11*x1+x2^2+9和x1*x2^2+x1-6*x2+8的图形'
）%修改标题
运行结果见图7－1：
图7－1
最后调用函数example7_2求方程组解：
x=example7_2（[0,0]）
x=
-1.00001.0000
这样就可以求得方程组
的一个解，从图7－1可以看出，方程组
还有一个解，为了求得它，还需要对原方程组作另外变形，这里不再作进一步的讨论。
7.3插值
插值计算在数据拟合和数据平滑等方面应用普遍。
MATLAB提供了用最近邻插值、线性插值、三次样条插值、三次插值和FFT插值法进行一维、二维、三维和高维插值的函数。
7.3.1一维插值
MATLAB中有两种一维插值，即多项式插值和基于FFT的插值。
⒈多项式插值
函数interp1进行一维插值。
一维插值是进行数据分析和曲线拟合的重要手段。
interp1函数使用多项式技术，用多项式函数拟合所提供的数据，并计算目标插值点上的插值函数值，其最常用的语法形式是：
yi=interp1（x,y,xi,method）
x和y为给定数据的向量，长度相同。
xi为包含要插值的点的构成的向量，method是一个可选的字符串，指定一种插值方法，包括：
●最近邻插值（method='
nearest'
）：
该方法将插值点的值设置为已知数据点中距离最近的点的值；
●线性插值（method='
linear'
该方法用线性函数拟合每对数据点，并返回xi处的相关函数值；
●三次样条插值（method='
spline'
该方法用三次样条函数拟合每对数据点，用spline函数在插值处进行三次样条插值；
●三次插值（method='
pchip'
或'
cubic'
该方法用pchip函数对向量x和y进行分段三次Hermite插值。
这几种方法在速度、内存和平滑性等方面有所差别，使用时可以根据需要进行选择，包括：
◣最近邻插值是最快的方法，但是，利用它得到的结果平滑性最差；
◣线性插值比最近邻插值要占用更多的内存，运行时间略长。
与最近邻法不同，它生成的结果是连续的，但是在顶点处有坡度变化；
◣三次样条插值的运行时间相对来说最长，内存消耗比三次插值略少。
它生成的结果平滑性最好。
但是，如果输入数据不很均匀，可能会得到意想不到结果；
◣三次插值需要更多内存，而且运行时间比最近邻插值和线性插值的长。
但是，使用此法时，插值数据及其导数都是连续的。
例7－4程序代码：
x=1:
10;
y=[1800778518500098058827995286700598];
x1=1:
0.1:
y1=interp1（x,y,x1,'
）;
%一维最近邻插值
plot（x1,y1）
y2=interp1（x,y,x1,'
%一维线性插值
plot（x1,y2,'
color'
'
r'
）
运行结果见图7－2。
图7－2 一维最近邻插值和线性插值
例7－5程序代码：
plot（x,y）;
%三次样条插值
plot（x1,y1,'
运行结果见图7－3。
图7－3 三次样条插值
⒉基于FFT的插值
函数interpft用基于FFT的方法进行一维插值。
本方法计算包含周期函数值的向量的傅里叶变换。
然后，它用更多的点计算逆傅里叶变换。
该函数的调用形式为：
y=interpft（x,n）
其中，x是一个包含周期函数值的向量，这些值在等间隔的点上采集。
n是样本大小。
7.3.2二维插值
二维插值在图像处理和可视化方面有着很重要的应用。
MATLAB用函数interp2进行二维插值。
该函数的一般形式为：
ZI=interp2（X,Y,Z,XI,YI,method）
其中，Z是一个矩阵数组，包含二维函数的值，X和Y为大小相同的数组，包含相对于Z的给定值。
XI和YI为包含插值点数据的矩阵，method表示插值方法，为可选参数。
MATLAB提供了三种不同的插值方法进行二维插值：
该方法用分区域常数曲面拟合数据，插值点的值是最近点的值；
●双线性插值（method='
该方法用双线性曲面拟合数据点，插值点的值是四个最近的值的组合。
本方法是分区域双线性的，比双三次插值法快，并且内存消耗更少；
●双三次插值（method='
该方法用双三次曲面拟合数据点，插值点的值是16个最近点的值的组合。
本方法是分区域三次的，结果的平滑性比前面两种的都好。
所有这些方法都要求X和Y数据是单调的，即从点到点，要么总是递增的，要么总是递减的。
应该用meshgrid函数准备这些矩阵。
例7－6程序代码：
%低分辨率的peaks函数图形
[x,y]=meshgrid（-4:
1:
4）;
z=peaks（x,y）;
surf（x,y,z）
运行结果见图7－4。
图7－4 低分辨率的peaks函数图形
例7－7程序代码：
[xI,yI]=meshgrid（-4:
0.25:
zI=interp2（x,y,z,xI,yI,'
%二维最近邻插值
surf（xI,yI,zI）
运行结果见图7－5。
图7－5 二维最近邻插值
例7－8程序代码：
%%二维双线性插值,'
可省略.
运行结果见图7－6。
图7－6 二维双线性插值
例7－9程序代码：
%二维双三次插值
运行结果见图7－7。
图7－7 二维双三次插值
7.3.3多维插值
MATLAB提供了几种多维数据的插值函数，如表7－1中所示。
表7－1 多维数据的插值函数
函 数
描 述
interp3
三维数据插值
interpn
多维数据插值
ndgrid
多维数据网格化
函数interp3进行三维插值。
计算三维样本集V中数据点之间的值。
VI=interp3（X,Y,Z,V,XI,YI,ZI,method）
其中，X，Y和Z指定数据点；
V包含与X，Y和Z对应的值；
XI，YI和ZI为interp3函数对V中数据进行插值的点。
对于超出范围的值，interp3函数返回NaN。
method表示插值方法，为可选参数。
对于三维数据，有三种不同的插值方法。
该方法选择最近点的值；
该方法基于最近的8个点进行分区域线性插值；
该方法基于最近的64个点进行分区域三次插值。
用interpn函数进行更高维数据的插值，该函数的常用形式为：
VI=interpn（X1,X2,X3,…,V,Y1,Y2,Y3,…,method）
其中，Y1,Y2,Y3,…为interpn函数对V中数据进行插值的点。
对于超出范围的值，interpn函数返回NaN。
高维数据插值，同样有最近邻插值、线性插值和三次插值三种方法。
用ndgrid函数为高维函数评价和插值生成数据数组。
该函数将一系列输入向量指定的图域转换为一系列输出数组。
第i维是输入向量xi的拷贝。
ndgrid函数的语法格式为：
[X1,X2,X3,…]=ndgrid（x1,x2,x3,…）
7.4曲线拟合
所谓曲线拟合，就是利用两个或多个变量的离散点，用平滑曲线来拟合它们之间的关系。
根据拟合方法的不同，有参数拟合和非参数拟合之分。
参数拟合，曲线不要求通过所有点，采用最小二乘法；
非参数拟合，要求曲线通过所有点，采用插值法。
由于曲线拟合是数据分析最常见的任务之一，MATLAB提供了多种函数和工具来进行曲线拟合，另外还有曲线拟合工具箱。
7.4.1最小二乘法
最小二乘法通过最小化残差的平方和来获得待定系数的估计。
第
个数据点的残差定义为测量响应值
和拟合响应值
之间的差值，即
，残差的平方和
常见的最小二乘法包括线性最小二乘、加权线性最小二乘、稳健最小二乘和非线性最小二乘等。
求解非线性最小二乘问题的Gauss-Newton法和Levenberg-Marquart法是老牌算法。
7.4.2多项式曲线拟合
用polyfit函数计算拟合数据集的多项式在最小二乘意义上的系数，调用形式为：
p=polyfit（x,y,n）
x和y是包含要拟合的
和
数据的向量，n是多项式的阶次。
例7－10程序代码：
y=[13311333818711723308151138011];
P=polyfit（x,y,4）%多项式曲线拟合，返回多项式的系数。
P=
1.0000-2.0000-0.00001.00001.0000
7.4.3相关工具
MATLAB支持用基本拟合界面进行拟合。
该拟合界面具有拟合快速，操作简便的优势。
它具有如下功能：
●使用样条插值、分段三次艾尔米特插值（PCHIP）或者是1到10阶的多项式插值进行数据拟合；
●利用一组数据可以同时作多条拟合曲线；
●可以绘制残差图；
●可以检查拟合结果；
●可以对拟合进行内插或外推；
●用拟合结果和标准残差在图中进行注释；
●可以将拟合和计算结果保存到MATLAB工作空间。
下面结合一个具体的例子加以说明。
例7－11 某商店某种产品的销售量如表7－2所示。
表7－2 某产品销售量资料
年份
2001
2002
2003
2004
2005
2006
2007
2008
2009
销售量（万件）
10.0
18.0
25.0
30.5
35.0
38.0
40.0
39.5
请用多项式曲线拟合上述数据。
按照下面的步骤进行操作。
⑴用上述数据绘散点图和线形图的组合图（见图7－8）；
%说明多项式曲线拟合界面的使用方法
x=2001:
2009;
y=[10182530.535384039.538];
scatter（x,y,'
filled'
）%绘散点图
plot（x,y）
图7－8 散点图和线形图的组合图
⑵在图形窗口的“Tools”菜单中单击“BasicFitting”菜单选项；
⑶两次单击“→”按钮。
打开的曲线拟合界面“BasicFitting”对话框如图7－9所示。
基本拟合界面中各选项的功能包括：
●Selectdata下拉式列表框：
该列表框中显示了图形窗口中图形用到的所有数据集的名称，在其中选择要拟合的数据。
一次只能选择一组数据，但在一组数据里可以同时拟合多条曲线。
●Centerandscalexdata单选钮：
选择此项以后，数据中心化为具有0均值，比例化为具有单位标准差。
对数据进行中心化和比例化，可以提高后面数值计算的精度。
●Plotfits：
使用本面板，可以用图形查看当前数据集的一种或多种拟合结果。
●Checktodisplayfitsonfigure：
选择当前数据集的拟合类型。
有两种拟合类型可供选择，即插值和多项式。
进行三次样条插值使用Spline函数，保形（shape-preserving）插值使用pchip函数（三次插值）。
多项式拟合使用polyfit函数。
可以选择多种拟合类型。
●Showequations单选钮：
选此项，在拟合图形上显示方程。
●Plotresiduals单选钮：
选此项，显示拟合曲线的残差，可用条形图、散点图或线形图显示。
●Shownormofresiduals单选钮：
选此项，显示残差的范数。
残差的范数是表示拟合优度的一个统计量，值越小，表示拟合程度越高。
用norm函数进行计算，即norm（V,2），其中V为残差向量。
图7－9 基本拟合界面
●Numericalresults方框：
使用该面板，可以在不绘拟合图的情况下探察对当前数据集进行单次拟合的数值结果。
●Fit：
选择拟合当前数据集的方程。
拟合结果显示在菜单下面的列表框中。
注意，在该菜单中选择方程并不影响“Plotfits”面板的状态。
所以。
如果试图在数据散点图中显示拟合曲线，可能需要在“Plotfits”面板中选择有关的核选框。
●Coefficientsandnormofresiduals：
显示“Fit”中选择的方程的计算结果。
●Savetoworkspace…：
打开一个对话框，使用它将拟合结果保存到工作空间变量中。
●FindY=f（X）：
对当前的拟合进行内插或外推。
●Entervalue（s）：
输入一个MATLAB表达式，进行拟合计算。
单击“Evaluate”按钮以后计算表达式，结果显示在有关的表格中。
当前拟合显示在“Fit”菜单中。
●Plotevaluatedresults：
选择此项，计算结果显示到图中的数据点上。
⑷在基本拟合界面中作以下设置（见图7－10）：
●用二次多项式拟合数据；
●在拟合图上显示方程；
●将拟合残差作为条形图显示，并将残差的图形作为子图显示；
●显示残差的范数。
当前数据集
用二次多项
式拟合数据
显示方程
绘残差图
显示残差
的范数
图7－10 进行选项设置
⑸利用“Plotfits”面板可以可视地探察当前数据集的多个拟合图形。
为了进行比较，通过选择合适的核选框来拟合数据的其他方程。
如果某个方程生成的结果在数值上不精确，MATLAB会显示一个警告信息。
此时，应该选择“Centerandscalexdata”核选框来改进数值精度。
⑹最后将拟合结果和残差显示在图7－11中。
图例显示了数据集和方程的名称。
如果图例覆盖了图形的一部分，可以通过单击和拖拉操作将它移到其他地方。
图7－11 拟合结果和残差图
⑺可以指定一个包含x值的向量，计算这些位置上的拟合值。
将该向量输入到“Evaluate”按钮左边的文本框中，然后单击“Evaluate”按钮。
例如，输入向量2010:
2015后，再单击“Evaluate”按钮，即可得到2010年到2015年的销售量的估计。
如图7－12所示。
⑻选择“Plotevaluatedresults”核选框，显示当前图形中对应当前数据的计算值，如图7－13所示
图7－12 显示x值和对应的拟合值
图7－13 拟合结果和对应的残差图
7.5数值微分
利用MATLAB函数，可以实现数值微分、梯度运算。
7.5.1数值微分运算
5.2一节中用diff函数实现了符号微分运算。
实际上，diff函数还可以实现数值微分运算，即求相邻元素之间的差值。
diff函数的调用格式如下：
Y=diff（X,n,dim）
n和dim均为大于或等于1的整数，n和dim的默认值均为1。
当dim大于X的维数时，将返回空数组；
当n大于或等于X的第dim维的长度时，也返回空数组。
●若X是向量，则可省略dim，diff（X）返回X的一阶差分；
diff（X,n）（n大于1且小于X的长度）返回X的n阶差分（等价于递归n次一阶差分）。
●若X是p行q列矩阵，则diff（X）返回X的每列元素的一阶差分；
diff（X,n）（n大于1且小于p）返回X的每列元素的n阶差分；
diff（X,n,2）（n大于1且小于q）返回X的每行元素的n阶差分。
●若X是
数组，则diff（X,n,dim）（
）返回X的第dim维元素的n阶差分。
例7－12 程序代码
x（:
:
1）=[123;
456;
789];
2）=[0.10.20.3;
0.40.50.6;
0.70.80.9]
D11=diff（x）
D12=diff（x,1,2）
D13=diff（x,1,3）
D21=diff（x,2）
1）=
123
456
789
2）=
0.10000.20000.3000
0.40000.50000.6000
0.70000.80000.9000
D11（:
333
0.30000.30000.3000
D12（:
11
0.10000.1000
D13=
-0.9000-1.8000-2.7000
-3.6000-4.5000-5.4000
-6.3000-7.2000-8.1000
D21（:
000
1.0e-015*
-0.11100.05550.0555
7.5.2数值梯度运算
用gradient函数进行数值梯度运算。
gradient函数的调用格式如下：
Fx=gradient（F）
[Fx,Fy]=gradient（F）
[Fx,Fy,Fz,...]=gradient（F）
[...]=gradient（F,h）
[...
升级会员 