正交多项式最小二乘法拟合
最小二乘法多项式拟合
最小二乘法多项式拟合对于给定的数据点N i y x i i ≤≤1),,(，可用下面的n 阶多项式进行拟合，即∑==+++=nk k k x a x a x a a x f 02210)(为了使拟合出的近似曲线能尽量反映所给数据的变化趋势，要求在所有数据点上的残差|)(|||i i i y x f -=δ都较小。
为达到上述目标，可以令上述偏差的平方和最小，即min ])([)(2121=-=∑∑==iiNi iN i y x f δ称这种方法为最小二乘原则，利用这一原则确定拟合多项式)(x f 的方法即为最小二乘法多项式拟合。
确定上述多项式的过程也就是确定)(x f 中的系数n k a k ≤≤0,的过程，根据最小二乘原则，则偏差平方和应该是这些系数的函数，即min ])([)(),,,(212110=-==∑∑==i i Ni i N i n y x f a a a S δ为使上式取值最小，则其关于n k a k ≤≤0,的一阶导数应该为零，即有∑∑∑∑=====⇒=-⇒=-=∂∂Ni i N i i i i N i i i N i y x f y x f y x f a S11110)(0])([0])([2 ∑∑∑∑=====⇒=-⇒=-=∂∂N i i i N i i i i i N i i i i N i i y x x f x y x f x y x f x a S11111)(0])([0])([2∑∑∑∑=====⇒=-⇒=-=∂∂N i i k i N i i ki i i N i k i i i N i k i k y x x f x y x f x y x f kx a S 1111)(0])([0])([2∑∑∑∑=====⇒=-⇒=-=∂∂N i i n i N i i ni i i N i n i i i N i n i n y x x f x y x f x y x f nx a S 1111)(0])([0])([2 将上面各等式写成方程组的形式可有∑∑∑∑∑∑=======++++⇒=Ni i N i n in N i iN i i Ni iN i iy x a x a x a N a yx f 1112211011)(∑∑∑∑∑∑==+=====++++⇒=Ni i i Ni n in Ni iNi ii Ni iiNi iiy x xa x a x a x a yx x f x 111132121011)(∑∑∑∑∑∑==+=+=+===++++⇒=Ni i k i Ni k n in Ni k iNi k ik iNi i k i Ni i k iy x xa xa xa x a y x x f x11122111011)(∑∑∑∑∑∑===+=+===++++⇒=Ni i n i Ni n in Ni n iNi n in iNi i n i Ni i n iy x xa xa xa x a y x x f x112122111011)(写成矩阵形式有⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⋅⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛∑∑∑∑∑∑∑∑∑∑∑∑∑∑∑∑∑∑∑======+=+==+==+===+=====N i i n i Ni ik i N i i i N i i n k N i ni Ni k n iNi n iNi ni N i k n i N i ki N i k i N i kiN i n i Ni k iNi i N i i Ni niNi k iNi iy x y x y x y a a a a x xxx x x x x x xxx x x x N 111110121111112111111121111上述方程组可以通过克莱姆法则来计算，从而解出各系数n k a k ≤≤0,得到拟合方程。
数值分析--32正交多项式与最小二乘拟合
j 0 n
{ j(x) = x j } 对应代数多项式 /* algebraic polynomial */
{ j(x) = cos jx }、{ j(x) = sin jx } { j(x), j(x) }对应三 角多项式 /* trigonometric polynomial */ { j(x) = e kj x , ki kj } 对应指数多项式 /* exponential
n
定理 Ba = c 存在唯一解 0(x), 1(x), … , n(x) 线性无关。
证明：若存在一组系数 {i } 使得 0 0 + 1 1 + ... + n n 0 则等式两边分别与0, 1, … , n作内积，得到：
0 ( 0 , 0 ) + 1 ( 1 , 0 ) + ... + n ( n , 0 ) 0 ( , ) + ( , ) + ... + ( , ) 0 1 1 1 n n 1 0 0 1 . . . ( , ) + ( , ) + ... + ( , ) 0 1 1 n n n n 0 0 n
1 2 49 3 y P( x) x + x 2 10 2
cond ( B ) 7623
|| B || 484,
|| B
1
||
63 4
§6 Orthogonal Polynomials & L-S Approximation
j 例：连续型拟合中，取 j ( x) x , ( x) 1, y( x) C[0, 1]
数据拟合的非线性模型多项式拟合的正交化方法超定方程组最小二乘
( x) a00 ( x) a11( x) ann ( x)
超定方程组:
0 ( x1 ) 1( x1 ) n ( x1 ) a0 y1
0
(
x2
)
1( x2 )
n
(
x2
)
a1
y2
数据拟合的非线性模型
观测数据
x
x1 x2 ·········· xm
f
y1 y2 ·········· ym
求拟合函数 f(x, c0, c1, ···, cn )满足 m [ f ( xi , c0 , c1 ,, cn ) yi ]2 min
i 1
例1.已知人口统计数据
年 1991 1992 1993 1994 1995 1996
数量 11.58 11.72 11.85 11.98 12.11 12.24
利用最小二乘法求指数拟合 y = c e a x
6
S(a, c) [c exp( ax j ) y j ]2 min j 1
2/16
指数函数拟合人口统计数据(单位：亿)
t 1991 1992 1993 1994 1995 1996
12/16
矩阵G的QR(正交三角)分解算法
G = ( gij)m×n m > n
将矩阵按列分块，记为
G = [g1, g2 , ······,gn ]
分解算法 ① q1=g1 ② q2 = g2 – r12q1, 其中 r12 = (g2，q1)/(q1，q1)
···········································
最小二乘法多项式拟合原理
最小二乘法多项式拟合原理最小二乘法多项式拟合原理最小二乘法是一种数学方法，用于寻找一个函数，使得该函数与已知数据点的残差平方和最小化。
尤其在数据分析和统计学中广泛应用，其中特别重要的应用是曲线拟合。
本文将介绍最小二乘法在多项式拟合中的原理。
多项式拟合多项式拟合是一种常见的曲线拟合方法，它将数据点逼近为一个固定次数的多项式。
假设有N个数据点(x1,y1),(x2,y2),…,(xN,yN)，希望找到一个关于x的M次多项式函数y=a0+a1x+a2x^2+...+aMx^M，最小化拟合曲线与数据点之间的残差平方和，即S(a0,a1,…,aM)=∑i=1N(yi−P(x))2其中P(x)=a0+a1x+a2x^2+...+aMx^M。
最小二乘法最小二乘法是一种优化方法，通过最小化残差平方和，寻找最优的拟合函数参数。
在多项式拟合中，残差平方和的最小值可以通过相应的求导数为零来计算拟合函数参数。
设残差平方和S的导数为零得到的方程组为∑xi0,…,xiMaM=∑yi⋅xi0,…,xiM，其中M+1个未知量为a0,a1,…,aM，共有M+1个方程，可以使用线性代数解决。
拟合错误与选择问题使用较高次数的多项式进行拟合，可能会导致过度拟合，使得拟合函数更接近每个数据点，因此更难以预测它们之间的关系。
另一方面，使用过低次数的多项式无法反映出数据点之间的较细节的关系。
因此，在实践中，我们需要权衡多项式次数和误差，以找到一个最合适的拟合结果。
总结最小二乘法是一种常用的曲线拟合方法，在多项式拟合中广泛应用。
通过最小化残差平方和，可以找到最优的拟合函数参数，权衡多项式次数和误差，可以得出最合适的拟合结果。
实验题目用正交多项式做小二乘曲线拟合
实验题目：用正交多项式做小二乘曲线拟合实验题目： 用正交多项式做最小二乘的曲线拟合 学生组号： 6 完成日期： 2011/11/27 1 实验目的针对给定数据的煤自燃监测数据中煤温与NO 22,之间的非线性关系，用正交多项式做最小二乘曲线拟合。
2 实验步骤2.1 算法原理设给定n+1个数据点：（yx kk,），k=0,1,···,n ，则根据这些节点作一个m 次的最小二乘拟合多项式pm（x ）=a+x a x a a mm x +++ (2)21=x a jmj j ∑=0①其中，m ≤n,一般远小于n.。
若要构造一组次数不超过ｍ的在给定点上正交的多项式函数系｛)(x Qj（ｊ＝０，１，．．．，ｍ）｝，则可以首先利用｛)(x Qj（ｊ＝０，１，．．．，ｍ）｝作为基函数作最小二乘曲线的拟合，即pm（x ）＝)(...)()(11x x x Qq Q q Q q mm+++ ②根据②式，其中的系数qj（ｊ＝０，１，．．．，ｍ）为∑∑===nk kjnk kjkjx Q x Q y q2)()(，ｊ＝０，１，．．．，ｍ ③将④代入③后展开就成一般的多项式。
构造给定点上的正交多项式)(x Qj（ｊ＝０，１，．．．，ｍ）的递推公式如下：⎪⎪⎩⎪⎪⎨⎧-=--=-==-+1,...,2,1),()()()()()(1)(11010m j x x x x x x x QQ Q Q Q j jj j j βαα ④其中αj＝dx x jk j=0，ｊ＝０，１，．．．，ｍ－１ ⑤βj＝dd j j1-，ｊ＝１，２，．．．，ｍ－１ ⑥∑==nk k jjx Q d2)(，ｊ＝０，１，．．．，ｍ－１ ⑦则实际计算过程中，根据⑤式逐步求出个正交多项式)(x Qj，并用公式④计算出q j，并将每次计算展开后累加到拟合多项式①中。
2.2 算法步骤用三个向量Ｂ,Ｔ,Ｓ,存放多项式)(1x Qj -，)(x Q j，)(1x Qj +的系数。
最小二乘拟合多项式
最小二乘拟合多项式最小二乘拟合多项式导言在数学和统计学中，最小二乘法是一种常见的数学优化和统计估计技术。
它被广泛应用于曲线拟合、参数估计和回归分析等领域。
其中，最小二乘拟合多项式是最常见和基础的应用之一。
本文将深入探讨最小二乘拟合多项式的原理、应用以及其在实际问题中的意义。
一、最小二乘法简介1.1 原理最小二乘法是一种通过最小化误差平方和来确定模型参数的方法。
在最小二乘法中，通过寻找最佳的参数估计使得模型预测值与观测值之间的差异最小化。
这样，我们可以得到一个最优的拟合曲线或函数，以便能够更好地描述观测到的数据。
1.2 应用最小二乘法在各个领域中都有广泛的应用。
在物理学中，最小二乘法常被用于拟合实验数据以确定物理定律的参数。
在工程学中，最小二乘法可用于估计信号的隐含参数，如音频信号处理中的频率分量估计。
在金融学、经济学和生物学等领域，最小二乘法也被用于回归分析、模式识别和图像处理等问题中。
二、最小二乘拟合多项式原理2.1 多项式拟合多项式拟合是最小二乘法的一种应用，用于构建一个多项式函数来拟合观测数据。
通过选择最适合的多项式次数，我们可以更好地逼近数据，并获得最优的拟合结果。
2.2 最小二乘拟合多项式最小二乘拟合多项式的目标是选择最佳的多项式来拟合给定的数据。
具体而言，它通过最小化残差平方和来确定最优的多项式系数，使得拟合曲线与观测数据之间的误差最小化。
这样，我们可以得到一个最优的拟合多项式，以便更好地描述数据的分布和趋势。
三、最小二乘拟合多项式的应用3.1 数据拟合最小二乘拟合多项式在数据拟合问题中有着广泛的应用。
通过拟合数据点，我们可以通过最小二乘法来估计数据的分布规律以及趋势。
这对于数据分析和预测具有重要意义，能够帮助我们更好地理解和利用数据。
3.2 预测与模型验证除了数据拟合，最小二乘拟合多项式还可以用于预测和模型验证。
通过构建拟合多项式，我们可以预测未来的数值或事件，并验证模型的准确性和可靠性。
正交多项式最小二乘法拟合
《MATLAB 程序设计实践》课程考核 一、编程实现以下科学计算算法，并举一例应用之。
（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年） “正交多项式最小二乘法拟合”正交多项式最小二乘法拟合原理正交多项式做最小二乘法拟合：不要求拟合函数y=f(x)经过所有点（x i ,y i )，而只要求在给定点x i 上残差δi=f(x i )-y i 按照某种标准达到最小，通常采用欧式范数||δ||2作为衡量标准。
这就是最小二乘法拟合。
根据作为给定节点x 0，x 1，…x m 及权函数ρ(x)>0,造出带权函数正交的多项式{P n （x ）}。
注意n ≤m,用递推公式表示P k （x ），即()()()()()()()01101111,,(1,2,,1)k k k k k P x P x x P x P x P x P x k n ααβ++-=⎧⎪=-⎨⎪=--=...-⎩ 这里的P k (x)是首项系数为1的k 次多项式，根据P k (x)的正交性，得()()()()()()()()()()()()()()()()()()()()2i 012i 02i 0211i 10x ,,x ,0,1,1,x ,0,1,1,x mi k i k ki k mk k k i i k k mk k k i k k i k mk k k i i x P x xP x P x a P x P x P x xP P k n P P P x P P k n P P P x ρρρβρ=+==---=⎧⎪⎪==⎪⎪⎪==⋅⋅⋅-⎨⎪⎪===⋅⋅⋅-⎪⎪⎪⎩∑∑∑∑ 根据公式(1)和(2)逐步求P k （x ）的同时，相应计算系数()()()02()(),(0,1,n (,)()mi j i k i k i kmk k iki i x x x f P a k P P x x ρϕϕρϕ=====⋅⋅⋅,∑∑）并逐步把*k a P k （x ）累加到S （x ）中去，最后就可得到所求的拟合函数曲线 ***0011n n y=S x =a P x +a P x ++a P x ⋅⋅⋅（）（）（）（）.流程图（2）（1）M文件function [p] = mypolyfit(x,y,n)%定义mypolyfit为最小二乘拟合函数%P = POLYFIT(X,Y,N)以计算以下多项式系数%P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1). if ~isequal(size(x),size(y))error('MATLAB:polyfit:XYSizeMismatch',...'X and Y vectors must be the same size.')end%检验X Y维数是否匹配x = x(:);y = y(:);if nargout > 2mu = [mean(x); std(x)];x = (x - mu(1))/mu(2);end%利用范德蒙德矩阵构造方程组系数矩阵V(:,n+1) = ones(length(x),1,class(x));for j = n:-1:1V(:,j) = x.*V(:,j+1);end% 对矩阵进行QR分解以求得多项式系数值[Q,R] = qr(V,0);ws = warning('off','all');p = R\(Q'*y);warning(ws);if size(R,2) > size(R,1)warning('MATLAB:polyfit:PolyNotUnique', ...'Polynomial is not unique; degree >= number of data points.')elseif condest(R) > 1.0e10if nargout > 2warning('MATLAB:polyfit:RepeatedPoints', ...'Polynomial is badly conditioned. Remove repeated data points.') elsewarning('MATLAB:polyfit:RepeatedPointsOrRescale', ...['Polynomial is badly conditioned. Remove repeated data points\n'...' or try centering and scaling as described in HELP POLYFIT.']) endendr = y - V*p;p = p.'; % 将多项式系数默认为行向量.5、运行流程图过程：clearx =[ 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000] y=[1.75 2.45 3.81 4.80 8.00 8.60]x1=0.5:0.05:3.0;p=mypolyfit(x,y,2)y1=p(3)+p(2)*x1+p(1)*x1.^2;plot(x,y,'*')hold onplot(x1,y1,'r')二、编程计算以下电路问题[例8-1-3]如图所示电路，已知R=5Ω，ωL=3Ω，C1ω=5Ω,Uc=100∠0,求I .R ,I .C ,I .和U .L ，U .S ，并画其相量图。
基于正交多项式拟合的气象数据处理方法
基于正交多项式拟合的气象数据处理方法＊孙宝京【摘 要】针对高空气象探测数据变化规律复杂、突变情况不可预测、数据量大等特点，提出了基于离散数据正交多项式最小二乘的分段拟合方法.首先，从计算稳定性角度，说明了采用离散数据的正交多项式最小二乘法的原因；其次，从曲线拟合的保形性角度，详细阐述了二次分段拟合原则；最后，以具有典型特征的气温探测数据为例，证明了采用基于离散数据正交多项式最小二乘的分段拟合方法能够获得明显优于传统内插方法的拟合精度，提高高空气象探测数据处理的精度和自动化处理程度.【期刊名称】吉首大学学报（自然科学版）【年(卷),期】2013(034)005【总页数】4【关键词】离散数据；正交多项式；最小二乘；分段拟合；气象探测当前气象探测站整理高空气象探测数据时，主要采用的方法是将整个探测空间分层，计算规定高度的中间数值，作为计算层数据，其他数据则采用依托计算层数据内插的方法计算.上述计算方法使用高度上的数据均由计算层数据插值得到，存在一定的截断误差和舍入误差；另外在较少的计算层数据中，必然会丢弃绝大多数探测数据，造成探测资源的极大浪费.随着气象探测及作业装备自动化水平和数据处理能力的不断提高，对高空气象探测数据处理精度和时效性提出了更高的要求，这也使得从根本上改变气象参数的计算方法成为可能［1］.由于气象要素随时间、空间在不断地变化，而且探测仪器本身也受到许多复杂因素的影响，因此各种高空气象探测数据普遍具有数据量大、变化规律复杂、突变情况较多且不可预测等特点［2］.以弹道气象诸元中最重要的高空实时气温／高度曲线为例，当探测高度达到2万m时，采用电子探空仪探测的数据超过2 000组，虽然高空实时气温总体呈下降趋势，但是在局部经常出现不可预测的无规则变化.针对高空气象探测数据特点，笔者基于离散数据正交多项式最小二乘拟合思想，提出根据数据变化率的变化情况，分段拟合各种高空气象探测数据曲线的方法.借助编程软件的自动处理，经过多次试验证明，此方法可直接得出气象通报所需要的气（虚）温－高度、气压－高度、风速－高度和风向－高度函数方程（高度为自变量），并且对高空气象探测数据处理的精确性、实效性和计算简捷性都有了很大的提高.因为高空实时气温／高度曲线在弹道气象诸元中最为重要，所以笔者将以此作为说明实例.1 离散数据分段最小二乘拟合方法用插值方法所找出的近似曲线，通过所有己知点，但是由观测或实验所获得的数据，不可避免地含有误差，这就使插值所找出的近似曲线保存所有观测误差，导致所得结果可能偏离了实际情况，出现曲线拟合的过度现象.而根据实际标准尽可能地通过节点的函数值近旁，以部分地抵消原始数据所包含的观测误差，则更具有实用价值，最小二乘拟合就是其中的一个重要方法［3］.1.1 离散数据的正交多项式最小二乘拟合在线性空间由Φ＝span（φ0，φ1，...，φn）求解最小二乘曲线拟合问题，经常会出现法方程组是病态方程组的情况，特别是当样本组数n偏大时更是如此，同时在求解法方程组时系数矩阵或右端项微小的扰动都可能导致解函数有很大的误差.为了确保计算过程的稳定性和精确性，笔者利用基于高空气象探测数据的正交多项式作为拟合函数系，使法方程组的系数矩阵变为对角矩阵，从而解得正交多项式的系数［4－5］.根据点集｛x1，x2，...，xm｝，构造出相应的正交多项式系｛p0（x），p1（x），...，pn（x）｝.由正交多项式三项递推关系，最高次项系数为1的正交多项式系｛pk（x）｝（k＝0，1，2，...，n）有如下递推关系：其中pk（x）为最高项系数为1的k次多项式.由｛pk（x）｝正交性可知，1.2 二次分段拟合原则因一次完整综合气象探测的数据量大，若对整组数据进行一次性拟合，则拟合精度和计算效果都比较差.同时，为克服传统分段曲线拟合方法中对数据点分段时经验成分较多的缺点，笔者根据数据曲线斜率的变化情况对曲线进行分段拟合，并且加入了使相邻曲线连续，即曲线边界点必须连续的约束条件［6－7］.考虑到高度是单调递增，为了适应高空气象数据曲线在局部随高度变化的无规则性，采取二次分段原则.一次分段：在气温和高度曲线中，设Δti＝（ti＋1－ti）／（hi＋1－hi），其中ti和hi分别为探测数据中第i个探测点的气温和高度，ti＋1和hi＋1分别为探测数据中第i＋1个探测点的气温和高度.若Δti－1为负，Δti，Δti＋1，Δti＋2，...，Δti＋a（a≥1）为正，Δti＋a＋1，...，Δti＋a＋b（b≥1）为负，Δti＋a＋b＋1为正，则｛ti，...，ti＋a＋b｝为第i个分段，第i＋1个分段从点ti＋a＋b开始.二次分段：在气温和高度曲线中，第i个分段已经由一次分段确定，设δi＋c＝Δti＋c－Δti＋c－1（1≤c≤a＋b－1）.若δi＋1，δi＋2，...，δi＋d（d≥1）均大于标准值δ0（为负数），δi＋d＋1≤δ0，且δi＋d＋2，...，δi＋d＋e（e≥1）小于或等于（－δ0），δi＋d＋e＋1＞（－δ0），则ti＋d＋e＋1为分段点，即｛ti，...，ti＋d＋e＋1｝为一个分段，下一个分段以ti＋d＋e＋1为起点，从δi＋d＋e＋2开始判断；若δi＋1≤δ0，且δi＋2，...，δi＋e（e≥1）小于或等于（－δ0），δi＋e＋1＞（－δ0），则ti＋e＋1为分段点，即｛ti，...，ti＋e＋1｝为一个分段，下一个分段以ti＋e＋1为起点，从δi＋e＋2开始判断［8－9］.为使在单个分段内拟合曲线在误差允许范围内最大程度地保持数据曲线的形状，并且减少数据处理程序的计算量，实际计算过程对原始探测数据中心化计算后，采用如下做法进行拟合：若a＋b≥5，则用四阶正交多项式为拟合函数系，即当一个分段包含4个以上的点时用四阶正交多项式拟合，若a＋b＝4，则用三阶正交多项式为拟合函数系，即当一个分段包含4个点时用三阶正交多项式拟合，若a＋b＝3，则用二阶正交多项式为拟合函数系，即当一个分段包含3个点时用二阶正交多项式拟合，同时，使参与第i个分段拟合的第i个分段中的最后一个数据点，作为第i＋1个分段的第1个数据点参与第i＋1个分段的拟合，这样第i个分段的拟合函数和第i＋1个分段的拟合函数在分界点上的值相等，从而实现所有分段拟合函数在整体上保持连续［10］，即2 计算实例以编号为087233的GZZ8型电子探空仪探测数据为例.图1是该组探测数据曲线中一部分，温度t∈［－61.31，－62.27］，高度H∈［15 125.89，1 5 494.56］.根据文中的分段规则，数据处理系统自动将此段数据曲线分为5段进行拟合.第1个分段为前13个点，h∈［15 125.89，15 230.23］，拟合函数为平方误差为0.08.第2个分段为后续7个点，h∈［15 230.23，15 289.47］，拟合函数为平方误差为0.03.第3个分段为后续10个点，h∈［15 289.47，15 397.85］，拟合函数为平方误差为0.02.第4个分段为后续5个点，h∈［15 397.85，15 439.32］，拟合函数为平方误差为0.01.第5个分段为最后7个点，h∈［15 439.32，15 494.56］，拟合函数为平方误差为0.03.根据传统高空气象探测数据处理方法，14 000～16 000m高度范围内的相关气象数据从标准高度层14 000m和16 000m处的实际探测数据进行线性内插得到.在此例中，14 000m和16 000m处的气温分别为－56.16℃和－62.37℃，从而内插得到高度为15 142.00，15 276.4，15 378.67，15 417.06，1 5 479.1m处的气温.各高度对应的实际探测值、内插结果和拟合结果见表1.由表1可见，采用基于离散数据正交多项式最小二乘的分段拟合法，可以获得比传统内插法更好的拟合精度.3 结语针对高空气象探测数据量大、变化规律复杂、突变情况较多且不可预测等特点，提出了基于离散数据正交多项式最小二乘的分段拟合法.经过实践证明，该方法能够根据算法自动分段，拟合出符合精度要求的曲线，在处理程序自动计算过程中确保计算的稳定性和高效性，避免传统高空气象探测数据处理方法中不断累计的舍入误差对结果的影响，提高高空气象探测数据处理结果的精度和自动化处理程度.参考文献：［1］ 孙宝京.炮兵防空兵气象信息系统研究［D］.沈阳：东北大学，2004.［2］ 曲延禄.外弹道气象学概论［M］.北京：气象出版社，1987：59－176.［3］ FRANK R.Scattered Data Interpolation：Tests of Some Methods［J］.Mathematical Computation，1982，38：181－200.［4］ MICCHELLI C.Interpolation of Scattered Data：Distance Matrix and Conditionally Positive Definite Functions［J］.Constructive Approximation，1986，2：11－22.［5］ 蔺小林，蒋耀林.现代数值分析［M］.北京：国防工业出版社，2004：212－255.［6］ WU Zong-min，SCHABACK R.Shape Preserving Properties and Convergence of Univariate Multi-Quadric Quasti-Interpolation［J］.ACTA Mathematice Application Sinica，1994，10：441－446.［7］ GREGORY J A.Shape Preserving Spline Interpolation［J］.Computer-Aided Design，1986，18（1）：53－57.［8］ KAUFMANN E，KLASS R.Smoothing Surfaces Using Reflection Lines for Families of Splines［J］.Computer-Aided Design，1988，20（6）：312－316.［9］ HOSCHEK J.Smoothing Curves and Surfaces［J］.Computer-Aided Geometric Design，1985，2（1－3）：97－105.［10］ JORG PRTERS.Local Smooth Surface Interpolation：A Classification［J］.Computer-Aided Geometric Design，1 990（7）：191－195.（责任编辑 向阳洁）。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《MATLAB 程序设计实践》课程考核一、编程实现以下科学计算算法，并举一例应用之。
（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“正交多项式最小二乘法拟合”正交多项式最小二乘法拟合原理正交多项式做最小二乘法拟合：不要求拟合函数y=f(x)经过所有点（x i,y i)，而只要求在给定点x i上残差δi=f(x i)-y i按照某种标准达到最小，通常采用欧式范数||δ||2作为衡量标准。
这就是最小二乘法拟合。
根据作为给定节点x0，x1，…x m及权函数ρ(x)>0,造出带权函数正交的多项式{P n（x）}。
注意n≤m,用递推公式表示P k（x），即()()()()()()()1101111,,(1,2,,1) k k k k kP xP x x P xP x P x P x k nααβ++-=⎧⎪=-⎨⎪=--=...-⎩这里的P k(x)是首项系数为1的k次多项式，根据P k(x)的正交性，得()()()()()()()()()()()()()()()()()()()()2i12i2i211i1x,,x,0,1,1,x,0,1,1,xmi k ik kik mk kk ii k kmk kk ik kik mk kk iix P xxP x P xaP x P xP xxP Pk nP PP xP Pk nP PP xρρρβρ=+==---=⎧⎪⎪==⎪⎪⎪==⋅⋅⋅-⎨⎪⎪===⋅⋅⋅-⎪⎪⎪⎩∑∑∑∑根据公式(1)和(2)逐步求P k（x）的同时，相应计算系数()()()2()(),(0,1,n(,)()mi j i k ik ik mk ki k iix x xf Pa kP Px xρϕϕρϕ=====⋅⋅⋅,∑∑）并逐步把*ka P k（x）累加到S（x）中去，最后就可得到所求的拟合函数曲线***0011n ny=S x=a P x+a P x++a P x⋅⋅⋅（）（）（）（）.流程图（2）（1）M文件function [p] = mypolyfit(x,y,n)%定义mypolyfit为最小二乘拟合函数%P = POLYFIT(X,Y,N)以计算以下多项式系数%P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1).if ~isequal(size(x),size(y))error('MATLAB:polyfit:XYSizeMismatch',...'X and Y vectors must be the same size.') end%检验X Y维数是否匹配x = x(:);y = y(:);if nargout > 2mu = [mean(x); std(x)];x = (x - mu(1))/mu(2);end%利用范德蒙德矩阵构造方程组系数矩阵V(:,n+1) = ones(length(x),1,class(x));for j = n:-1:1V(:,j) = x.*V(:,j+1);end% 对矩阵进行QR分解以求得多项式系数值[Q,R] = qr(V,0);ws = warning('off','all');p = R\(Q'*y);warning(ws);if size(R,2) > size(R,1)warning('MATLAB:polyfit:PolyNotUnique', ...'Polynomial is not unique; degree >= number of data points.') elseif condest(R) > 1.0e10if nargout > 2warning('MATLAB:polyfit:RepeatedPoints', ...'Polynomial is badly conditioned. Remove repeated data points.')elsewarning('MATLAB:polyfit:RepeatedPointsOrRescale', ...['Polynomial is badly conditioned. Remove repeated data points\n'...' or try centering and scaling as described in HELP POLYFIT.'])endendr = y - V*p;p = p.'; % 将多项式系数默认为行向量.5、运行流程图过程：clearx =[ 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000] y=[1.75 2.45 3.81 4.80 8.00 8.60]x1=0.5:0.05:3.0;p=mypolyfit(x,y,2)y1=p(3)+p(2)*x1+p(1)*x1.^2;plot(x,y,'*')hold onplot(x1,y1,'r')二、编程计算以下电路问题[例8-1-3]如图所示电路，已知R=5Ω，ωL=3Ω，C1ω=5Ω,Uc=100∠0,求I .R ,I .C ,I .和U .L ，U .S ，并画其相量图。
理论分析：根据电路分析Z=R+j*(Xl-Xc)Ic=Uc/Z3;Z3=-j*XcIr=Ur/Z2=Uc/Z2;Z2=RI=Ir+IcUl=I*Z1;Z1=j*XLUs=Ul+Ur计算得Ir =2；Ic =2.00iI =2.00 + 2.00iUl =-6.00 + 6.00iUs =4.00 + 6.00iM文件clearR=5;XL=3;XC=5;UC=10;UR=UC;%为给定元件赋值Z1=j*XL;Z2=R;Z3=-j*XC;%定义各电抗disp('电流')IR=UR/Z2%计算IRIC=UC/Z3%计算ICI=IR+IC%计算Idisp('电压')UL=I*Z1%计算ULUS=UC+UL%计算USdisp('IR IC I UL US')disp('幅值');disp(abs([IR,IC,I,UL,US]))disp('相角');disp(angle([IR,IC,I,UL,US])*180/pi)ph=compass([IR,IC,I,UL,US]);%分别画出IR,IC,I,UL,US相量图set(ph,'linewidth',3)运行流程图：运行图开始读取已知数据构造电抗Z1，Z2，Z3计算Ir ,Ic, I ,Ul ,UsIR=UR/Z2；IC=UC/Z3；I=IR+IC；UL=I*Z1；US=UC+UL结束输出Ir ,Ic, I ,Ul ,Us并画出对应的相量图三、编程解决以下问题求自然三次样条曲线，经过点（-3,2），（-2,0），（1,3），（4,1），而且自由边界条件S ''(-3)=0,S ''(4)=0。
算法说明：三次样条也是分片三次插值函数，它是在给定的区间[a ，b]上的一个划分：a=x 0<x 1<…<x n =b,已知函数f(x)在xj 上的函数值为 f(x j )=y j ，（j=0,1,2,3...,n)如果存在分段函数[][][]⎪⎪⎩⎪⎪⎨⎧∈⋯⋯∈∈=-n n n x x x x S x x x x S x x x x S x S ,)(,)(,)()(1212101满足下述条件：（1）S(x)在每一个子区间 [ x j-1，x j ]，（j=0,1,2,3...,n)上是三次多项式； （2）S(x)在每一个内接点 （j=0,1,2,3...,n)具有直到二阶的连续导数； 则成为节点x 0，x 1，…，x n 上的三次样条函数。
若 S(x) 在节点 x 0，x 1，…，x n 上海满足插值条件：（3）S(x j)=y j （j=0,1,2,3...,n)则称S(x)为三次样条插值函数。
由（1）知，S(x)在每一个小区间[ x j-1，x j ]上是一三次多项式，若记为S j(x)，则可设：S j(x)=a j x3+b j x2+c j x+d j要确定函数S(x)的表达式，需确定4n个未知数{aj,bj ,cj, dj}（j=0,1,2,3...,n）。
由（2）知S(x)，S'(x)，S''(x)在内节点x1, x2,....,x n-1上连续，则Ⅰ型S(x j-0)=S(x j+0)Ⅱ型S'(x j-0)=S'(x j+0)Ⅲ型S''(x j-0)=S''(x j+0) j=0,1,2,3...,n-1可得3n-2个方程，又由条件（3）S(x j)=y j j=0,1,2,3...,n得n个方程，共得到4n-2个方程。
要确定4n个未知数，还差俩个方程。
通常在端点x0=a. x n=b处各附加一个条件称边界条件，常见有三种：（1）自然边界条件：S''(x0)=S''(x n)=0（2）固定边界条件：S'(x0)=f'(x0),S'(x)=f'(x n)（3）周期边界条件：S'(x0)=S'(x n)，S''(x0)=S''(x n)共4n个方程，可唯一的确定4n个未知数流程图：（2）运行流程图源代码：function s=myspline2(x0,y0,y21,y2n,x)%s=myspline(x0,y0,y21,y2n,x)%x0,y0 为已知插值点,y21,y2n为二型边界条件n=length(x0);km=length(x);a(1)=-0.5;b(1)=3*(y0(2)-y0(1))/(2*(x0(2)-x0(1)));for j=1:(n-1)h(j)=x0(j+1)-x0(j);%h(j)为第j个插值区间的长度endfor j=2:(n-1)alpha(j)=h(j-1)/(h(j-1)+h(j));beta(j)=3*((1-alpha(j))*(y0(j)-y0(j-1))/h(j-1)+alpha(j)*(y0(j+1)-y0(j))/h(j));a(j)=-alpha(j)/(2+(1-alpha(j))*a(j-1));b(j)=(beta(j)-(1-alpha(j))*b(j-1))/(2+(1-alpha(j))*a(j-1)); end%alpha(j),beta(j)为插值基函数m(n)=(3*(y0(n)-y0(n-1))/h(n-1)+y2n*h(n-1)/2-b(n-1))/(2+a(n-1)); for j=(n-1):-1:1m(j)=a(j)*m(j+1)+b(j);endfor k=1:kmfor j=1:(n-1)if ((x(k)>=x0(j))&(x(k)<x0(j+1)))l(k)=j;endendendfor k=1:kmsum=(3*(x0(l(k)+1)-x(k))^2/h(l(k))^2-2*(x0(l(k)+1)-x(k))^3/h(l(k) )^3)*y0(l(k));sum=sum+(3*(x(k)-x0(l(k)))^2/h(l(k))^2-2*(x(k)-x0(l(k)))^3/h(l(k) )^3)*y0(l(k)+1);sum=sum+h(l(k))*((x0(l(k)+1)-x(k))^2/h(l(k))^2-(x0(l(k)+1)-x(k))^ 3/h(l(k))^3)*m(l(k));s(k)=sum-h(l(k))*((x(k)-x0(l(k)))^2/h(l(k))^2-(x(k)-x0(l(k)))^3/h (l(k))^3)*m(l(k)+1);end举例：x=[-3 -2 1 4]y=[2 0 3 1]x0=[-3:0.15:4];y0=myspline(x,y,0,0,x0);plot(x0,y0)hold onplot(x,y,'or')legend('计算值','实验值')。