最小二乘建模

最小二乘建模（共7篇）

最小二乘建模 篇1

0 引言

地壳形变与各种精密工程变形监测,不论垂直位移监测,还是水平位移监测,常做定点、定期重复观测,要处理的数据是一组时间序列,处理目的是确定位移量与时间推移的动态变化规律,发现变形发生的原因、变形特征,以便做出科学预测、预报[1,2]。

时间序列方法分析定点复测数据的成果有很多,其讨论的问题主要集中在根据定点实测数据选用AR(p)自回归模型、MA(q)可逆滑动平均模型或者ARMA(p,q)自回归—可逆滑动平均混合模型;模型的识别方法;确定模型的阶数以及模型参数的估计等。事实上,由于变形机理的复杂性,定点观测数据随时间变化过程中会包含各种因素综合的系统误差和随机误差的干扰,使所建立的AR模型、MA模型或者ARMA模型只是实际问题的近似表达,不可避免地附加了模型误差,误差的存在势必影响所建模型及预测结果的准确性[2]。

近年来,研究人员从信号处理的角度出发,将定点观测数据分解成趋势项、周期项和随机项三部分,根据不同部分的特点,分别采用不同的模型进行模拟和逼近,取得了一系列的成果[3,4]。然而,趋势项与周期项相互影响,如果趋势有误差,对于提取周期项和随机项会造成较大的影响。因此,为了更全面地反映变形特性,有效滤除模型误差,本文以AR(p)模型为例,考虑模型误差对变形预测的影响,提出半参数补偿最小二乘精化AR模型,将模型误差用非参数项加以弥补,并将其运用于形变时间数据序列分析中。

1 顾及模型误差S的AR(p)建模

目前广泛采用的AR(p)模型为:

式(1)中:xt为观测数据序列,φi为自回归系数,at为白噪声序列。式(1)的平差模型可以写成:

L为n维观测向量;B为列满秩设计矩阵;x为t维参数向量(对应式(1)中φi);t为必要观测数;Δ为n维观测误差向量。此模型中,通常假定Δ是期望为0的偶然误差,也就是说除去观测误差,观测值L完全表示为未知参数x的函数。

考虑模型误差,设AR(p)的模型误差为S,当模型误差S与偶然误差Δ相比是一个微小量时,忽略模型误差S不会对参数估值产生太大的影响,即可采用式(1)对形变时间序列进行模拟;而当S比较大时,就会对参数的估计产生较大的影响(这一点可以从本文后面的计算结果中得到验证)。此时,式(1)、式(2)并不能严格成立,式(1)应改写为:

式(2)改写为:

式(3)即为考虑模型误差条件下的AR(p)模型,S=[s1,s2,Λsn]是一个描述模型误差的n维未知向量。考虑一般的情形,可认为模型误差或观测值的系统误差的性态非常复杂,无法用少数参数表示,因此给每个观测方程增加一个待定量,也就是所谓的非参数分量。这样在观测方程中既有参数分量又有非参数分量,因此式(4)称为半参数模型[2,3,5]。

2 补偿最小二乘估计下AR(p)的解

据式(4)写出误差方程:

式(5)的求解是在最小二乘的目标函数基础上增加一个非参数部分的补偿项,即求解如下条件极值问题:

式(6)即为补偿最小二乘平差准则。α是一个给定的纯量因子,在极小化过程中对V和S起平衡作用,称为平滑因子(平滑参数)。R是一个适当给定的正定矩阵,又称为正规化矩阵,是被用来估计参数S的某种函数关系类型,二次型STRS反映对向量S的某种度量。根据文献[5,6,7,8,9,10],得到补偿最小二乘模型的解:

其中N=BTPB,M=P+αR-PBN-1BTP。

这样就可通过式(7)、式(8)及式(5)计算AR(p)模型参数分量的估值、非参数分量的估值S及观测值改正数V。通过分析得到的非参数分量,就可以重新认识所选的数学模型,从而实现对模型的精化。补偿最小二乘法的关键是选择适当的正规化矩阵R和平滑参数α,R和α是事先给定的量[5,6,7,8]。本文计算中,根据变形数据时间序列的特点,取相邻两观测点模型误差之差的平方和来确定R,即R=TTT,α的确定采用作者提出的一种Xu函数法。

3 Xu函数确定平滑参数α

补偿最小二乘估计中,选择合适的R后,需确定平滑参数α,已有的研究文献可以通过信噪比值法、L-曲线法以及交叉核实法和广义交叉核实法等[6,7,8,9,10]来实现。作者在比较经典最小二乘与补偿最小二乘估计基础上,提出一种新的α确定方法,Xu函数(关于Xu函数的有效性及与其它确定α的方法比较见文献[10])。

偶然误差与模型误差相比是一个非常小的量,经典最小二乘平差中,不考虑模型误差的存在,事实上模型误差与偶然误差都看成误差来处理,补偿最小二乘估计则将模型误差看成一个独立变量来处理,如果我们选择了合适的正规化矩阵R和平滑参数α,那么补偿最小二乘求出来的模型误差与经典平差理论求出来的误差之间的差值应该是个很小的值,据此得出Xu函数:

式(10)中:e为经典平差求得的误差,R为正规化矩阵,S为求得的模型误差,V为补偿最小二乘求得的误差,P为观测量的权阵。

4 补偿最小二乘估计下AR(p)建模步骤

(1)根据定点变形观测数据,按式(1)建立AR模型,依据动态建模思想和F统计检验确定模型阶数p;

(2)按式(3)建立顾及模型误差条件下的AR模型,由补偿最小二乘估计理论确定AR(p)模型的参数分量φi,以及非参数分量(模型误差)的估值S;

(3)根据(2)中求出的参数值,可建立预测方程进行形变监测预报;预报方程为式(3)。

5 算例分析

为了说明补偿最小二乘估计AR(p)模型在变形数据处理中的应用及其处理模型误差的可行性,同时比较本文方法与现有方法的优劣性,利用某建筑物的36次沉降观测序列[11]进行建模计算,采用如下方案:

方案1:利用1～30期观测数据建立补偿最小二乘AR(p)模型;方案2:利用1～30期观测数据建立AR(p)模型;方案3:利用1～30期观测数据建立GM(1,1)模型;方案4:利用灰色–神经网络(GM+BP)[4]模型。

