内容发布更新时间 : 2024/8/27 21:28:53星期一 下面是文章的全部内容请认真阅读。
中国民航大学《数值分析》课程报告
基于对插值法与最小二乘法以及数值积分的MATLAB实训
刘超佳董巧丽
（中国民航大学机场学院交通运输工程（道铁方向）-2018082040）
摘要：插值是指已知某函数的在若干离散点上的函数值或者导数信息，通过求解该函数中待定形式的插值函数以及待定系数，使得该函数在给定离散点上满足约束。插值函数又叫作基函数，如果该基函数定义在整个定义域上，叫作全域基，否则叫作分域基。如果约束条件中只有函数值的约束，叫作Lagrange插值，否则叫作Hermite插值。从几何意义上将，拟合是给定了空间中的一些点，找到一个已知形式未知参数的连续曲面来最大限度地逼近这些点而插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。本文着重对复合梯形求积公式、复合Simpson公式、复合Cotes公式及最小二乘法数据拟合结合MATLAB进行了练习，以巩固课堂所学，为实验应用做好准备。 关键词：数值分析插值MATLAB函数SimpsonCotes最小二乘法
MATLAB training based on interpolation, least square method and
numerical integration Liu Chao-jiaDong Qiao-li
(School of airport engineering,2018082040,Civil Aviation University of China)
Abstract:Interpolation refers to the function value or derivative information of a given function at a number of discrete points. By solving the undetermined interpolation function and undetermined coefficient of the function, the function satisfies the constraint at a given discrete point.An interpolation function is also called a basis function, and if the basis function is defined on the entire domain, it's called a global basis, otherwise it's called a subdomain basis.If the constraint condition only contains the constraint of function value, it is called Lagrange interpolation; otherwise, it is called Hermite interpolation.In the geometric sense, fitting is to find a continuous surface of known form and unknown parameters to approximate the points to the greatest extent given some points in space, while interpolation is to find a continuous surface (or several pieces of smooth) to pass through the points.This paper focuses on the practice of the compound trapezoidal quadrulation formula, the compound Simpson formula, the compound Cotes formula and the least squares data fitting combined with MATLAB, so as to consolidate what we have learned in class and prepare for the experimental application.
Key words: numerical analysis; interpolation; MATLAB function; Simpson Cotes; least square method
Author resume: LIU Chao-jia(1996-),Master,Research direction: airport pavement structure
亚齐发现了第一颗小行星谷神星。经过40
1.引言
天的跟踪观测后，由于谷神星运行至太阳
1801年，意大利天文学家朱赛普·皮背后，使得皮亚齐失去了谷神星的位置。
中国民航大学《数值分析》课程报告
随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星，但是根据大多数人计算的结果来寻找谷神星都没有结果。时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥伯斯根据高斯计算出来的轨道重新发现了谷神星。
高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中，而法国科学家勒让德于1806年独立发现“最小二乘法”，但因不为世人所知而默默无闻。两人曾为谁最早创立最小二乘法原理发生争执。1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，见高斯-马尔可夫定理。
2.结合MATLAB应用几种插值法
2.1对于函数f(x) = sin(x) ，用复化梯形公式、复化Simpson公式及复化Cotes公式计算积分，并比较其误差。
首先，在MATLAB输入以下程序定义函数：
function y=f(x)
y=sin(x); 2.1.1复合梯形求积公式的代码 2.1.2复合Cotes公式的代码 function Cn = Cn(a,b,n) format long h = (b-a)/n; sum1 = 0; sum2 = 0; fori=0:n-1
sum1=sum1+32*f(a+(i+1/4).*h)+12* f(a+(i+1/2).*h)+32*f(a+(i+3/4).*h); end
for j = 1:n-1
sum2 = sum2 + 14*f(a+j.*h);
endCn=h/90*(7*f(a)+sum1+sum2+7*f(b));
2.1.3不同区间和分割次数下的精确值和近似值的比较代码
y=[0.150920964949986,0.106044306328975,1.260000791649389,1.122733727
function Tn=Tn(a,b,n) format long h=(b-a)/n; sum=0; for k=1:n-1
sum=sum+f(a+k.*h);
end
Tn=(f(a)+2*sum+f(b))*h/2; end
2.2 复合Simpson公式的代码 function Sn = Sn(a,b,n) format long h = (b-a)/n; sum1 = 0; sum2 = 0;
for i = 0:n-1
sum1=sum1+f(a+(i+1/2).*h); end
for j = 1:n-1
sum2 = sum2 + f(a+j.*h); end
Sn=h/6*(f(a)+4*sum1+2*sum2+f(b)); 670800,0.185427920752058,0.955618503172607];
Tn=[0.144707438957434,0.103357516942886,0.809037206720437,1.074587530720150,0.176653381231435,0.874739367319743];
Sn=[0.150933735388946,0.106047691917755,1.267931587831384,1.122836754412495,0.185448626084829 ,0.956086587141625];
Cn=[0.150920958200157,0.106044290463326,1.259953005777527,1.122733636186552,0.185427890841067,0.955744078179222];
plot(y,'-k*'); hold on; plot(Tn,'-g.');
中国民航大学《数值分析》课程报告
hold on;
plot(Sn,'-bo'); hold on;
plot(Cn,'-r+');
legend('精确值','Tn','Sn','Cn'); 2.2.1图像
2.2.2结果分析及误差分析
在上面的数值计算中，复合Cotes公式计算所得的误差最小，其次是复合Simpson公式，误差最大的是复合梯形公式。
2.2.3精确值比较
它们的精度有很大的差别：与精确值比较，复合梯形公式的结果只有一位有效数字，复合Simpson公式的结果有四位有效数字，而复合Cotes公式的结果却有七位有效数字。
由此可知，复合Cotes公式的代数精度最高，复合Simpson其次，复合梯形公式的代数精度最低。 2.2.4收敛性分析
2.2.5稳定性分析
可知三个求积公式的求积系数都为正，所以复合梯形公式、复合Simpson公式和复合Cotes公式均是稳定的。 2.3用辛普森公式求积分 ??
???????????并估计误差。
function [y e]= Simpson(f,a,b,M) % f被积函数；a积分下限；b积分上限；M子区间个数（将x分为多少个区间）
h=(b-a)/(2*M); s1=0; s2=0; for i=1:M x=a+(2*i-1*h; s1=s1+feval(f,x); end
for j=1:(M-1) x=a+2*j*h;
s2=s2+feval(f,x); end
y=h/3*(feval(f,a)+2*s2+4*s1+feval(f,b));
e=quad(f,0,1)-y;%误差（运行后不显
中国民航大学《数值分析》课程报告
示，把这行命令的分号去掉就运行可以显示误差）
我的结果是
>> Simpson(f,a,b,M) e =-4.4409e-016 ans =2.2183
2.3.1复合simpson求积算例 clc
format long clear
a=0;b=1;p=0; n=4;m=0; h=(b-a)/n; for i=1:n-1 x=0; x=a+h*i;
m=m+sin(x)/x; end
h=(b-a)/(2*n); for k=1:2:2*n
p=p+sin(a+k*h)/(a+k*h); end
Sn=h/3*[1+4*p+2*m+sin(b)/b] %辛普森公式结果
result=sinint(1) %计算精确值 format short
四．复核梯形公式求积算例
x=0:0.00001:1;%x用来储存积分点 y=(x+1).*sin(x);%y用来求解积分点x处的函数值 I=trapz(x,y)
I =0.7608663730793 syms x
y=(x+1)*sin(x);%被积函数表达式 II=int(y,0,1)
II =sin(1) - 2*cos(1) + 1 %II 即为该被积函数的解析解 II E=eval(II)
II_E =0.760866373071617 %II的数值解
%可以看出梯形求积公式在步长等于0.00001的情况下，数值积分的解与解析解的数值能达到小数点后11位保持一致。
2.4 Romberg算法算例 clear;clc;close all; format long
% %被积函数为f(x)=4/(1+x^2);积分区间为[0,1]
% b=1;a=0;h=b-a;eps=10^(-5);
%误差界eps%被积函数为f(x)=(x^3+sin(x))/x;积分区间为[0.3,0.8] %误差界eps=10^(-5)
b=0.8;a=0.3;h=b-a;eps=10^(-5);%误差界eps
kmax=10;%最大递推次数
T1=h/2*((a^3+sin(a))/a+(b^3+sin(b))/b)
S1=0;C1=0;C2=0;R1=0;R2=0; for k=1:kmax h=(b-a)/2^k; i=1:2^(k-1); x=a+(2*i-1)*h;
fx=sum((x.^3+sin(x))./x); T2=T1/2+fx*h S2=T2+(T2-T1)/3 if(k<3)
if k==2
C2=S2+(S2-S1)/15 end else
C2=S2+(S2-S1)/15 R2=C2+(C2-C1)/63 if abs(R2-R1) T1=T2;S1=S2;C1=C2; end fprintf('所求积分I=%9.8f\\n',R2); 2.5 牛顿（Newton)插值算法 拉格朗日插值的优点是格式整齐和 中国民航大学《数值分析》课程报告 规范，有误差估计方式，它的缺点是没有承袭性，当需要增加节点时，必须重新计算插值的基函数。 2.5.1牛顿插值公式 2.5.2牛顿插值法的MATLAB的综合程序 编写M文件 %newton.m %求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序 %输入的量:X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量，Y是纵坐标向量， %x是以向量形式输入的m个插值点，M在[a,b]上满足｜f~(n+1)(x)｜≤M %注：f~(n+1)(x)表示f(x)的n+1阶导数 %输出的量：向量y是向量x处的插值，误差限R，n次牛顿插值多项式L及其系数向量C， %差商的矩阵A function[y,R,A,C,L] = newton(X,Y,x,M) n = length(X); m = length(x); for t = 1 : m z = x(t); A = zeros(n,n); A(:,1) = Y'; s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0; for j = 2 : n for i = j : n A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1)); end q1 = abs(q1*(z-X(j-1))); c1 = c1 * j; end C = A(n, n); q1 = abs(q1*(z-X(n))); for k = (n-1):-1:1 C = conv(C, poly(X(k))); d = length(C); C(d) = C(d) + A(k,k);%在最后一维，也就是常数项加上新的差商 end y(t) = polyval(C,z); R(t) = M * q1 / c1; end L = poly2sym(C); 例：用Lagrange插值来求sinx在某点的值，并估计其误差，已知sin0° = 0, sin30° = 0.5, sin45° = 0.7071, sin60° = 0.8660, sin90° = 1. X Y 0 pi/6 0 0.5 pi/4 0.7071 pi/3 0.8660 pi/2 1 解： >> X = [0 pi/6 pi/4 pi/3 pi/2]; >> Y = [0 0.5 0.7071 0.8660 1]; >> x = linspace(0,pi,50); >> M = 1; >> [y,R,A,C,L] = newton(X, Y, x, M); >> y1 = sin(x); >> errorbar(x,y,R,'.g') >> hold on