单位代码
学号SY1401138
分类号
密级
应用数理统计（第一个论文）
山东省旅游发展影响因素多元回归分析
院(系)名称材料科学与工程学院
专业名称材料科学与工程
学生姓名李新杰
任课教师冯伟
2014年12月
摘要
本文主要通过对山东省旅游收入的多因素分析，建立以山东省旅游总收入为因变量，以国内旅游人数、接待入境人数、旅行社总数、旅游饭店数量以及旅游社职工人数等为自变量的多元线性回归模型，并利用SPSS统计软件建立逐步回归模型，找到影响山东省旅游业发展的显著性变量，并对所得的模型给予合理的经济解释。
关键词：逐步回归法山东省旅游发展SPSS 相关性显著性主成分
目录
摘要 (1)
1 引言 (1)
2 数据收集 (2)
3统计数据的初步分析 (3)
3.1 变量间的相关性分析 (3)
3.2一元线性模型的验证 (3)
4 回归分析 (5)
4.1 回归模型的建立 (5)
4.2 回归模型参数的估计 (6)
4.2.1构建模型 (6)
4.2.2剔除变量分析 (7)
4.2.3回归系数分析 (7)
5 回归模型的验证与修正 (8)
5.1 方差分析 (8)
5.2 回归方程的拟合度检验 (9)
5.3 残差检验 (9)
5.3多重共线性检验 (11)
5.4 回归模型的修正 (12)
6 结果 (13)
参考文献 (15)
山东省旅游发展影响因素多元回归分析
学号：SY1401138
姓名：李新杰1 引言
随着社会经济快速发展，生活节奏加快，人们的压力变得越来越大，为减轻压力，既能放松自己，又能拓展自我视野的旅游就成为了人们的首要选择。从我国近5 年的统计数据来看，我国每年的旅游收入正在逐年递增，旅游消费已成为中国人们日常支出中的重要部分。山东省地处黄海之滨和黄河入海口，有着秀丽的自然风光，众多的人文景观，旅游资源十分丰富。全省拥有旅游景区、景点509处，其中泰山和曲阜“三孔”列入世界遗产名录，青岛烟台、威海代表了中国海滨旅游的一大片。全省旅游资源品位高，种类全，分布广，综合条件好，旅游业发展和旅游总收入位于全国前列，为了更好地了解山东省旅游业的发展，对山东省旅游业发展的影响因素建立回归模型分析，找出其核心影响因素。
在应用回归分析去处理实际问题时，必须通过合理经济的方法建立最优回归方程。建立最优回归方程时要注意两个方面：(1)方程中要包含所有的显著作用的自变量，不能遗漏；(2)希望变量个数尽可能少，不含有无意义的变量，而且还应该使这类方程的S达到最小。目前最常用的是逐步回归分析方法，即利用自变量和因变量的一系列同步观测数据，通过对相关矩阵的变换和数理统计的假设检验，逐步把显著性的自变量选入回归方程中，同时也把非显著性的自变量从回归方程中剔除，最终建立一个最优回归方程。
2 数据收集
表2-1 2000-2013年山东省旅游总收入、国内旅游人数、入境旅游人数、旅游社总数、旅游饭
店总数、旅行社职工人数
注：以上数据根据《山东省统计年鉴2000-2013年》整理所得
3统计数据的初步分析
3.1 变量间的相关性分析
为了知道旅游总收入具体和哪些变量有较大的关系，并将这些变量加入到线性模型中，首先要对旅游总收入和5个变量进行相关性分析，得到各个数据之间的相关系数表：
表3-1 各个变量之间的相关系数表
从表3-1可以看出旅游总收入Y和其他变量之间的相关系数，其中旅游总收入（亿元）和X1：国内旅游人数（万人次）的线性正相关程度最高，其次是X2：入境旅游人数（万人次），而旅行社职工人数等相关程度相对较小，所以需要对变量进行一元线性模型验证，以确定是否需要排除掉变量。
3.2一元线性模型的验证
以上我们通过相关性分析确定了各相应变量对旅游总收入Y的影响，为了确定是否需要将所有的变量都加入到线性模型中，下面将通过做出旅游总收入Y
分别和其他5个变量的散点图来进行验证：
（a）(b)
(c)
(d)
（e）
图3-1 因变量和自变量间的散点图：
(a)为旅游总收入Y和国内旅游人数的散点图，(b)为旅游总收入Y和入境旅游人数的散点图，(c)为旅游总收入Y和旅行社总数的散点图，(d)为旅游总收入Y和旅游饭店总数的散点图，
(e)为旅游总收入Y和旅行社职工人数的散点图
从图3-1中的因变量旅游总收入和5个自变量的散点图来看，旅游总收入和5个自变量都有很好的线性关系，这说明通过相关性分析得到的这5个和旅游总收入有关系的自变量都是正确的，而旅游社总数、饭店总数、旅行社职工人数与旅游总收入的相关性差不多，故无需对数据进行删除，因此在接下来进行多元逐步回归分析的时候会将这5个变量都加入到多元线性模型中进行模型建立和分析。
4 回归分析
4.1 回归模型的建立
采用线性回归分析建立的模型为：Y=a+b1X1+b2X2+ … +bnXn；其中Y为因变量的预测值或估计值；X1，X2……Xn为自变量。a和b1、b2…… bn为回归系数。若使以上线性回归分析方法达到最优，就要求自变量满足以下两个条件：（1）在线性回归分析模型中，要包含所有对Y影响显著的自变量，消除对Y 影响不显著的自变量。
（2）模型包含的各自变量之间不存在多重共线性，即各自变量之间不存在线性关系或近似线性关系。
为了解决以上两个问题，最有效的方法是采用逐步回归分析方法。其基本思想是在所考虑的全部因素中，按其对Y作用显著程度的大小，由大到小地逐个引入回归方程。那些对Y作用不显著的变量可能自始至终都未被引人回归方程。另一方面，已被引人回归方程的变量在引入新变量后也可能因为变成对Y作用不显著而从回归方程中剔除。
在回归分析中，对自变量的选择很重要。逐步回归法能使回归式子保留几个最为显著的自变量经过分析，影响山东省旅游收入的主要因素有国内旅游人数、接待入境人数、旅行社总数、旅游饭店数量以及旅游社职工人数，为此设定以下多元线性回归模型：
Y=a + b1X1 + b2X2 + b3X3 + b4X4 + b5X5
其中Y 为山东省旅游总收入（亿元），X1为国内旅游人数（万人次）、X2为入境游客人数（万人次）、X3为旅行社总数(家)、X4为旅游饭店总数（家）、X5为旅行社职工人数（人次）。其中：b i =（1，2，3，4，5）分别表示各变量系数，表示各解释变量对被解释变量Y的影响程度。
4.2 回归模型参数的估计
4.2.1构建模型
通过利用SPSS 软件的线性回归分析，将国内旅游人数（万人次）、入境游客人数（万人次）、旅行社总数(家)、旅游饭店总数（家）、旅行社职工人数（人次）作为自变量，将山东省旅游总收入作为因变量，进行逐步分析法，得到表4-1。
表4-1 输入或者移出到模型中的变量表
从表4-1中可以看到最终模型中存在的自变量是国内旅游人数和入境旅游人数。选择的判据是变量进入回归方程的F的概率不大于0.05，剔除的判据是变量进入回归方程的F的概率不小于0.10。最先进入构建模型1的变量是国内旅游人
数，之后分别进入的是入境旅游人数，构成了模型2。
4.2.2剔除变量分析
表4-2中给出了2种模型中分别不处在其中的相关变量的有关统计量，包括标准化回归系数Beta、回归系数显著性检验的t值、P(Sig)值、偏相关系数、共线性统计容差。
表4-2逐步回归过程中不在模型中的变量
4.2.3回归系数分析
表4-3中给出了各个模型的偏回归系数B、标准差、常数、标准化回归系数、回归系数显著性检验的t值和P值。按照表格数据，最终得到的线性回归的结果为：
Model1：Y= -352.656+0.100X1
Model2：Y= -305.679+0.117X1－1.893X2
两个模型经t检验的P值分别都是0.000，按照=0.10水平，均有显著性意义。
表4-3逐步回归过程的各方程的系数表
5 回归模型的验证与修正
5.1 方差分析
对回归方程的显著性检验就是要看自变量从整体上对随机变量Y是否有明显的影响，主要检验方法有F检验法和t检验法，在这里通过t检验法来进行，在=0.05的水平下得到表5-1所示的方差分析表：
表5-1各个模型的方差分析表
表5-1为方差分析表，从图中可以看出统计量F(1)= 9141.508，F(2)= 13209.207，相伴P值均小于0.000。其明显小于0.05的水平值，拒绝原假设，所以得到的四个模型都是具有显著性水平的，表明Y与X1，X4之间存在线性回归关系。另外，sum of squares 栏中给出了回归平方和，残差平方和和总平方和，
df 为自由度。模型2的回归平方差和残差平方差等均小于第1个模型，说明第2
个模型更好。
5.2 回归方程的拟合度检验
拟合度用于检验回归方程对样本观测值的拟合程度。其值越接近1，表明回
归拟合的效果越好。表5-2展示了逐步回归过程中2个模型的相关系数值。
表5-2 回归过程的拟合度表
从表格中可以看出模型1为只用X1：国内旅游人数（万人次）表示旅游总收
入Y时的相关系数R=0.999，判定系数R2=0.999，调整判定系数为0.999，回归估
计的标准方差为59.60734。
模型2为用X1：国内旅游人数（万人次），X4：旅游饭店总数表示游总收入
Y时的相关系数R=1.000，判定系数R2=1.000，调整判定系数为1.000，回归估
计的标准方差为35.07921。
表中列出了回归方程常用统计量，可以看出，随着模型中自变量个数的增加，
系数的值也在增加，估计的标准误差下降。一般的，修正的值可以比较正确的反
映拟合度。可以看出模型2的回归效果较好，故选择模型2为线性回归方程。
5.3 残差检验
如表5-3为残差统计表，表中显示预测值、残差、标准化预测值、标准化残差的最小值、最大值、均值、标准差及样本容量。准化残差的绝对值最大值<3，说明样本数
据中没有奇异数据。
表5-3残差统计表
图5-1给出了观测累计概率的P-P图：
图5-1 回归标准化残差的标准P-P图图5-2给出了回归残差和标准预测值的散点图：
5.3多重共线性检验
多元线性回归模型的基本假设中要求设计矩阵X中列向量之间不存在密切
的线性关系。若自变量的观测值之间存在线性关系，则称它们之间有着多重共线性。当自变量存在多重共线性时，利用最小二乘法得到的参数估计值很不稳定，回归系数的方差随着共线性强度的增加而加速增长，会造成回归方程高度显著的情况下，所有回归系数都通不过显著性检验，甚至会出现回归系数的正负号无法得到合理解释的的情况。
对于多重共线性的诊断，这里采用的方法是方差扩大因子法（VIF），通过SPSS计算出各变量的VIF值如下表5-4所示：
表5-4 VIF值检验结果
一般认为VIF>10时即说明自变量之间有严重的多重共线性，由表5-4中可以看出，对于模型2，自变量国内旅游人数（X1）和入境旅游人数（X2）的VIF值均为29.486>10，故自变量X1和X2之间存在严重的多重共线性，需要进行修正。
5.4 回归模型的修正
由之前的检验结果表明所建立的多元线性回归模型中自变量之间存在严重的多重共线性，这里采用用主成分分析法对多重共线性进行修正。
主成分分析是一种通过降维方法把多项相关指标化为少数几个不相关综合指标的多元数据处理技术，这里考虑将多维变量（X1，X2）综合成一维指标y1，利用SPSS软件得到主成分分析的结果如下表5-5和表5-6所示：
表5-5 成分矩阵
表5-6 解释的总方差
由表5-5可以得到第一主成分y1=0.991X1+0.991X2，由表5-6可知第一主成分y1的累计方差贡献率为99.145%，即第一主成分y1从原始变量X1，X2中提取了99.145%的信息量。因此，在累计方差贡献率如此高的情况下可以认为用主成分分析法对模型进行修正是合理的。
建立主成分y1和因变量旅游总收入（y）的线性回归模型，分析的结果如下表5-7所示：
表5-7 主成分分析模型
因此，主成分分析方法最终所确定的线性回归模型如下式所示：
y= -358.611+0.100y1
其中y1=0.991X1+0.991X4，即为第一主成分，累计方差贡献率为99.145%。
6 结果
（1）本文以山东省旅游业发展为例，利用统计知识里面的多元线性回归分析的方法对旅游总收入、国内旅游人数、入境旅游人数、旅行社总数、旅游饭店总数、旅行社职工人数之间的关系进行了分析，得到了山东省旅游发展的多元线性回归方程为：
Y= -305.679+0.117X1－1.893X2
式中：Y为旅游总收入，X1为国内旅游人数，X2为入境旅游人数。
（2）对拟合出的多元线性回归模型的优劣性进行了检验，包括方差诊断、残差检验、多重共线性检验。检验结果表明，建立的模型中自变量X1和X2之间存在严重的多重共线性。
（3）针对自变量之间出现的严重的多重共线性现象，采用主成分分析法对模型进行了修正，最终得到的线性回归模型为：
y=0.100y1
其中y1=0.991X1+0.991X4，即为第一主成分，累计方差贡献率为99.145%。可以认为建立了合适的线性回归模型。
（4）建立的回归模型表示山东旅游总收入和国内旅游人数以及入境旅游人数呈显著的线性相关关系，而把旅游社总数、旅游饭店总数、旅游社职工人数排除在
外，这表明国内旅游人数和入境旅游人数对旅游业发展有较大贡献的经济指标。
参考文献
[1] 张文彤，闫洁. SPSS统计分析基础教程[M]. 北京：高等教育出版社，2004
[2] 《山东省统计年鉴》，2014.
[3] 韩於羹. 应用数理统计[M]. 北京：北京航空航天大学出版社，2002.
[4] 孙海燕周梦李卫国冯伟. 应用数理统计[M]. 北京：北京航空航天大学出
版社，2010.
《数值分析B》大作业一 SY1103120 朱舜杰 一．算法设计方案： 1.矩阵A的存储与检索 将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] . 由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是： A的带内元素a ij=C中的元素c i-j+2,j 2.求解λ1，λ501，λs ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs； 如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。 ②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求 出对应的按摸最大的特征值λ，max， 如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。 3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik （k=1,2，…，39）。 使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。 4.求解A的（谱范数）条件数cond（A）2和行列式d etA。 ①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和 最小特征值。
②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。 二．源程序 #include 《数值分析》计算实习题目 第一题： 1. 算法设计方案 （1）1λ，501λ和s λ的值。 1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。 2)使用反幂法求λs ，其中需要解线性方程组。因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。 （2）与140k λλμλ-5011=+k 最接近的特征值λik 。 通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。 （3）2cond(A)和det A 。 1）1=n λλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。 2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。 由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。 2.全部源程序 #include 数值分析第三次大作业 一、算法的设计方案： （一）、总体方案设计： x y当作已知量代入题目给定的非线性方程组，求（1）解非线性方程组。将给定的(,) i i 得与(,)i i x y 相对应的数组t[i][j],u[i][j]。 （2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=(,)i i f x y 。 （3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。 （4）观察和(,)i i p x y 的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,)i i f x y ，再与对应的(,)i i p x y 比较即可，这里求解 (,)i i p x y 可以直接使用（3）中的C[r][s]和k 。 （二）具体算法设计： （1）解非线性方程组 牛顿法解方程组()0F x =的解* x ，可采用如下算法： 1）在* x 附近选取(0) x D ∈，给定精度水平0ε>和最大迭代次数M 。 2）对于0,1, k M =执行 ① 计算() ()k F x 和()()k F x '。 ② 求解关于() k x ?的线性方程组 () ()()()()k k k F x x F x '?=- ③ 若() () k k x x ε∞∞ ?≤，则取*()k x x ≈，并停止计算；否则转④。 ④ 计算(1) ()()k k k x x x +=+?。 ⑤ 若k M <，则继续，否则，输出M 次迭代不成功的信息，并停止计算。 （2）分片双二次插值 给定已知数表以及需要插值的节点，进行分片二次插值的算法： 设已知数表中的点为： 00(0,1,,) (0,1,,)i j x x ih i n y y j j m τ=+=???=+=?? ，需要插值的节点为(,)x y 。 1) 根据(,)x y 选择插值节点(,)i j x y ： 若12h x x ≤+ 或12 n h x x ->-，插值节点对应取1i =或1i n =-， 目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)() 1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10) 算法： 以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法： ()[]()() []()[]()111111I 00000 i n n n B A I gause i n Q A I u Bu u λλ-?-?-=-?-?? ?-=????→=??????→= ?? ? 选主元的消元 检查知无重特征值 由于=0i A I λ- ，因此在经过选主元的高斯消元以后，i A I λ- 即B 的最后一行必然为零，左上方变 为n-1阶单位矩阵[]()()11I n n -?-，右上方变为n-1阶向量[]()11n Q ?-，然后令n u 1=-，则 ()1,2,,1j j u Q j n ==???-。 这样即求出所有A所有实特征值对应的一个特征向量。 #include 数值分析第二次大作业 史立峰 SY1505327 一、 方案 （1）利用循环结构将sin(0.50.2)() 1.5cos( 1.2)() {i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值，得到需要变换的 矩阵A ； （2）然后，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵A (n-1)。 对A 拟上三角化，得到拟上三角矩阵A (n-1)，具体算法如下： 记A(1)=A ，并记A(r)的第r 列至第n 列的元素为()n r r j n i a r ij ,,1,;,,2,1) ( +==。 对于2,,2,1-=n r 执行 1. 若 ()n r r i a r ir ,,3,2) ( ++=全为零，则令A(r+1) =A(r),转5；否则转2。 2. 计算 () ∑+== n r i r ir r a d 1 2 )( ()( )r r r r r r r r r r d c a d a c ==-=++则取,0sgn ) (,1)(,1若 )(,12r r r r r r a c c h +-= 3. 令 () n T r nr r r r r r r r r R a a c a u ∈-=++) ()(,2)(,1,,,,0,,0 。 4. 计算 r r T r r h u A p /)(= r r r r h u A q /)(= r r T r r h u p t /= r r r r u t q -=ω T r r T r r r r p u u A A --=+ω)()1( 5. 继续。 （3）使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值，也是A 的全部特征值，具体算法如下： 1. 给定精度水平0>ε和迭代最大次数L 。 2. 记n n ij n a A A ?-==][) 1()1()1(，令n m k ==,1。 《数值分析A》 一、算法设计方案 整个程序主要分为四个函数，主函数，拟上三角化函数，QR分解函数以及使用双步位移求解矩阵特征值、特征向量的函数。因为在最后一个函数中也存在QR分解，所以我没有采用参考书上把矩阵M进行的QR分解与矩阵Ak的迭代合并的方法，而是在该函数中调用了QR分解函数，这样增强了代码的复用性，减少了程序长度；但由于时间关系，对阵中方法的运算速度没有进行深入研究。 1.为了减少QR分解法应用时的迭代次数，首先对给定矩阵进行拟上三角化处理。 2.对经过拟上三角化处理的矩阵进行QR分解。 3.注意到计算特征值与特征向量的过程首先要应用前面两个函数，于是在拟上三角化矩阵的基础上对QR分解函数进行了调用。计算过程中，没有采用goto语句，而是根据流程图采用其他循环方式完成了设计，通过对迭代过程的合并，简化了程序的循环次数，最后在计算特征向量的时候采用了列主元高斯消去法。 二、源程序代码 #include 北京航空航天大学 数值分析大作业八 学院名称自动化 专业方向控制工程 学号 学生姓名许阳 教师孙玉泉 日期2014 年11月26 日 一．题目 关于x , y , t , u , v , w 的方程组(A.3) ???? ?? ?=-+++=-+++=-+++=-+++79 .0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0y w v u t x w v u t y w v u t x w v u t (A.3) 以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z =f (x , y )。 表A-1 二维数表 t z u 0 0.4 0.8 1.2 1.6 2 0 -0.5 -0.34 0.14 0.94 2.06 3.5 0.2 -0.42 -0.5 -0.26 0.3 1.18 2.38 0.4 -0.18 -0.5 -0.5 -0.18 0.46 1.42 0.6 0.22 -0.34 -0.58 -0.5 -0.1 0.62 0.8 0.78 -0.02 -0.5 -0.66 -0.5 -0.02 1.0 1.5 0.46 -0.26 -0.66 -0.74 -0.5 1. 试用数值方法求出f (x , y ) 在区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的近似表达式 ∑∑===k i k j s r rs y x c y x p 00 ),( 要求p (x , y )以最小的k 值达到以下的精度 ∑∑==-≤-=10020 7210)],(),([i j i i i i y x p y x f σ 其中j y i x i i 05.05.0,08.0+==。 2. 计算),(),,(* ***j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值，以观察p (x , y ) 逼 近f (x , y )的效果，其中j y i x j i 2.05.0,1.0**+==。 《数值分析A》计算实习题目第一题 一．算法设计方案： 1.矩阵A的存储与检索 将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] . 由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是： A的带内元素a ij=C中的元素c i-j+2,j 2.求解λ1，λ501，λs ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs； 如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。 ②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ，max， 如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。 3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik （k=1,2，…，39）。 使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。 4.求解A的（谱范数）条件数cond（A）2和行列式d etA。 ①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。 ②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有 对角线上元素的乘积。 二．源程序(VS2010环境下，C++语言) #include 一、算法设计方案 1.使用牛顿迭代法，对原题中给出的i x i 08.0=，j y j 05.05.0+=， (010 ,020i j ≤≤≤≤)的11*21组j i y x ,分别求出原题中方程组的一组解，于是得到一组和i i y x ,对应的j i t u ,。 2.对于已求出的j i t u ,，使用分片二次代数插值法对原题中关于u t z ,,的数表进行插值得到 ij z 。于是产生了z=f(x,y)的11*21个数值解。 3.从k=1开始逐渐增大k 的值，并使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，得到每次的σ,k 。当7 10-<σ时结束计算，输出拟合结果。 4.计算)5,,2,1,8,,2,1)(,(),,(* ***???=???=j i y x p y x f j i j i 的值并输出结果，以观察),(y x p 逼近),(y x f 的效果。其中j y i x j i 2.05.0,1.0* *+==。 二、算法实现方案 1、求(,)f x y ： （1）Newton 法解非线性方程组 0.5cos 2.670.5sin 1.07(1)0.5cos 3.740.5sin 0.79 t u v w x t u v w y t u v w x t u v w y +++-=??+++-=? ? +++-=??+++-=?， 其中，t, u, v ,w 为待求的未知量，x, y 为代入的已知量。 设(,,,)T t u v w ξ=，给定精度水平12110ε-=和最大迭代次数M ，则解该线性方程组的迭代格式为： *(0)(0)(0)(0)(0)(k+1) ()()1()(,,,)()()0,1,T k k k t u v w F F k ξξξ ξξξ-?=?'=-??= ? 在附近选取初值， 迭代终止条件为()(1) () 1/k k k ξξ ξε-∞ ∞ -≤，若k M >时仍未达到迭代精度，则迭代计算失 败。 其中，雅可比矩阵 0.5*cos(t) + u + v + w - x - 2.67t + 0.5*sin(u) + v + w - y - 1.07()0.5*t + u + cos(v) + w - x - 3.74t + 0.5*u + v + sin(w) - y - 0.79F ξ???? ? ?=?????? ， 北京航空航天大学2020届研究生 《数值分析》实验作业 第九题 院系：xx学院 学号： 姓名： 2020年11月 Q9:方程组A.4 一、 算法设计方案 （一）总体思路 1.题目要求∑∑=== k i k j s r rs y x c y x p 00 ),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。 ),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。 2.),(* * j i y x f 与1使用相同方法求得，),(* * j i y x p 由计算得出的p(x,y)直接带入),(* * j i y x 求得。 1. ),(i i y x 与),(i i y x f 的数表的获得 对区域D ={ （x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ，yj=1+0.008j ，得到),(i i y x 共31×21组。将每组带入A4方程组，即可获得五个二元函数组，通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。再将 ),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。 2.乘积型最小二乘曲面拟合 2.1使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下： ????? ??=k i i k x x x x B 0000 ， ????? ??=k j j k y y y y G 0000 数表矩阵如下： ???? ? ? ?=),(),(),(),(0000j i i j y x f y x f y x f y x f U 记C=[rs c ]，则系数rs c 的表达式矩阵为： 11-)(-=G G UG B B B C T T T ）（ 通过求解如下线性方程，即可得到系数矩阵C 。 UG B G G C B B T T T =)()( 2.2计算),(),,(* ***j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值 ),(**j i y x f 的计算与),(j i y x f 相同。将),(**j i y x 代入原方程组，求解响应) ,(* *ij ij u t 进行分片双二次插值求得),(**j i y x f 。),(* *j i y x p 的计算则可以直接将),(**j i y x 代入所求p(x,y)。 二、 源程序 ********* 第三次数值分析大作业Q9************ integer::i, j, K1, L1, n, m dimension X(31), Y(21), T(6), U(6), Z(6, 6), UX(11, 21), TY(11, 21), FXY(11, 21), C(6, 6) dimension z1(31, 21), z2(31, 21), z3(31, 21), z4(31, 21), z5(31, 21) dimension X1(8), Y1(5), FXY1(8, 5), PXY1(8, 5), UX1(8, 5), TY1(8, 5) 《数值分析》计算实习作业 （第二题） 算法设计方案： 1、对矩阵A 赋值，取计算精度ε=1×10-12； 2、对矩阵A 进行拟上三角化，得到A (n-1)，并输出A (n-1)； 对矩阵A 的拟上三角化，通过直接调用子函数inftrianglize(A)来实现；拟上三角化得到的矩阵A (n-1)输出至文件solution.txt 中。 3、对A (n-1)进行QR 分解并输出Q 、R 及RQ 矩阵； QR 分解通过直接调用子函数QRdescom(A,Q,R, n)实现。 4、运用QR 方法求所有的特征值，并输出； (1)初始时令m=n ，在m>2的条件下执行； (2)判断如果|A mm-1|<ε，则得到一个特征值，m=m-1，转(4)；否则转(3)； (3)判断如果|A m-1m-2|<ε，则得到两个特征值，m=m-2，转(4)； (4)判断如果m ≤2，转(6)；否则转(5)； (5)执行相似迭代，转(2)； k k T k k k k k k k k k k Q A Q A R Q M I D A D tr A M ==+-=+1)2)det(( (6)求出最后的一个或两个特征值； (7)输出全部的特征值至文件solution.txt 中。 5、输出QR 分解法迭代结束之后的A (n-1)至文件solution.txt 中； 6、通过反幂法求出所有实特征值的特征向量并输出。 首先令B=(A-λi I)，其中λi 是实特征值；反幂法通过调用子函数Bpowmethod(B,x1)实现，最终λi 对应的特征向量就是x1；最后将所有的实特征值的特征向量输出。 数值分析第二题 梁进明 SY0906529 算法设计方案。 一．矩阵的QR 分解。把矩阵A 分解为一个正交矩阵Q 与一个上三角矩阵R 的乘积，称为矩阵A 的 正三角分解，简称QR 分解。QR 分解的算法如下： 记1A A =，并记[]r ij n n Ar a ?=，令1Q I =（n 阶单位矩阵） 对于r=1,2,…,n-1执行 (1) 若(1,2,...,)r ir a i r r n =++全为零，则令1r r Q Q +=，1r r A A +=转(5);否则转(2) (2) 计算 2r d = ()sgn()r r rr r c a d =-(若() 0r rr a =,则取r r c d =) 2() r r r r rr h c c a =- (3) 令()()()1,(0,...,0,,,...,)r r r T n r rr r r r nr u a c a a R +=-∈ (4) 计算 11//r r r T r r r r r T r r r r T r r r r Q u Q Q u h p A u h A A u p ωω++==-==- (5) 继续 当此算法执行完后就得到正交矩阵n Q Q =和上三角矩阵n R A =且有A QR =。 二．矩阵的 拟上三角化。对实矩阵A 的拟上三角化具体算法如下： 记(1) A A =，并记()r A 的第r 列到第n 列的元素为(1,2,...,;,1,...,)r ij a i n j r r n ==+。 对于1,2,...,2r n =-执行 （1） 若() (2,3,...,)r ir a i r r n =++全为零，则令(1) ()r r A A +=,转(5);否则转(2)。 （2） 计算 《数值分析B》 第三次数值分析大作业 院系：04 能源与动力工程学院 姓名：王开逍 学号：SY1104207 一．算法设计方案： 1、使用牛顿迭代法，对原题中给出的X(I)=0.08*I，Y(J)=0.5+0.05*J的11*21组X(I),Y(J)分别求得原题非线性方程组的一组解，于是得到一对和X(I),Y(J)对应的T(I),U(J). 2、对于已经求出来的U(I)，T(J)，使用分片二次代数插值法对原题中关于Z,T,U的数表进行插值，得到Z(I,J).于是产生了Z=F(X,Y)的11*21个数值解. 3、从K=1开始逐渐增大K的值，并使用最小二乘曲面拟合法对Z=F(X,Y)进行拟合，得到每次的K和精度TAO，当TAO 北航数值分析大作业第二题(f o r t r a n) -CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN “数值分析“计算实习大作业第二题 ——SY1415215 孔维鹏一、计算说明 本程序采用带双步位移的QR方法求解矩阵A的所有特征值，然后采用反幂法求解矩阵A的实特征值对应的特征向量。 在采用带双步位移的QR方法求解特征值时，对教材上所提供的具体算法作稍微的改动，以简化程序，具体算法如下所示： 1、计算出A拟上三角化后的矩阵，给定精度水平和最大迭代次数L； 2、记，令k=1，m=n； 3、如果，则可直接计算出最后1或2个特征值，转8，否则转4； 4、如果，则可得一个特征值，置m=m-1；转3，否则转5； 5、如果，则可得两个特征值，置m=m-2；转3，否则转6； 6、记，计算 7、k=k+1,转3 8、A的全部特征值已经求出，停止计算。 二、计算源程序（FORTRAN） PROGRAM SY1415215_2 PARAMETER (N=10) DIMENSION A(N,N),A1(N,N),A2(N,N),C(2,N),Q(N,N),R(N,N),CR(N),CM(N)!C为存储特征值的数组,1为实部，为虚部 REAL(8) A,A1,A2,C,Q,R,CM E=1E-12 !精度水平 L=1000 !迭代最大次数 OPEN(1,FILE='数值分析大作业第二题计算结果.TXT') DO I=1,N DO J=1,N IF(I==J) THEN A(I,J)=*COS(I+*J) ELSE A(I,J)=SIN*I+*J) ENDIF ENDDO ENDDO A1=A WRITE(*,"('矩阵A为：')") WRITE(1,"('矩阵A为：')") DO I=1,N DO J=1,N WRITE(*,"(2X,,2X,\)") A(I,J) WRITE(1,"(2X,,2X,\)") A(I,J) 数值分析第三次大作业 一、算法的设计方案 1、求解非线性方程组 将题目中给出的(,)i i x y 当作已知量代入题目给定的非线性方程组，求出与 (,)i i x y 相对应的数组te[i][j]，ue[i][j]，此处采用的是牛顿法解非线性方程组，其 算法如书上91页所示。 2、分片二次代数插值 对所求出的数组te[i][j]，ue[i][j]，通过分片二次代数插值运算，得到与数组te[11][21]，ue[11][21]对应的数组ze[11][21]，从而得到二元函数z=(,)i i f x y ，此处采用如书上101页例2中所示的分片二次代数插值。 3、曲面插值 利用x[11]，y[21]，ze[11][21]建立二维函数表，进行曲面插值计算，逐步提高k 值，计算其精度，看其是否满足要求，以此来确定循环结束的时刻，并得到曲面拟合的系数矩阵C[r][s]，此处的算法如书142页所示，只需将所需矩阵给出，然后按公式进行计算即可。 4、比较 观察和),(j i y x p 逼近(,)i i f x y 的效果。观察逼近效果只需要利用新给的点列 (,)i i x y 重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,) i i f x y ， 再与对应的(,)i i p x y 比较即可，这里求解(,)i i p x y 可以直接使用（3）中的C[r][s]和k 。 5、几点说明 分片二次插值的结果x[i]，y[j]，ze[i][j]输出到一个文件shubiao.txt 中，方便结果的复制与粘贴。 曲面插值的结果输出到一个文件xishu.txt 中，包括循环中每一次的k 值以及误差平方和sigma 的值，还有最后满足误差要求时曲面插值的系数C[r][s]。 观察逼近效果的结果输出到一个文件shubiao1.txt 中，方便结果的复制与粘贴。 一、算法设计方案： 按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。 总体思路： 1）初始化矩阵 首先需要将需要求解的矩阵输入程序。为了防止矩阵在后面的计算中被破坏保存A[][]。 2）对给定的矩阵进行拟上三角化 为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。由于矩阵的QR 分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。这里用函数 void QuasiTriangularization()来实现，函数形参为double型N维方阵double a[][N]。 3）对拟上三角化后的矩阵进行QR分解 对拟上三角化的矩阵进行QR分解会大大减小计算量。用子程序void QR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。 4）对拟上三角化后的矩阵进行带双步位移的QR分解 为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。这里用函数intTwoStepDisplacement_QR()来实现。这是用来存储计算得到的特征值的二维数组。考虑到特征值可能为复数，因此将所有特征值均当成复数处理。此函数中，QR分解部分用子函数void QR_decompositionMk()实现。这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。 5）计算特征向量 得到特征值后，计算实特征值相应的特征向量。这里判断所得特征值的虚数部分是否为零。求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。因此先初始化矩阵Array，计算(A-λI)，再求解方程组。 方程组的求解采用列主元的高斯消去法，由函数intGauss_column(double a[][N],double b[],double X[])实现。由于此给定矩阵的特殊性，其没有重根，所有对应于每一特征值只有一个特征向量，因此可以用高斯消去法求解此奇异的线性方程组。首先用高斯消去将矩阵(A-λI)化为上三角阵，其最后一行必定全为零。因此在反代的过程中，只要令x[]的最后一个元素为“1”，即可得到方程组的一个基础解系，从而也就是一个特征向量。 6）输出有关结果 根据题目要求，需要输出拟上三角化后的矩阵、QR分解结束后的矩阵、矩阵全部特征值及对应实特征值的特征向量。由于输出结果要求都要保留12位有效数字，所以将结果输出到文件result.txt中更有利于数据的打印。程序中构造了两个输出函数专门来解决不同输出结果的问题，void print_lambda(complex lambda[]);void print_matrix(double mat[][N])。 《数值分析A》计算实习题目二 姓名 学号 数值分析第二次大作业 一、算法设计方案 首先构造矩阵A，利用Householder矩阵对矩阵A作相似变换，把A化为拟上三角矩阵A（n-1），算法如课本P61。 使用QR分解法对矩阵A（n-1）进行QR分解，算法如课本P59, 进而求出所得矩阵的Q、R、RQ矩阵。 然后对A（n-1）进行带双步位移的QR分解求矩阵的全部特征值，采用以下几步进行： 第一步：判断是否a m,m-1(k)<=ε ,若不是,则进入第四步。若是，则a m,m-1(k)为特征值，m=m-1，若m=1，则进入第二步，若m=2进入第三步，否则转第四步。 第二步：m=1，则a11(k)为特征值，转向结束步。 第三步：m=2,则可以求出A的两个特征值s1和s2，转向结束步。 第四步：判断是否a m-1,m-2(k)<=ε，若不是，进入第五步。若是，则得到A的两个特征值s1和s2，令m=m-2，若m=1，进入第二步，若m=2进入第三步，否则进入第一步。 第五步：判断是否达到循环上限，若达到，则结束，否则进入第六步。 第六步：对A进行双步位移QR分解，这里的算法如课本P64，分解后转向计数步。 计数步：对循环次数进行计数，并转向第一步。 结束步：显示所求得的特征值。 最后对实特征值利用列主元高斯消元法求解其对应的特征向量，算法如课本p17. 二、源程序代码 #include 北航数值分析大作业第一题幂法与反幂法 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII 《数值分析》计算实习题目 第一题： 1. 算法设计方案 （1）1λ，501λ和s λ的值。 1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。 2)使用反幂法求λs ，其中需要解线性方程组。因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。 （2）与140 k λλμλ-5011=+k 最接近的特征值λik 。 通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。 （3）2cond(A)和det A 。 1）1=n λλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。 2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。 由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。 2.全部源程序 #include 一、算法设计方案 1、解非线性方程组 将各拟合节点（x i ，y j ）分别带入非线性方程组，求出与 (,)i i x y 相对应的数组 te[i][j]，ue[i][j]，求解非线性方程组选择 Newton 迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle 分解法。 2、二元二次分偏插值 对数表z(t,u)进行分片二次代数插值，求得对应（t ij ，u ij ）处的值，即为),(j i y x f 的值。根据给定的数表，可将整个插值区域分成 16 个小 的区域，故先判断t ij , u ij 所在，的区域，再作此区域的插值，计算 z ij ，相应的Lagrange 形式的插值多项式为： °11 2211 (,)()()(,)m n k r k r k m r n p t u l t l u f t u ++=-=-= ∑∑ 其中 1 1()m w k w m k w w k t t l t t t +=-≠-= -∏ (k=m-1, m, m+1) °1 1()n w r w n r w w r y y l u y y +=-≠-= -∏ (r=n-1, n, n+1) 3、曲面拟合 从k=1开始逐渐增大k 的值，使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，当710-<σ时结束计算。拟合基函数φr (x)ψs (y)选择为φr (x)=x r ，ψs (y)=y s 。拟合系数矩阵c 通过 连续两次解线性方程组求得。[]rs c * =C ，11()()T T T --=C B B B UG G G 其中 01 110101 1 [()]1 k k r i k x x x x x x x ?????? ?==?? ??????B L L M M M M L ，0 01 1101011 [()]1 k k s j k y y y y G y y y ψ????? ? ==???????? L L M M M M L [(,)]i j f x y =U 4、观察比较 计算)5,,2,1,8,,2,1)(,(),,(****???=???=j i y x p y x f j i j i 的值并输出结果，以观察),(y x p 逼近),(y x f 的效果。其中j y i x j i 2.05.0,1.0**+==。 二、全部源程序 // hean.cpp : 定义控制台应用程序的入口点。 // #include "stdafx.h" #include 北航数值分析大作业第一题幂法与反幂法
北航数值分析报告第三次大作业
北航数值分析大作业第二题精解
北航数值分析大作业第二题
北航数值分析第二次大作业--QR分解
北航数值分析报告大作业第八题
北航数值分析课程第一次大作业讲解
北航数值分析大作业3
BUAA数值分析大作业三
北航数值分析大作业第二次
北航数值分析大作业 第二题 QR分解
北航数值分析大作业3
北航数值分析大作业第二题(fortran)
北航数值分析第三次大作业
北航 数值分析第二次大作业(带双步位移的QR方法)
北航数值分析第二次大作业--QR分解
北航数值分析大作业第一题幂法与反幂法
北航数值分析A大作业3
相关主题
文本预览