当考虑模型误差影响时,用补偿最小二乘估计解算AR(p)模型,采用Xu函数确定平滑参数α,α=0.001,其中参数分量φi为φ=(-0.924967,-0.589398,0.040616),模型误差估值见图1。当不考虑模型误差时,根据最小二乘解算的AR(p)模型参数φ'i,φ'=(0.041086,0.327809,0.635059),AR(p)模型、GM(1,1)模型及(GM+BP)模型计算结果见图2、图3。

比较补偿最小二乘估计下的AR(p)模型与经典最小二乘下参数AR(p)模型(图1、图2),发现半参数AR(p)模型得到的观测值平差值与其真值的拟合精度明显优于常规AR模型,其特点是能够在补偿最小二乘准则下分离出平差模型中的模型误差,削弱了对预测结果的影响,从而获取比较可(上接第53页)靠的参数估值。

从图2可以看出GM(1,1)模型计算拟合曲线,其只反应了变形观测数据序列的均值变化,相当于实测曲线的平均光滑曲线,究其原因,该模型建立在指数回归模型基础上,建议使用该模型时应结合形变观测数据的实际情况。

采用文献[4]时间序列的分解方法进行变形预测(图3),考虑变形数据的特点,在拟合预测精度方面比单一模型有所提高,但其精度明显低于半参数AR(p)模型。

6 结论

(1)由算例分析可知,对于形变时间序列数据建模,应考虑模型误差对预测结果的影响,函数模型存在模型误差时,经典最小二乘平差很难发现和识别,若忽略模型误差,将给参数估值带来不利影响;采用补偿最小二乘估计解算,同时兼顾了参数和非参数两类因素,精度有大幅度提高,可以达到毫米级,说明补偿最小二乘方法应用于形变时间序列分析是有效的,同时说明本文提出的Xu模型确定平滑参数α是一种合理的方法。

(2)在补偿最小二乘模型中,正则矩阵R和平滑参数α选取的合适与否将直接影响模型的精度,由于篇幅有限,本文计算过程中并未讨论选取不同正则矩阵R和平滑参数α对变形预测精度的影响,需进一步深入研究。

最小二乘建模 篇2

针对牛粪高温厌氧发酵过程的特点,本文主要研究了将最小二乘支持向量机(Least Squares Support Vector Machines,LS—SVM)这种基于支持向量机(Support Vector Machine,SVM)改进后的快速算法,应用到牛粪高温厌氧发酵这一生化过程中来。基于实验数据的仿真结果表明,该方法能够通过多批次的样本数据的学习,建立沼气产量的在线预估模型。

1 最小二乘支持向量机建模方法

众所周知,最小二乘法(Least Squares,LS)是解决多元函数回归的经典方法。Suykens等人[3,4]在SVM的基础上引进了LS,将不等式约束改为等式约束,避免了求解耗时的QP问题,并将求解优化问题转化为求解线性方程,从而很大程度上降低了运算时间,为在线估计创造了有利条件。而且相对于常用的线性不敏感损失函数,LS-SVM不再需要指定逼近精度ε,只是LS-SVM的解不具有稀疏性。

用于函数估计最小二乘支持向量机的算法推导如下:

设undefined,其中xk∈Rn为输入数据,yk∈R为输出类比。在权w空间中最小二乘支持向量机分类问题可以描述如下:

undefined

约束条件:yk[wTφ(xk)+b]=1-ekk=1,…,N 。

根据式(1),定义拉格朗日函数:

undefined

其中,拉格朗日乘子αk∈R。对上式进行优化,即对w,b,ek,αk的偏导数等于0,可以得到矩阵方程:

undefined

其中,undefined。

同时将Mercer条件代入到Ω=ZZT,可得:

Ωkl=ykylφ(xk)Tφ(xl)=ykylΨ(xk,xl) (4)

因此,式(1)的分解可以通过解式(4)和式(5)获得最小二乘支持向量机函数为:

undefined

其中:Ψ(x,xk)是核函数,目的是从原始空间抽取特征,将原始空间中的样本映射为高维特征空间中的一个向量,以解决原始空间中线性不可分的问题。

2 厌氧发酵沼气产量预估模型

本文用于建模分析的数据来自牛粪高温厌氧发酵工程,采集发酵温度、PH值、固体浓度、进料量和出料量等数据,将其中一部分样本数据作为训练集,剩余的样本作为测试集。所有的仿真计算实验在1台Pentium850MHz内存2GB的计算机上进行,编程语言使用Matlab7.0。

2.1 参数调整与分析

LS-SVM方法中有两个重要的参数γ和σ2需要调整[5]。γ为正则化参数即惩罚系数,用于控制函数的拟合误差,γ值越大,拟合误差越小,相应的训练时间也就越多,但γ过大会导致过拟合;σ2是核函数参数,代表着带宽,随着σ2的变小,拟合误差会变小,相应的训练时间也就变长,但σ2过小会导致过拟合[6]。目前关于对参数的选择还没有一套系统有效的方法,本文以牛粪高温厌氧发酵这一具体过程,对拟合误差、预测误差及训练时间随γ和σ2的变化趋势进行了仿真,仿真实验结果如图1、2所示。

从图1、2中我们可以很清楚地得到其基本方法,首先将σ2固定,让γ不断增大,因为随着γ的增大,拟合误差和预测误差都会先变小,当预测误差首次出现极小值时即停止训练,相应的γ值作为选定值。然后将此γ值作为初值,用同样的方法训练得到σ2,于是LS-SVM的参数(γ,σ2)即可确定。这样得到的(γ,σ2)不仅保证了预测误差小即模型的推广能力好,而且训练时间短,可以适应实时在线预报的要求。

2.2 预估结果与分析

选取4批次数据,对其中2-4批次的样本数据分别建立沼气产量在线预报模型。当一共有i批次的样本时,取前i- 1批次进行训练并计算拟合误差,第i批次用来检验模型的推广能力并计算预测误差。表1列出了所建预报模型的预测误差。

从表1中可以看出预报模型的精度很高,而且通过批次的增加预测误差还将慢慢变小,说明模型还有一定的自学习能力。为了全面检验模型的推广能力,再采用交叉验证法进行仿真,即任取4批数据中的3批按上述方式训练预报模型,第4批次作为检验,分别计算每次的拟合误差、预测误差和训练时间,列于表2。

从表2得知虽然训练样本改变,但所建预报模型的预测误差变化并不大,说明此方法有一定的鲁棒性,而且训练时间都很短,可以对相关变量进行实时预报。图3列出了对样本数据的在线预报沼气产量模型和实际生产数据的比较结果。

由图3可知,所建立的预报模型有较高的精度,这表明在发酵过程操作条件下,仅通过几个输入变量就可以实现对沼气产量进行准确的在线预报。

