《数值分析实验报告--实验3--函数逼近与曲线拟合.pdf》由会员分享，可在线阅读，更多相关《数值分析实验报告--实验3--函数逼近与曲线拟合.pdf（13页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、数值分析实验三三：函数逼近与曲线拟合 1 / 13 数值分析实验数值分析实验三三：函数逼近与曲线拟合：函数逼近与曲线拟合 1 曲线逼近方法的比较 1.1 问题描述 曲线的拟合和插值， 是逼近函数的基本方法， 每种方法具有各自的特点和特定的适用范围，实际工作中合理选择方法是重要的。 考虑实验 2.1 中的著名问题。下面的 MATLAB 程序给出了该函数的二次和三次拟合多项式。 x=-1:0.2:1; y=1./(1+25*x.*x); xx=-1:0.02:1; p2=polyfit(x,y,2); yy=polyval(p2,xx); plot(x,y,o,xx,yy); xlabel(x);
2、 ylabel(y); hold on; p3=polyfit(x,y,3); yy=polyval(p3,xx); plot(x,y,o,xx,yy); hold off; 实验要求： (1) 将拟合的结果与拉格朗日插值及样条插值的结果比较。 (2) 归纳总结数值实验结果，试定性地说明函数逼近各种方法的适用范围，及实际应用中选择方法应注意的问题。 1.2 算法设计 对于曲线拟合，这里主要使用了多项式拟合，使用 Matlab 的 polyfit 函数，可以根据需要选用不同的拟合次数。 然后将拟合的结果和插值法进行比较即可。 本实验的算法比较简单，此处不再详述，可以参见给出的 Matlab 脚本
3、文件。 1.3 实验结果 1.3.1 多项式拟合 1.3.1.1 多项式拟合函数 polyfit 和拟合次数 N 的关系 数值分析实验三三：函数逼近与曲线拟合 2 / 13 首先使用 polyfit 函数对 f(x)进行拟合。为了便于和实验 2.1 相比较，这里采取相同的参数，即将拟合区间-1,1等分为 10 段，使用每一段区间端点作为拟合的数据点。分别画出拟合多项式的次数 N=2、3、4、6、8、10 时，f(x)和多项式函数的图像，如图 1 所示。Matlab脚本文件为 Experiment3_1_1.m。 Figure 1 多项式拟合与拟合次数 N 的关系 可以看出，拟合次数 N=2 和
4、 3 时，拟合效果很差。增大拟合次数，N=4、6、8 时，拟合效果有明显提高，但是 N 太大时，在区间两端附近会出现和高次拉格朗日插值函数类似的龙格现象。所以，对于多项式插值来说，适当地选择插值次数，对于插值结果的准确度至关重要。通常，在使用多项式拟合处理科学实验得到的数据时，我们需要通过理论分析和推导，事先确定拟合的次数，这样会大大提高拟合的准确度和效率。 1.3.1.2 多项式拟合与插值法比较 取多项式拟合的次数 N=10，朗格朗日插值和样条函数插值区间个数也都为 10。画出对于函数 f(x)，三者的结果如图 2 所示。Matlab 脚本文件为 Experiment3_1_2.m。 可以看
5、出拉格朗日插值和多项式拟合近乎完全重合， 二者在区间端点附近都存在龙格现象。而样条函数插值逼近原函数的效果较好。 1.3.2 函数逼近方法总结 函数逼近的数值方法包括插值和拟合。 插值法经过数据点， 而拟合方法通常不会经过所有的数据点。在处理实验得到的数据时，如果我们预先知道数据服从怎样的分布，如线性和指数函数，可以选用多项式拟合和指数拟合。当然我们应该首先观察数据的分布，然后选择使用什么样的你和方法。如本实验中利用 f(x)得到的数据点，明显不符合多项式函数，所以选择多项式拟合的结果很不理想，应该考虑选用高斯曲线模型进行拟合。 使用插值法时，要注意高次拉格朗日插值时，会出现龙格现象，我们应该
6、尽量增大插值数值分析实验三三：函数逼近与曲线拟合 3 / 13 的区间，使得我们关注的部分尽可能落在区间中部，这样可以尽量减小龙格现象的影响。三次样条插值函数理论上是收敛的， 因此可以获得较好的插值效果， 但是由于它在每一个小区间上的表达式都不同，有时会显得繁琐。 Figure 2 多项式拟合和插值法的比较 1.4 分析总结 本次实验比较了多项式拟合与几种插值方法的结果。 不同的拟合方法和插值方法各有优点和缺点，实际应用时应该根据需要，合理地进行选择。 2 最小二乘拟合的经验公式和模型 2.1 问题描述 2.1.1 已知经验公式 某类疾病发病率为y和年龄段x(每五年为一段，例如 05 岁为第一
7、段，610 岁为第二段)之间有形如bxyae的经验关系，观测得到的数据表如下： x 1 2 3 4 5 6 7 8 9 y 0.898 2.38 3.07 1.84 2.02 1.94 2.22 2.77 4.02 x 10 11 12 13 14 15 16 17 18 19 y 4.76 5.46 6.53 10.9 16.5 22.5 35.7 50.6 61.6 81.8 实验要求： （1）用最小二乘法确定模型bxyae中的参数a和b(提示函数：lsqcurvefit，lsqnonlin)。 数值分析实验三三：函数逼近与曲线拟合 4 / 13 （2）利用 MATLAB 画出离散数据及拟
8、合函数bxyae图形。 （3）利用 MATLAB 画出离散点处的误差图，并计算相应的均方误差。 2.1.2 最小二乘拟合模型未知 某年美国轿车价格的调查资料如表，其中ix表示轿车的使用年数，iy表示相应的平均价格。 实验要求：试分析用什么形式的曲线来拟合表中的数据，并预测使用 4.5 年后轿车的平均价格大致为多少？ ix 1 2 3 4 5 6 7 8 9 10 iy 2615 1943 1494 1087 765 538 484 290 226 204 2.2 算法设计 主要使用 Matlab 的 lsqcurvefit 函数进行最小二乘法拟合，算法设计比较简单，可以参考实验结果部分给出的
9、Matlab 脚本，此处不再详述。 2.3 实验结果 2.3.1 已知经验公式 使用 Matlab 的 lsqcurvefit 函数进行拟合， 由于已知经验公式，因此拟合函数选用bxyae形式的指数函数即可。Matlab 脚本文件为 Experiment3_2_1.m。 Figure 3 疾病发病率拟合结果 数值分析实验三三：函数逼近与曲线拟合 5 / 13 Figure 4 拟合结果误差图 拟合得到的表达式为 0.30900.2369xye 画出离散数据点和拟合函数的图形，如图 3 所示。误差图如图 4 所示。均方误差为 1.9857 2.3.2 最小二乘拟合模型未知 使用最小二乘法进行拟合
10、的时候， 如果数据分布的模型未知， 首先应该画出数据散点图，通过经验判断出数据的大概分布，然后选择经验公式进行拟合。本实验的 Matlab 脚本文件为 Experiment3_2_2.m。 将汽车的平均价格和使用年数画出散点图，如图 5 所示。 通过观察，可以看出该数据大致上服从指数分布，我们不妨设iy和ix服从以下函数关系： bxyaec 然后使用 Matlab 的 lsqcurvefit 函数进行拟合。从散点图可以看出，当 x=0 时，y 大概在3000 附近，所以我们设 a，b，c 的初始值分别为 3000,0,0。最后拟合的结果如图 6 所示。 其表达式为 0.295373544.21
11、1613.0295xye x=4.5，即使用年数为 4.5 年时，汽车的价格为 925.11 美元。 数值分析实验三三：函数逼近与曲线拟合 6 / 13 Figure 5 汽车平均价格和使用年数的散点图 Figure 6 汽车平均价格和使用年数拟合结果 2.4 分析总结 本实验主要对已知模型和未知模型良两种情况， 使用最小二乘法进行了拟合。 当数据模型未知时，需要根据数据的散点图和经验判断出数据的分布模型，并不断尝试，最终找到最优的拟合方法。通过实验掌握了 Matlab 的 lsqcurvefit 函数的使用，深刻理解了最小二乘法数值分析实验三三：函数逼近与曲线拟合 7 / 13 拟合的思想。
12、 3 研究最佳平方逼近多项式的收敛性质 3.1 问题描述 取函数( )xf xe，在-1，1上以勒让德多项式为基函数，对于0,1,10n构造最佳平方逼近多项式( )np x， 令( )( )( )nnxf xp x， 将( ) nxx的曲线画在一个图上。 令11max( )( )nnxEf xp x ，画出nEn的曲线。做出nEn之间的最小二乘曲线，能否提出关于收敛性的猜测。 3.2 算法设计 课本 73 页的（2.9）式给出了勒让德多项式的递推公式，又 01P x ， 1P xx，根据递推公式可以依次求出 P2，P3，P10，在 Matlab 中通过循环语句实现。 课本 87 页的（4.9）
13、式给出了基于勒让德多项式展开的最佳平方逼近函数的表达式。 本实验主要使用了 Matlab 的符号变量、符号表达式和定积分 int 函数进行运算，使用循环语句遍历 n=0 到 n=10 的不同情况，然后使用绘图函数得到各种情况下的图形。详细算法可以 Matlab 脚本文件 Experiment3_3.m。 3.3 实验结果 按照题目中的要求，对于 n=0,1，10 时，分别构造出 Pn(x)，并计算误差函数( )( )( )nnxf xp x，画出( ) nxx曲线，如图 7 所示，由于随着 n 的增加，( ) nxx在急剧降低，因此为了方便观察，这里 y 轴使用对数坐标刻度。nEn曲线如图 8
14、 所示。 从图 7 中可以看出，误差函数的分布随着 n 的增加，急剧降低。从图 8 中可以看出，误差限随着 n 增加，以指数速度衰减。所以，我们可以判断，对函数 f(x)，以勒让德多项式作为基函数进行展开，得到的最佳平方逼近函数，是收敛的。图 8 中误差限与 n 的最小二乘拟合的关系式为 -1.294( )5.6589 e-0.0055293nE n 数值分析实验三三：函数逼近与曲线拟合 8 / 13 Figure 7 n 变化时误差函数曲线 Figure 8 误差限与 n 最小二乘拟合曲线 3.4 分析总结 本次实验针对函数 f(x)=exp(x)，在-1,1以勒让德多项式作为基函数展开，计
15、算得到最佳平方逼近函数，并针对不同的 n，计算了误差函数。通过画出观察和计算误差限与 n 的关系， 得出了误差函数会随着 n 的增加不断减小的结果， 进而验证了最佳平方逼近多项式的收敛性。 数值分析实验三三：函数逼近与曲线拟合 9 / 13 附件：源码 Experiment3_1_1.m % Experiment3_1_1.m % Programmer: Zhang Weizhen % Date: 2016-10-27 clear; format long; x=-1:0.2:1; y=1./(1+25*x.*x); subplot(2,3,6); plot(x,y,o); hold on;
16、grid on; xlabel(x); ylabel(y); xx=-1:0.01:1; yy=1./(1+25*xx.*xx); plot(xx,yy,k); N=10; x_polyfit=-1:0.02:1; pf=polyfit(x,y,N); y_polyfit=polyval(pf,x_polyfit); plot(x_polyfit,y_polyfit,-r); title(N=,num2str(N); legend(fitting points,f(x),polyfit ,num2str(N); hold off; Experiment3_1_2.m % Experiment3
17、_1_2.m % Programmer: Zhang Weizhen % Date: 2016-10-28 clear; format long; n = 10; % number of points 数值分析实验三三：函数逼近与曲线拟合 10 / 13 xi = ones(1,n+1); for i=0:1:n xi(i+1) = -1+2*i/n; % 节点 x 值 end yi = 1./(1+25*xi.*xi); % 节点 y 值 x_point = -1:0.01:1; % 画图点 x 值 f = 1./(1+25*x_point.*x_point); % 画图点原函数 f 值 y
18、_lagrange = f_lagrange(xi,yi,x_point); y_spline = spline(xi,yi,x_point); plot(x_point, f,k, x_point, y_lagrange, -m,x_point, y_spline,-r); hold on; grid on; xlabel(X); ylabel(Y); % 多项式拟合 x=-1:0.2:1; y=1./(1+25*x.*x); plot(x,y,o); x_polyfit=-1:0.02:1; pf=polyfit(x,y,10); y_polyfit=polyval(pf,x_polyfi
19、t); plot(x_polyfit,y_polyfit,-b); legend(f(x),Lagrange 10,Spline 10,fitting points,f(x),polyfit 10); hold off; Experiment3_2_1.m % Experiment3_2_1.m % Programmer: Zhang Weizhen % Date: 2016-10-28 clear; format long; xdata=1:1:19; ydata=0.898 2.38 3.07 1.84 2.02 1.94 2.22 2.77 4.02 4.76 5.46 6.53 10.
20、9. 16.5 22.5 35.7 50.6 61.6 81.8; x0=1 1; x = lsqcurvefit(x,xdata) x(1)*exp(x(2)*xdata), x0, xdata,ydata); xx=1:0.01:19; yy=x(1).*exp(x(2).*xx); 数值分析实验三三：函数逼近与曲线拟合 11 / 13 figure(1); plot(xdata,ydata,sb,MarkerSize,6,MarkerFaceColor,b); hold on; plot(xx,yy,-r); grid on; xlabel(X); ylabel(Y); legend(f
21、itting points,y=,num2str(x(1),*exp(,num2str(x(2),*x) ); hold off; fprintf(a=,num2str(x(1),n); fprintf(b=,num2str(x(2),n); figure(2); y_error = x(1).*exp(x(2).*xdata) - ydata; plot(xdata, y_error,sb,MarkerSize,6,MarkerFaceColor,b); grid on; xlabel(X); ylabel(Y); legend(Y-error); sigma = sqrt(sum(y_er
22、ror.2)/length(xdata); fprintf(sigma=,num2str(sigma),n); Experiment3_2_2.m % Experiment3_2_2.m % Programmer: Zhang Weizhen % Date: 2016-10-28 clear; format long; % ? xdata=1:1:10; figure(1); ydata=2615 1943 1494 1087 765 538 484 290 226 204; plot(xdata,ydata,sb,MarkerSize,6,MarkerFaceColor,b); grid o
23、n; xlabel(X); ylabel(Y); legend(Data points) % ? x0=3000 0 0; x= lsqcurvefit(x,xdata) x(1)*exp(-x(2)*xdata)+x(3), x0, xdata, ydata); 数值分析实验三三：函数逼近与曲线拟合 12 / 13 xx=1:0.01:10; yy=x(1).*exp(-x(2).*xx)+x(3); figure(2); plot(xdata,ydata,sb,MarkerSize,6,MarkerFaceColor,b); hold on; plot(xx,yy,-r); grid on
24、; xlabel(X); ylabel(Y); legend(fitting points,y=,num2str(x(1),*exp(,num2str(-x(2),. *x)+(,num2str(x(3),); hold off; fprintf(a=,num2str(x(1),n); fprintf(b=,num2str(x(2),n); fprintf(c=,num2str(x(3),n); Experiment3_3.m % Experiment3_3.m % Programmer: Zhang Weizhen % Date: 2016-10-29 clear; format long;
25、 % epsilonx 曲线 figure(1); hold on; NN=10; M=cell(1,NN); E=zeros(1,NN+1); for i=0:1:NN N=i; P=zeros(1,N+1,sym); syms x; P(1)=1; P(2)=x; for n=2:1:N+1 P(n+1)=1/(n)*(2*n-1)*x*P(n)-(n-1)*P(n-1); end S=0; for n=1:1:N+1 S=S+int(exp(x)*P(n),x,-1,1)/int(P(n)*P(n),x,-1,1)*P(n); % S 是符号函数 数值分析实验三三：函数逼近与曲线拟合 1
26、3 / 13 end epsilon=abs(exp(x)-S); x=-1:0.01:1; y=eval(epsilon); E(1,i+1)=max(y); plot(x,y); M1,i+1=strcat(e,num2str(N),(x); end grid on; xlabel(X); ylabel(epsilon(n, x); legend(M); set(gca,yscale,log); title(不同 n 时，得到的 n(x)x 图像); hold off; % Enx 曲线，及最小二乘拟合 figure(2); xdata=1:NN+1; plot(xdata,E,sb,Ma
27、rkerSize,6,MarkerFaceColor,b); % 散点图 hold on; x0=0 0 0; % 拟合 x=lsqcurvefit(x,xdata) x(1)*exp(-x(2)*xdata)+x(3), x0, xdata, E); xx=1:0.01:NN+1; yy=x(1).*exp(-x(2).*xx)+x(3); plot(xx,yy,-r); grid on; xlabel(n); ylabel(En); legend(fitting points,E(n)=,num2str(x(1),*exp(,num2str(-x(2),. *n)+(,num2str(x(3),); title(E(n)n 曲线); fprintf(a=,num2str(x(1),n); fprintf(b=,num2str(x(2),n); fprintf(c=,num2str(x(3),n); hold off;