需要指出的是,实际生产中一般是无法得到沼气产量的准确估计的,必须考虑更多的过程因素。但本文所提出的方法仍旧适用,在实际应用中,只需将输入实际工况所要估计的相关数据,而且可以选取多个可测变量作为LS-SVM的输入来实现对沼气产量的估计。

3 结论

牛粪发酵过程是一个复杂的生化工程,建立其沼气产量精确预测模型对生化过程的优化控制具有十分重要的作用。本文采用LS-SVM建模方法,建立了牛粪厌氧发酵过程变量在线预报模型,在对数据的要求方面远远低于神经网络和多元回归等传统方法。通过仿真实验获得了LS-SVM参数调整策略。仿真结果表明此模型具有很强的拟合、泛化能力和学习能力,对样本依赖程度低,可以较合适地描述时变和非线性的特性,适合复杂的发酵预估,具有广泛的应用前景。本研究是具有实用价值的预测方法,进一步指导大型沼气工程的智能控制,是我们下一步研究的工作。

参考文献

[1]张智焕.复杂系统预测控制算法及其应用研究[D].浙江大学,2002:17-20.

[2]王建林,于涛.发酵过程生物量软测量技术的研究进展[J].现代化工,2005,25(6):22-25.

[3]Xiong ZH,Zhang J.Neural network model based on-line re-op-timization control of fed-hatch processes using a modified itera-tive dynamic programming algorithm[J].Chemical Engineering and Processing,2005,44(4):477-484.

[4]Hamrita T.K and Wang S.Patten recognition for modeling and on-line diagnosis of bioprocesses,IEEE Trails,On Industry Applica-tions,2000,36(5),1295-1299.

[5]刘国海,周大为,徐海霞,等.基于SVM的微生物发酵过程软测量建模研究[J].仪器仪表学报,2009,30(6):1228-1232.

最小二乘建模 篇3

在实验数据的处理、分析时, 数据拟合是经常采用的一种方法。数据拟合的目标是找到能反映变量之间关系的一种表达式, 使其在某种准则下最佳地接近已知数据。其原理有最小二乘法、契比雪夫法等, 且以最小二乘法最为常见和重要。

随着人类认识能力的不断进步以及计算技术的快速发展, 对于变量之间的未知关系, 应用曲线拟合的方法揭示其内在规律具有重要的理论和现实意义。针对最小二乘曲线拟合的有关问题以及相应的MATLAB实现进行探讨。

1 最小二乘曲线拟合

曲线拟合是指:已知n个数据点 (xi, yi) , i=1, 2, ……, n, 其中xi不全相同, 寻求函数f (x;a1, a2, …, am) 的待定参数a1, a2, …, am的一组取值, 使得在这组取值之下, 函数f (x;a1, a2, …, am) 与已知n个数据点整体上最为接近。

最小二乘曲线拟合方法根据已知数据, 首先构造出能够反映含有待定参数的函数f (x;a1, a2, …, am) 与n个数据点 (xi, yi) , i=1, 2……, n偏离程度的函数:

然后应用数学方法求函数J (a1, a2, …, am) 的最小值a1, am2, i…n, am=J (a1, a2, …, am) , 此时a1, a2, …, am的取值就是所求的待定值。这样一组取值使得函数f (x;a1, a2, …, am) 与n个数据点在二次平方和意义下最为接近。

2 最小二乘曲线拟合的MATLAB实现

2.1 线性最小二乘曲线拟合的MATLAB实现

线性最小二乘曲线拟合的含有待定参数的函数形式为:

其中rk (x) 为事先选定的一组已知函数, ak为待定系数, k=1, 2, …, n, m

该方程组可简化为RTRA=RTY, 其中:

当r1 (x) , …, rm (x) 线性无关时, RTR可逆, 方程组RTRA=RTY有唯一解。

根据线性最小二乘拟合理论可得关于待求参数构成的矩阵A的解, 利用MATLAB下命令A=RY, 可以直接求得待求参数, 下面举例说明。

例1下表给出了一组实测数据

已知数据 (xi, yi) 的函数原型为:

试用已知数据进行曲线拟合, 求出待定参数的值。

在MATLAB中输入以下程序:

运行程序, 输出待定参数结果为:

下面画出拟合得到的曲线及已知的数据散点图:

由图1可见, 曲线拟合效果非常好。事实上, 所给数据就是由已知曲线:

2.2 非线性最小二乘曲线拟合

在最小二乘曲线拟合中, 如果待定参数不能全部以线性形式出现, 如指数拟合函数f (x) =a0+a1e-a2·x, 这便是非线性最小二乘曲线拟合。

非线性最小二乘曲线拟合与线性最小二乘曲线拟合的原理没有什么区别, 但是最小二乘的解常常难以通过人的手算实现, 从而制约了该方法的应用。随着计算机技术的进步、专业软件的不断涌现, 这一问题的求解已不再困难。但是, 非线性曲线拟合中初值的选取是一个重要的问题, 目前为止还没有固定的理论或方法给出一般性的结论。

在MATLAB中实现非线性最小二乘曲线拟合有三个常用的命令。

(1) lsqcurvefit () 命令, 其使用格式为x=lsqcurvefit (fun, x0, xdata, ydata) , 其中fun是要拟合的非线性函数, x0是初始参数, xdata, ydata是拟合点的数据, 该函数最终返回系数矩阵。

(2) nlinfit () 命令, 其应用格式为beta=nlinfit (x, y, fun, beta0) , 其中x和y是拟合点数据, fun是被回归 (拟合) 的函数, beat0是初始参数。

(3) lsqnonlin () 命令, 其应用格式为x=lsqnonlin (fun, x0) , 其中fun的定义与lsqnonlin () 函数中fun的定义有差别, x0仍为初始参数向量, 将输出的系数结果放在变量x中。

例2假设已知

并已知该函数原型为y (x) =a1exp (a2x) +0.54exp (a3x) .sin (a4x) , 其中ai为待定系数。

在MATLAB中输入以下命令:

%建立函数原型, 则可以根据他来进行下面的求取系数的计算

运行结果为:

所求得的系数与原式中的系数相近。

如果想进一步提高精度, 则需修改最优化的选项, 函数的调用格式也将随之改变。

在MATLAB中输入以下命令:

%修改精度, 即是修改其限制条件

%两个空矩阵表示系数向量的上下限a', res运行结果为:

下面画出由拟合得到的曲线及已知的数据散点图, 如图2。

3 结论

首先介绍了最小二乘法的原理, 其次对线性最小二乘曲线拟合和非线性最小二乘曲线拟合应用MATLAB进行了具体实现, 使得相关曲线拟合理论更加生动易懂, 这使我们可以在今后的研究和工作中建立曲线模型对相应对象进行曲线拟合, 并应用MATLAB来实现, 从而找到更好更形象的反映变量之间关系的曲线, 也就找到最合适的模型了。

参考文献

[1]薛定宇, 陈阳泉.高等应用数学问题的MATLAB求解[M].北京:清华大学出版社, 2004.

[2]刘琼荪, 何中市, 傅鹂等.数学实验[M].北京:高等教育出版社, 2000.

[3]梁国业, 廖健平.数学建模[M].北京:冶金工业出版社, 2004.

轴承沟道形状误差的最小二乘评定 篇4

轴承主要用于确定旋转件与固定件相对运动,是起支承或导向作用的典型零部件,在机器及各类机械中占有重要地位。因此,研究轴承沟道形状误差评定对提高轴承的精度及整个机器或机械的性能有重要的意义。目前多使用圆度、线轮廓度来评价轴承沟道形状误差,即用一条线轮廓和底圆形成的几何要素的误差对轴承沟道的形状误差进行评定,而对轴承沟道整体形状误差的评定国内外很少报道。本文依据最小二乘形状误差的定义,提出了轴承沟道形状误差的最小二乘评定。

1 最小二乘算法过程及步骤

1.1 提取测量点坐标

将工件水平放置在工作台上,利用三坐标测量机对轴承沟道的每条线轮廓进行采样取点,从中取得的测量点为Pij(xij,yij,zij)(i表示测量线轮廓数,i=1,2,…,M;j表示每条线轮廓上的测量点数,j=1,2,…,N)。

1.2 求空间每条线轮廓的中心点

通过三坐标测量机对轴承沟道每条线轮廓进行测量并得到采样点,用最小二乘法对每条线轮廓求其半径及圆心的位置,但轴承沟道轮廓是非整圆,不能直接用整圆的最小二乘法公式,需要坐标变换求空间每条线轮廓的中心点。

1.2.1 坐标变换

在空间坐标系中,直接用最小二乘法求每条线轮廓的中心点比较复杂,为了便于计算,将空间内的线轮廓绕z轴转θi角度至yoz平面。其中:

undefined。

空间线轮廓上各点坐标通过坐标(xij,yij,zij)变换至yoz平面上的新坐标为:undefined。

1.2.2 求所有平面线轮廓的圆心

设平面线轮廓的最小二乘拟合方程为:

(y-a′i)2+(z-b′i)2=rundefined。 (1)

其圆心为(a′i,b′i),半径为ri,平面线轮廓最小二乘拟合方程见图1。对式(1)展开并整理得:

y2+z2-2ai′y-2bi′z+(ai′2+bi′2-rundefined)=0 。 (2)

令c=ai′2+bi′2-rundefined,由式(2)得:

y2+z2-2ai′y-2bi′z+c=0 。 (3)

以每条线轮廓上测点坐标(y′ij,z′ij)代入式(3),方程式明显不等于零,而为:

Δij=y′2ij+z′ij2-2a′iy′ij-2b′iz′ij+c 。 (4)

其中:Δij为实际线轮廓上各点与理论线轮廓上对应点的函数偏差值。设:

当偏差值为最小值时,理论线轮廓逼近实际线轮廓效果最佳。令:

依据求最小值的方法,对系数a′i、b′i、c求偏导数后得出:

其中:undefined。

再由式(3)得出:undefined,在求出待定系数a′i,b′i,c和ri后,可以得出平面线轮廓上所有中心点O′i(a′i,b′i)。

1.2.3 坐标变换

将上一步得到的平面线轮廓的中心点代入下式可以得到空间线轮廓上的中心点Oi(ai,bi,ci):

1.3 最小二乘法拟合空间圆

将每条线轮廓的中心点用最小二乘法拟合一个空间圆,但空间圆没有特定的方程, 只能由空间圆和空间平面的方程联立来表示。距差与z轴的关系如图2所示,由于拟合的空间圆不在一个平面内,其圆心(中心点)为z向的一系列点,r′1,r′2,r′0分别为被测轮廓的测点到中心点的距离,通过式r′21 -r′20 =Δz2,证明空间坐标系中被测轮廓的测点与中心点的距离r′相对于z轴的距差(Δz)大得多,故Δz对r′的影响可以忽略不计。因此可得出中心点在空间平面上,而空间圆也在空间平面上且平行于xoy平面。

设空间圆的方程为:

(x-a0)2+(y-b0)2+(z-c0)2=R2。

其中:圆心为(a0,b0,c0),半径为R。展开并整理得:

x2+y2+z2-2a0x-2b0y-2c0z+aundefined+bundefined+cundefined=R2 。 (10)

设A=-2a0,B=-2b0,C=-2c0,D=aundefined+bundefined+c20-R2。由式(6)得:

x2+y2+z2+Ax+By+Cz+D=0 。 (11)

将最小二乘法拟合每条线轮廓得到的中心点坐标Oi(ai,bi,ci)代入式(11),方程式明显不等于零,而为:

δi=aundefined+bundefined+cundefined+Aai+Bbi+Cci+D 。 (12)

其中:δi为空间圆上各点的函数偏差值,则:

当undefined为最小值时,空间圆拟合的效果最佳。令:

undefined。 (14)

依据求最小值的方法,对系数A、B、C和D求偏导数后得出关于方程的系数A、B、C和D的线性方程如下:

解得A、B、C、D后,代入计算得:

undefined。

即可得到圆心坐标(a0,b0,c0)和圆半径R。

1.4 求解空间圆与每条线轮廓平面的交点坐标

空间圆与线轮廓平面的方程为:

如式(17)所示,每条线轮廓所在平面与z轴平行且共面,那么该平面的方向数C等于零;最小二乘拟合的空间圆与xoy平面平行且共面,所以每条线轮廓所在平面与空间圆相交,可以得到所有的交点坐标值Li(xi,yi,zi)。而交点坐标值Li(xi,yi,zi)既在平面内又在空间圆内,并且交点与所对应线轮廓的测点在一个象限内。

1.5 轴承沟道形状误差

轴承沟道形状误差模型见图3。

由图3可知,每条线轮廓上的测量点Pij(xij,yij,zij)至所对应的每个交点Li(xi,yi,zi)的距离为:

由轴承沟道形状误差的定义可以知道,每条线轮廓上各测量点至所对应交点的距离中,存在最大值、最小值和极差(最大值与最小值的差值)。因而可以得到M个最大值max(dij)、最小值min(dij)和极差Δdijcha。

以最大值、最小值所在的线轮廓绕z轴转360°形成圆环区域,比较所有的极差值,其中最大值为max(Δdijcha),即得到包容整个轴承沟道的最小二乘形状误差值TLSM为:

TLSM=max{Δdijcha} 。 (19)

2 结论

采用最小二乘算法进行轴承沟道形状误差评定,计算简单,完全满足最小二乘法的评定标准。另外,评定方法不仅可以给出轴承沟道形状误差,还能适用于现代机械的精密设计和精密测量。这种评定方法简单正确,收敛速度快,易于计算机程序实现,具有较好的推广应用价值。

参考文献

[1]王伦.轴承基本知识[M].北京:机械工业出版社,2002.

[2]《数学手册》编写组.数学手册[M].北京:高等教育出版社,2004.

[3]沈世德,徐辛伯.最小二乘圆弧法在图像分析中的应用[J].机械设计与制造,1999,10(5):46-47.

广义逆的计算与最小二乘估计 篇5

在系统理论、信号处理和矩阵方程求解等许多应用领域中,往往需要给出形如:

的矩阵的某种类型的广义逆,其中表示样本数目的k通常都很大,而且随时间增长而增大,因而计算量也迅速增加,这在实时应用中是应该避免的。

本文给出了某些类型的广义逆的算法,并用所得结果,给出了求解由多个矩阵方程所构成的系统的最小二乘估计,从而推广了[1]的工作,并简化了[2]中有关求多个方程公共解的算法,文中使用广义逆的记号均引自[3]。

1广义逆的计算

设k≥1为一整数,为具有相同行数的矩阵且满足

再设的某种广义逆,我们将证明,存在着满足

的矩阵Mk和Gk,使得

为的同种广义逆。

1°,如果(·)-表{1}-逆;

2°,如果(·)-表{1,2}-逆;

1°的证明:注意到Qk-1为对称幂等矩阵且Ak=(AkTAk)(1)AkT=(AkTAk)(1)HkTQk-1[3],我们有

2°的证明:注意到对任意A(1,3)∈A{1,3},总存在着矩阵Z,使得A(1,3)=A-+(I-A-A)Z,因此AA-=AA+,于是由(1),(2)和(3)可知

为对称阵,另外由1°可知:

3°的证明:由2°,我们只需再进一步证明

引理3[5]如果Gk如下确定:

从引理1到引理3知,为获得的广义逆,必须求出

容易看出,以上两式的计算量将随着k的增大而迅速增大,为避免这一困难,注意到

或等价地,

显然,(5)～(8)的计算量与k无关.

综上所述,我们证明了以下两定理。

定理1设Mk和Zk由(1),(2)和(3)所定义且对k=0,1,…,令

Q-1=I,Ak=Qk-1Hk,Qk=Qk-1-AkGk

则

定理2设Mk和Zk由(1),(2)和(3)所定义且对k=0,1,…,令

下面,我们将定理1与定理2应用于

其中和Hk为具有相同列数的矩阵,且注意到如下事实[6]

我们有

定理3设由(9)所定义,且令

则

定理4设和Zk由(9)和(10)所定义且对k=0,1,…,令

则

2最小二乘估计

考虑观测系统

其中Yk∈Rm×p,Hk∈Rm×n为已知,Vk∈Rm×p为零均值随机矩阵,X为待估未知矩阵,对任何k,称X赞k为X的最小二乘估计(LSE),如果X赞k为最优化问题

的最小范数解,而

及加权矩阵wi∈Rm×n皆正定,i=0,1,…,k.

首先考虑无加权情形,即的情形。

为使估计误差的方差最小[6,7],我们有

其中的Moore-Penrose逆.

为得到(14)的递推式,注意到定理4以及(9)、(10)和(13),我们有

定理5系统(11)的X的无加权的极小范数最小二乘估计可如下递推地给出

其中Gk由定理4递推而得。

证明:

如果m=p=1,则定理5便退化为[1]的结果。定理5还有如下简单的推论:

若(11)中的Vk=0,坌k,则定理5给出多个方程所构成的系统

的公共解,亦即我们有

则(16)的公共解可按

递推地求得。

该推论中要求计算Ak+,考虑广义逆的性质和定理3,如果Gk=Qk-1Ak(1),则(17)给出(16)的公共解,从而只需计算A(1)∈A{1}而使计算量大为减少。类似地,如果Gk=Ak(1,4),则(17)给出(16)的具有极小范数的公共解。

现在再考虑一般加权最小二乘估计,因

故可通过无加权问题(18)的右端来处理加权问题。其中

或等价地,

应用定理5于无加权问题(18),便有

定理6令P-1=0,Q-1=I,且对k=0,1,…令

H′kTCk,则系统(11)的加权极小范数二乘估计可如下递推地求得:

为使(11)的加权最小二乘估计是最优估计[8],其中加权矩阵应如何选择?答案如下:

设而其方差设为

则由[4]中的推论6.4.4知,当时,为(11)在给定Yo,…,Yk时的最优最小二乘估计。

关于定理3和定理4的更多应用,可进一步参考[6,9,10]等。

3结论

应用定理1、2、3的结果,将[1]中关于单变量情形下的最小二乘估计推广到了多变量的情形,给出了后者的最优最小二乘估计的算法,同时还简化了[2]中求多个相容矩阵方程公共解的算法。

摘要：给出了某些广义逆的计算方法,将这些结果应用于观测系统的参数估计,给出了系统无加权和加权最小二乘的最新估计结果,推广了[1]的工作,同时简化了[2]中求多个相容矩阵方程公共解的算法。

关键词：算法,最小二乘估计,广义逆,方差,加权

参考文献

[1]Albert A,Sittle R.W.A method for computing leastsquares estimators that keep up with data SIAM[J].Control,1965,3:384-417.

[2]Morris G.L,Odell P.L.Common solutions for a matrixequations with applications[J].Assoc comput.math,1968(15):272-274.

[3]A.B.Israel,T.N.E.Greville.Generalized inverses:Theory and Applications[M].New York:John Wiley,1974.

[4]Campbell S.L,Meyer,C.D.Jr.Generalized Inverses of lin-ear transformations[M].London:Pitman publishing Ltd,1979:104-115.

[5]Cline R.E.Representations for the generalized inverses ofa partitioned matrix[J].Soc.Indust.Appl.Math.1964(12):588-600.

[6]X.Liu(刘轩黄).Recursive computation of generalizedinverses with applications to optimal state estimation[J].Control theory and advanced technology.1995(10):1485-1497.

[7]Rao G.R.Linear statistical inference 2nd its applications(2nd edtion)[M].New York:John Wiley&Sons,1973.

[8]刘轩黄.无初始状态统计知识时多变量系统的最优平滑与预报[J].海南大学学报,1998,16(3):210-214.

[9]X.Liu(刘轩黄).Optimal and minimum energy optimaltracking problems[J].苏州科技学院学报,2005,22(1):1-9;22(2):9-18.

最小二乘建模 篇6

卫星导航定位在现代军事和经济领域扮演着不可或缺的重要角色, 卫星导航系统除了提供导航功能外, 还必须具有在系统不能使用时及时向用户发出告警的能力, 这种能力叫做系统的完好性。完好性和定位精度是卫星导航中两个重要指标。当卫星导航用于航空、武器系统等生命相关场合时, 完好性成为决定性指标。因此, 卫星导航系统的完好性问题已成为卫星导航领域重点解决的问题, 并将存在于卫星导航整个生命周期。

引起接收机伪距观测量偏差的主要因素有较大的卫星钟漂、导航电文的不正确上传、各种欺骗与干扰以及卫星组成部分的故障等, 这里统称为卫星故障, 在GPS系统中这种故障发生的概率大约为10～4/小时。对这种故障, 导航系统本身做出告警反应的时间较长, 如GPS系统的反应时间为15分钟到数小时, 对于生命相关场合显然太慢, 不能满足需求。因此, 卫星故障的快速检测只有在用户端进行, 因而出现了接收机自主完好性检测技术 (Receiver Autonomous Integrity Monitoring, 英文简写RAIM) 。本文针对导航接收机常用的位置解算方法, 即最小二乘法, 设计了基于最小二乘残差法的RAIM检测方案, 并对方案进行了仿真分析。

1 基于最小二乘残差法的RAIM检测方案

导航接收机位置解算常用最小二乘法, 针对最小二乘位置解算方法, 设计了相应的RAIM检测方案, 即基于最小二乘残差法的RAIM检测方案。基于最小二乘残差法的RAIM方案除了对当前可视卫星的几何结构有很强的依赖性, 还对可视卫星数目提出了很高要求。如果接收机要完成故障检测, 最少需要5颗可视卫星;如果接收机要完成故障排除, 则最少需要6颗卫星。

1.1 最小二乘位置解算法[5,6]

对于伪距定位方程:

undefined. (1)

在X0= (x0, y0, z0, B0) ′ (B=cΔt, 看作一个未知数, 以米为单位) 处进行泰勒级数展开[4,5], 略去二次项及高次项得到:

undefined

令undefined

则:undefined. (3)

undefined. (4)

undefined. (5)

将式 (3) (4) (5) 代入式 (2) , 方程可整理为:

ρi-ρ0=αi (x-x0) +βi (y-y0) +γi (z-z0) + (B-B0)

=αiΔx+βiΔy+γiΔz+ΔB. (6)

测量四颗卫星的伪距, 可获得四个类似式 (6) 的方程, 联立这四个方程得到一个四元线性方程组:

令

undefined

ΔX=ΔX ΔY ΔZ ΔBT.

Y=T.

则式 (7) 变为: GΔX=Y. (8)

解此方程得: ΔX=G-1Y. (9)

若测量伪距的卫星数n大于4, 则G为n×4的矩阵, 解此方程得:

ΔX= (GTG) -1GTY. (10)

所以 X=X0+ΔX. (11)

再以新的X作为X0, 继续上述解算, 迭代数次后, 当满足设定的误差时, 迭代结束, 得到所求的解, 实践经验表明, 迭代次数为5左右。

1.2 故障检测与识别

基于最小二乘位置解算原理, 利用最小二乘迭代解算过程的中间结果, 就可以完成RAIM故障检测与识别。引入RAIM故障检测与识别的最小二乘位置解算过程如图1所示, 其中完好性处理部分的详细流程如图2所示。

完好性处理部分首先对观测量Y进行故障检测后, 观测量Y才能进入最小二乘迭代运算过程。完好性处理主要完成可视卫星数目的判断、水平定位误差保护级HPL、检验统计量及其限值的计算, 并对故障进行检测、识别等。具体过程如下:在RAIM检测之前, 有两个重要的限值需要计算, 它们是水平定位误差保护限值HAL和故障检测门限T。水平定位误差保护限值和故障检测门限可通过文献中给出的算法计算得到, 或者可通过文献提出的基于最小风险准则来确定。

在RAIM检测开始时, 首先根据获取的观测信息, 判断可视卫星数目是否达到故障检测与识别要求, 判断可视卫星数目的流程如图3所示。判断可视卫星的主要依据是计算过境卫星的仰角, 如果仰角大于5度, 认为可视。理论上卫星导航系统中最多可视卫星为12颗, 所以判断可视卫星时最多进行12次循环判断。

如果前面判断的可视卫星数目符合要求, 对水平定位误差保护级HPL进行计算, 并与HAL比较, 看是否超限, 以保证RAIM检测的有效性。当HPL小于HAL时, 对检验统计量进行计算。如果检验统计量小于检测门限值T, 则认为无故障, 否则, 进入故障识别与隔离阶段。如果可视卫星数目只有5颗, 那么只能完成故障的检测, 并对检测到的故障向用户发出告警;如果可视卫星数目大于等于6颗, 则可以对检测到的故障进行隔离。最后, 输出完好性处理结果。

2 方案仿真

为了验证基于最小二乘残差法的RAIM方案对故障检测与识别的性能, 利用STK5.0软件产生星座数据, 采用MATLAB软件对该方案进行仿真。由于GPS系统已有较为准确的误差统计数据, 因此仿真采用的卫星星座是GPS卫星星座。利用了在GPS系统时间1 Jan 2010 13:15:00.00至1 Jan 2010 13:17:00.00之间, 在某地区所有的可视卫星, 它们分别是GPS星座系统第1、2、3、6、9、11、15、20颗卫星。仿真条件如下:故障检测误警率为0.005, 由此计算得到故障检测门限值为10.343 m。设观测噪声均值为0, 标准差为5 m。在观测时段内第30秒时刻, 对其中一颗卫星对应的伪距加入一个100 m突变, 模拟此时卫星出现故障, 相应的伪距出现异常。

在上面仿真条件下, 仿真结果如图4所示。横轴为观测时间, 单位为秒。左边纵轴表示检验统计量, 单位为米, 右边纵轴为观测伪距残差, 单位为米。从仿真结果图上可以看出, 在30秒时刻观测伪距发生异常时, 相应的伪距残差出现突变, 从而导致了检验统计量的突变。图中所示的检验统计量在30秒时刻已经超过了故障检测门限值, 因此, 此时的故障完全能够被检测出来。仿真结果表明, 本文所设计的基于最小二乘残差法的RAIM故障检测方案是可行的, 有效的。

3 结束语

完好的系统故障检测性能、快速的故障反馈能力是现代卫星导航定位系统接收机自主完好性检测的发展趋势。本文提出的基于接收机自主完好性检测方案具有告警时问短、运算量小的特点。接收机位置解算如果采用最小二乘法, 就可以很容易地将该方案嵌入其中。我国二代卫星导航定位系统正在建设当中, 建议在系统建设中考虑好系统的完好性检测扩展能力, 在接收机的设计中加入自主完好性检测功能。

摘要：为了使导航接收机的自主完好性检测算法与最小二乘位置解算法结合起来, 提供准确可靠的定位定时服务, 对最小二乘位置解算法和基于最小二乘残差法的接收机自主完好性检测算法进行了研究。针对导航接收机最小二乘位置解算法, 提出了基于最小二乘残差法的卫星导航接收机自主完好性检测方案, 并对方案进行了仿真分析, 仿真结果表明方案是可行的, 有效的。

关键词：导航接收机,自主完好性,最小二乘法

参考文献

[1]Brown R G.Global Positioning System Theory and Appli-cations[M].American:American Institute of Aeronauticsand Astronantics, 1996:143-164.

[2]Xie G, Pullen S, Luo M, et al.Integrity Design and Up-dated Test Results for the Stanford LAAS Integrity MonitorTestbed[C].Albuquerque, NM:Proceedings of ION 57thAnnual Meeting and CIGTF 20th Biennial Guidance TestSymposium, 2001:681-693.

[3]冯庆堂, 沈林成, 常文森.一种监测GPS完好性的新方法[J].导航与控制, 2004, 3 (1) :4-5.

[4]Bradford W, Parkinson, Penina Axelrad.AutonomousGPS Integrity Monitoring Using the Pseudorange Residual[J].Journal of the Institute of Navigation, 1988, 35 (2) :456-460.

[5]王惠南.GPS导航原理与应用[M].北京:科学出版社, 2003.

[6]刘基余.GPS卫星导航定位原理与方法[M].北京:科学出版社, 2003.

[7]黄晓瑞, 田巍, 李波.GPS接收机的自主完善性监测算法研究[J].遥测遥控, 2003, 24 (1) :1-3.

最小二乘建模 篇7

考虑下列二阶抛物问题

$\{\begin{array}{l}u{}_t(x,t)-\nabla \cdot (a(x)\nabla u(x,t))=f(x,t),\\ (x,t)\in \Omega \times J,\\ u(x,t)=0,(x,t)\in \Gamma \times J,\\ u(x,0)=u{}_0(x),x\in \Omega ‚\end{array}(1)$

(1)式中J=[0,T],Ω⊂Rn为有界凸域,边界Γ是Lipschitz连续的。假定a(x)有界,即存在常数α*和α*,使得

0≤α*≤a(x)≤α*,∀x∈Ω (2)

对上述问题人们提出了许多数值模拟格式,如有限元格式[1],混合有限元格式[2],扩展混合有限元格式[3],最小二乘混合有限元格式[4]等。现将最小二乘思想和扩展混合有限元方法相结合,提出并分析了最小二乘扩展混合有限元格式。该格式继承了最小二乘和扩展混合元方法的优点,即格式能同时高精度逼近未知函数,未知函数的梯度以及流体的通量,所形成的有限元方程组是对称正定的,有限元空间无须满足LBB相容性条件。数值算例说明了该格式的有效性。

1 最小二乘扩展混合元格式

定义空间Λ=(L2(Ω))n,V=H${}_0^1$(Ω),W=H(div;Ω)={q∈(L2(Ω))n,∇·q∈L2(Ω)},相应范数为$\parallel \mu \parallel {}_{0,\Omega }=(\mu ,\mu ){}^{\frac{1}{2}},\parallel v\parallel {}_{1,\Omega }=(\parallel v\parallel {}_{0,\Omega }^2+\parallel \nabla v\parallel {}_{0,\Omega }^2){}^{\frac{1}{2}},\parallel q\parallel {}_{{\rm H}(\text{d}\text{i}\text{v};\Omega )}=(\parallel q\parallel {}_{0,\Omega }^2+\parallel \nabla \cdot q\parallel {}_{0,\Omega }^2){}^{\frac{1}{2}}$,其中(·,·)表示空间(L2(Ω))n中的内积。作时间剖分τ=T/N,N为正的常数,在时间tn=nτ,n=1,2,…,N处逼近方程的解。$\text{Ω}={\displaystyle \underset{{\rm K}\in \tau {}_h}{\cup }{\rm K}}$是Ω的一个剖分,h=max{diam(K):K∈Th}。设Vh,Λh和Wh分别是V,Λ和W的有限维子空间,具有通常逼近性质。Vh,Λh和Wh的一个标准的选法是k,r和s次的分片多项式空间。

记λ=-∇u,σ=aλ,并用一阶差分近似ut,则式(1)可变成离散形式I

$\{\begin{array}{l}\overline{\partial }{}_{t}u{}^n=\nabla \cdot \sigma {}^n=f{}^n+R{}_1^n,x\in \Omega ,\\ \lambda {}^n+\nabla u{}^n=0,x\in \Omega ,\\ a\lambda {}^n-\sigma {}^n=0,x\in \Omega ,\\ u{}^n(x)=0,x\in \Omega ,\\ u{}^0(x)=u{}_0(x)x\in \Omega \end{array}(3)$

(3)式中$\overline{\partial }{}_{t}u{}^n=(u{}^n-u{}^{n-1})/\tau ,R{}_1^n=\overline{\partial }{}_{t}u{}^n-u{}_t^n$。

对于(u,λ,σ)∈V×Λ×W和(v,u,q)∈V×Λ×W,定义双线性形式

$\begin{array}{l}a{}_{{\rm N}}(u,\lambda ,\sigma ;\upsilon ,\mu ,q)=(u+\tau \nabla \cdot \sigma ,\upsilon +\tau \nabla \cdot q)+\\ \tau (\lambda +\nabla u,\mu +\nabla \upsilon )+\tau (a{}^{-1}(a\lambda -\sigma ),a\mu -q)(4)\end{array}$

根据离散形式式(3),省略时间截断误差定义最小二乘扩展混合元格式

格式Ⅰ 给定初始值u${}_h^0$∈Vh,对于n=1,2,…,N,求(u,λ,σ)∈Vh×Λh×Wh,使得

$\begin{array}{l}a{}_{{\rm N}}(u{}_n^n,\lambda {}_n^n,\sigma {}_n^n;\upsilon {}_h,\mu {}_h,q{}_h)=(u{}_h^{n-1}+\tau f{}^n,\upsilon {}_h+\\ \tau \nabla \cdot q{}_h),\forall (\upsilon {}_h,\mu {}_h,q{}_h)\in V{}_h\times \Lambda {}_h\times W{}_h。(5)\end{array}$

2 解的存在唯一性

为证明格式Ⅰ解的唯一性,首先说明双线性形式aN(·;·)是强制的。

引理1 对任意的(v,μq)∈V×A×W存在c>0,使得aN(v,μ,q;v,μq)≥c(‖v‖${}_{1,n}^2$+‖μ‖${}_{0,n}^2$+‖q‖H(div n))。 (6)

证明 引入适当选取的常数β>0,式(4)可整理为

$\begin{array}{l}a{}_{{\rm N}}(\upsilon ,\mu ,q;\upsilon ,\mu ,q)=(\tau \nabla \cdot q+(1-\beta )\upsilon ,\\ \tau \nabla \cdot q+(1-\beta )\upsilon )+\tau (q-a(\mu +\beta \upsilon ),\\ a{}^{-1}(q-a(\mu +\beta \upsilon )))+\\ 2\tau (\mu +(1-a\beta )\nabla \upsilon ,\mu +(1-a\beta )\nabla \upsilon )+\\ (2\beta -\beta {}^2)(\upsilon ,\upsilon )。(7)\end{array}$

由于υ(x)=0,x∈Γ,所以由Poincaré-Friedrichs不等式知,存在正常数CF,使得

‖v‖${}_{0,n}^2$≤CF‖∇v‖${}_{0,n}^2$ (8)

结合式(2)有

$\begin{array}{l}a{}_{{\rm N}}(\upsilon ,\mu ,q;\upsilon ,\mu ,q)\ge \beta (2(\tau \alpha {}_{*}+C{}_F)-\\ \beta (\tau \alpha {}^{*2}+\tau \alpha {}^{*}+C{}_F))\parallel \nabla \upsilon \parallel {}_{0,\Omega }^2。(9)\end{array}$

取β=τα*/(τα*2+τα*+CF),则

$a{}_{{\rm N}}(\upsilon ,\mu ,q;\upsilon ,\mu ,q)\ge C\parallel \nabla \upsilon \parallel {}_{0,\Omega }^2\ge C\parallel \upsilon \parallel {}_{1,\Omega }^2。(10)$

从而由式(10)和式(4),知

$\begin{array}{l}\parallel \mu \parallel {}_{0,\Omega }^2\le 2\parallel \mu +\nabla \upsilon \parallel {}_{0,\Omega }^2+2\parallel \nabla \upsilon \parallel {}_{0,\Omega }^2\le Ca{}_{{\rm N}}\\ (\upsilon ,\mu ,q;\upsilon ,\mu ,q),(11)\\ \parallel q\parallel {}_{0,\Omega }^2\le 2\parallel a\mu -q\parallel {}_{0,\Omega }^2+2\parallel a\mu \parallel {}_{0,\Omega }^2\le Ca{}_{{\rm N}}\\ (\upsilon ,\mu ,q;\upsilon ,\mu ,q),(12)\\ \parallel \nabla \cdot q\parallel {}_{0,\Omega }^2\le \frac{2}{\tau }\parallel \upsilon +\tau \nabla \cdot q\parallel {}_{0,\Omega }^2+\frac{2}{\tau }\parallel \upsilon \parallel {}_{0,\Omega }^2\le Ca{}_{{\rm N}}\\ (\upsilon ,\mu ,q;\upsilon ,\mu ,q),(13)\end{array}$

结合式(10)—式(13)便证得式(6)成立。

定理1 设f∈L2(Ω),则格式Ⅰ存在唯一解(μ${}_h^n$,λ${}_h^n$,σ${}_h^n$)∈Vh×Λh×Wh。

证明 考虑空间V×Λ×W,其范数定义为‖υ‖1,Ω+‖μ‖0,Ω+‖Q‖H(div;Ω)。由引理1知双线民生形式αN(·;·)在空间Vh×Λh×Wh中是强制的;而且由Cauchy不等式易证αN(·;·)在Vh×Λh×Wh中是有界的。因此由Lax-Milgram引理知格式Ⅰ解存在唯一。

3 数值算列

为验证所提格式的有效性,在区域[0,1]×[0,1]上考虑方程

$\{\begin{array}{l}u{}_t-(au{}_x){}_x=f,(x,t)\in (0,1)\times (0,1],\\ u(0,t)=u(1,t)=0,t\in [0,1],\\ u(x,0)=\sin \pi x,x\in [0,1]\end{array}(14)$

(14)式中$a=\frac{1}{100},f=\left(1+\frac{\pi {}^2}{100}\right)\text{e}{}^t\sin \pi x$。

容易验证方程的真解为$u=\text{e}{}^t\sin \pi x,\lambda =-\pi \text{e}{}^t\cos \pi x,\sigma =-\frac{1}{100}\pi \text{e}{}^t\cos \pi x$。取Λh为分片常数多项式空间,Vh和Wh为分片线性多项式空间,用MATLAB编程计算近似解uh,λh和σh,绘出真解和最小二乘扩展混合元解的比较图,见图1—图3。

摘要：对二阶抛物问题提出了最小二乘扩展混合元数值模拟格式。数值算例验证了格式的有效性。

关键词：二阶抛物问题,最小二乘,扩展混合有限元方法

参考文献

[1]Thome V.Galerkin finite element methods for parabolic problems.Lecture Notes in Mathematics,Springer,1984

[2]Johnson C V.Thome,Error estimates for some mixed finite element methods for parabolic type problem.RAIRO Numer Anal,1981;15:41—78

[3] Guo Ling,Chen Huanzhen.An expanded characteristic mixed finite element method for a convection-dominated transport problem.Journal of Coputational Mthematics,2005;23:479—490