一、算法设计方案：
按题目要求，本程序运用带双步位移的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])。
二、程序源代码：
#include "stdafx.h"
#include "stdlib.h"
#include "math.h"
#define L 100 //定义最大迭代次数
#define N 10 //定义矩阵大小
#define err 1e-12 //定义误差限
//定义一个结构体来表示复数
typedefstruct complex{
double rpart;
double ipart;
};
FILE * pReadFile;
//主函数
int _tmain(intargc, _TCHAR* argv[])
{
inti,j,t;
double y[N],X[N],a[N][N],A[N][N],B[N][N],Q[N][N],R[N][N],RQ[N][N];
struct complex lambda[N];
//声明要调用的函数
void QuasiTriangularization(double a[][N]);
void QR_decomposition(double A[][N],double Q[][N],double R[][N]);
void QR_decompositionMk(double mk[][N],int m, double a[][N]);
void print_lambda(complex lambda[]);
void print_matrix(double mat[][N]);
void multi_matrix(double a[][N],double b[][N],double c[][N]);
intTwoStepDisplacement_QR(double a[][N],complex lambda[]);
intGauss_column(double a[][N],double b[],double X[]);
//矩阵初始化
for (i = 0; i < N; i++){
for (j = 0; j < N; j++){
if (i == j){
a[i][j] = 1.5 * cos((i+1) + 1.2 * (j+1));
A[i][j] = a[i][j];
}
else{
a[i][j] = sin(0.5 * (i+1) + 0.2 * (j+1));
A[i][j] = a[i][j];
}
}
}
for (i = 0; i < N; i++){
y[i] = 0;
}
//对矩阵a[][]拟上三角化
QuasiTriangularization(a);
//打开文件result.txt
pReadFile = fopen( "result.txt", "w" );
//打印结果到文件result.txt中
fprintf(pReadFile,"拟上三角化之后的矩阵A[%d][%d]:\n",N,N);
//printf("拟上三角化之后的矩阵A[%d][%d]:\n",N,N);
print_matrix(a);
//对拟上三角化后的矩阵a[][]，进行QR分解
QR_decomposition(a,Q,R);
fprintf(pReadFile,"Q[%d][%d]:\n",N,N);
//printf("Q[%d][%d]:\n",N,N);
print_matrix(Q);
fprintf(pReadFile,"R[%d][%d]:\n",N,N);
//printf("R[%d][%d]:\n",N,N);
print_matrix(R);
multi_matrix(R,Q,RQ);
fprintf(pReadFile,"RQ[%d][%d]:\n",N,N);
//printf("RQ[%d][%d]:\n",N,N);
print_matrix(RQ);
//双步位移QR分解求全部特征值
TwoStepDisplacement_QR(a,lambda);
//利用校正的列主元素高斯消元法求每个实特征值对应的特征向量
for (t = 0; t < N; t++){
if (lambda[t].ipart == 0){
for (i = 0; i < N; i++){
for (j = 0; j < N; j++){
if (i == j)
B[i][j] = A[i][j] - lambda[1].rpart;
else
B[i][j] = A[i][j];
}
}
Gauss_column(B,y,X);
fprintf(pReadFile,"\n矩阵与特征值λ=%.12e对应的特征向量为X[%d]：\n",lambda[t].rpart,N);
//printf("\n矩阵与特征值λ=%.12e对应的特征向量为X[%d]：\n",lambda[t].rpart,N);
for (i = 0; i < N; i++){
fprintf(pReadFile,"%.12e ",i,X[i]);
//printf("X[%d]: %.12e ",i,X[i]);
}
}
}
fclose( pReadFile );
scanf("%d",&i);
return 0;
}//主函数
/*************************************
矩阵的拟上三角化
输入n阶方阵a[][]，将a[][]拟上三角化
无返回值
**************************************/
void QuasiTriangularization(double a[][N]){ intr,i,j;
double tr,hr,cr,dr,sum;
double ur[N],pr[N],qr[N],wr[N];
for (r = 0; r < N-2; r++){
//判断a[i][r](i=r+2,r+3,...,n-2)是否全为零
sum = 0;//变量sum使用前清零
for (i = r+2; i < N; i++){
sum = sum || a[i][r];
}
//如果不是全部都是零，则计算
if (sum){
//计算dr
sum = 0;
for (i = r+1; i < N; i++){
sum += a[i][r] * a[i][r];
}
dr = sqrt(sum);
//计算cr
if (a[r+1][r] > 0)
cr = -dr;
else
cr = dr;
//计算hr
hr = cr * cr - cr * a[r+1][r];
//取向量ur[]
for (i = 0; i < N; i++){
if (i < r+1)
ur[i] = 0;
else
if (i == r+1)
ur[i] = a[i][r] - cr;
else
ur[i] = a[i][r];
}
//计算向量qr[]
for (i = 0; i < N; i++){
sum = 0;
for (j = r+1; j < N; j++)
sum += a[i][j] * ur[j];
qr[i] = sum / hr;
}
//计算向量pr[]
for (i = 0; i < N; i++){
sum = 0;
for (j = r+1; j < N; j++)
sum += a[j][i] * ur[j];
pr[i] = sum / hr;
}
//计算tr
sum = 0;
for (i = r+1; i < N; i++)
sum += pr[i] * ur[i];
tr = sum / hr;
//计算wr[]
for (i = 0; i < N; i++){
if (i < r+1)
wr[i] = qr[i];
else
wr[i] = qr[i] - tr * ur[i];
}
//计算新产生的矩阵a[][]
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
a[i][j] = a[i][j] - wr[i] * ur[j] - ur[i] * pr[j];
}
}
}
/*************************************
矩阵的QR分解
将A[][]分解为Q[][]*R[][]
无返回值
**************************************/
void QR_decomposition(double A[][N],double Q[][N],double R[][N]){ intr,i,j;
double tr,hr,cr,dr,sum;
double ur[N],pr[N],wr[N];
《数值分析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 水平位移监测方案 一、精度选择 按照设计要求，对照《工程测量规范》（GB 50026-2007），选用三等水平位移监测网进行检测，可以满足精度要求。 表1-2 水平角方向观测法的技术指标 （1）观测原理：如下图所示，如需观测某方向上的水平位移PP′，在监测区域一定距离以外选定工作基点A，水平位移监测点的布设应尽量与工作基点在一条直线上。沿监测点与基准点连线方向在一定远处（100~200m）选定一个控制 （2）精度分析： 由小角法的观测原理可知，距离D和水平角β是两个相互独立的观测值，所以由上式根据误差传播定律可得水平位移的观测误差： 水平位移观测中误差的公式，表明： ①距离观测误差对水平位移观测误差影响甚微，一般情况下此部分误 差可以忽略不计，采用钢尺等一般方法量取即可满足要求； ②影响水平位移观测精度的主要因素是水平角观测精度，应尽量使用 高精度仪器或适当增加测回数来提高观测度； ③经纬仪的选用应根据建筑物的观测精度等级确定，在满足观测精度 要求的前提下，可以使用精度较低的仪器，以降低观测成本。 优点：此方法简单易行，便于实地操作，精度较高。 不足：须场地较为开阔，基准点应该离开监测区域一定的距离之外，设在不受施工影响的地方。 由此可知，对仪器测角精度的要求，取决于监测点距离站点的远近。距离越远，则要求测角精度越高。根据现场踏勘布点，最远监测点距离站点不超过50m，对照《工程测量规范》，选用三等或四等水平位移监测网进行检测，可以满足精度要求。本次实习采用测小角法测量三等水平位移监测网进行检测。 二、作业流程 1.选点选取两个监测点P1，P2、一个测站点（工作基点）A、一个后视点B。 2.观测按照测回法水平角观测水平夹角。在A点安置全站仪，在B点和P1,P2点设置瞄准标志，按下列步骤进行测回法水平角观测。 （1）在全站仪盘左位置瞄准目标B，将度盘置零，读得水平度盘读数并记录。（2）瞄准目标P1，读得水平度盘读数并记录。盘左位置测得半测回水平角。（3）倒转望远镜成盘右位置，瞄准目标B，将度盘置零，读得水平度盘读数并记录。 （4）瞄准目标P1，读得水平度盘读数并记录。盘右位置测得半测回水平角。（5）用盘左、盘右两个位置观测水平角取平均值作为一测回水平角观测的结果。 《数值分析》计算实习题目 第一题： 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 水平位移监测方案 文稿归稿存档编号：[KKUY-KKIO69-OTM243-OLUI129-G00I-FDQS58- 水平位移监测方案 一、精度选择 按照设计要求，对照《工程测量规范》（GB 50026-2007），选用三等水平位移监测网进行检测，可以满足精度要求。 表1-1 水平位移基准网的主要技术指标 表1-2 水平角方向观测法的技术指标 （1）观测原理：如下图所示，如需观测某方向上的水平位移PP′，在监测区域一定距离以外选定工作基点A，水平位移监测点的布设应尽量与工作基点在一条直线上。沿监测点与基准点连线方向在一定远处（100~200m）选定一个控制点B，作为零方向。在B点安置觇牌，用测回法观测水平角BAP，测定一段时间内观测点与基准点连线与零方向间角度变化值，根据δ=△β*D/ρ（式中D为观测点P至工作基点A的距离，ρ=206265）计算水平位移。 （2）精度分析： 由小角法的观测原理可知，距离D和水平角β是两个相互独立的观测值，所以由上式根据误差传播定律可得水平位移的观测误差： 水平位移观测中误差的公式，表明： ①距离观测误差对水平位移观测误差影响甚微，一般情况下此部分误差可以忽 略不计，采用钢尺等一般方法量取即可满足要求； ②影响水平位移观测精度的主要因素是水平角观测精度，应尽量使用高精度仪 器或适当增加测回数来提高观测度； ③经纬仪的选用应根据建筑物的观测精度等级确定，在满足观测精度要求的前 提下，可以使用精度较低的仪器，以降低观测成本。 优点：此方法简单易行，便于实地操作，精度较高。 不足：须场地较为开阔，基准点应该离开监测区域一定的距离之外，设在不受施工影响的地方。 由此可知，对仪器测角精度的要求，取决于监测点距离站点的远近。距离越远，则要求测角精度越高。根据现场踏勘布点，最远监测点距离站点不超过50m，对照《工程测量规范》，选用三等或四等水平位移监测网进行检测，可以满足精度要求。本次实习采用测小角法测量三等水平位移监测网进行检测。 数值分析第三次大作业 一、算法的设计方案： （一）、总体方案设计： 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。 水平位移观测法、垂直位移观测法的种类,特点和适用条件 水平位移监测：对水工建筑物的顺水流方向或顺轴线方向的水平位移变化进行监测常用观测方法分两大类。一类是基准线法,基准线法是通过一条固定的基准线来测定监测点的位移,常见的有视准线法、引张线法、激光准直法、垂线法。 另一类是大地测量方法,大地测量方法主要是以外部变形监测控制网点为基准,以大地测量方法测定被监测点的大地坐标,进而计算被监测点的水平位移,常见的有交会法、精密导线法、三角测量法、GPS观测法等。 一、视准线法：通过视准线或经纬仪建立一个平行或通过坝轴线的铅直平面作为基准面，定期观测坝上测点与基准面之间偏离值的大小即为该点的水平位移。 适用于直线形混凝土闸坝顶部和坝面的水平位移观测。当采用这一方法时，主要的是要求它们的端点稳定，所以必须要作适当的布置，只能是定期地测定端点的位移值，而将观测值加以改正。视准线观测方法特点是速度快，精度较高，原理简单、方法实用、实施简便、投资较少的特点, 在水平位移观测中得到了广泛应用。 不足是对较长的视准线而言, 由于视线长, 使照准误差增大, 甚至可能造成照困难。当即准线太长时，目标模糊，照准精度太差且后视点与测点距离相差太远，望远镜调焦误差较大，无疑对观测成果有较大影响。 小角法：是水平位移监测中常用的方法，该方法最早应用于水库大坝的变形监测，其基本原理是一通过大坝轴线的固定不变的铅直平面为基准面，通过测定基准线方向之间的微小角度从而计算观测点相对予基准线的偏离值，根据偏离值在各观测周期中的变化确定位移量。由于所需测定的位移通常很细微，因此对位移的观测精度要求很高，需要采取各种提高观测精度的措施，观测过程中需要对各作业环节严格把握，哪怕仅仅是一个小环节的失误，都可能导致最终监测精度不能满足要求。 二、引张线法：利用张紧在两工作基点之间的不锈钢丝作为基准线，测量沿线测点和钢丝之间的相对位移，以确定该点的水平位移。 适用于大型直线形混凝土的廊道内测点的水平位移观测。主要用于测定混凝土建筑物垂直于轴线方向的（顺水流方向）水平位移。 活动觇牌法: 主要用于短距离视准线观测中，活动觇牌多用于水工建筑物、桥梁、码头和滑坡等水平位移观测，可满足坝内精密导线测量的近坝区水平位移监测网等各种场合的测量需要，活动觇标是被安置在位移标点上，供经纬仪照准，从而在觇标的游标尺上读出位移标点的偏离值。主要特点传动灵活、隙动差小，可精确到0.1mm . 《数值分析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**+==。 航空航天大学 《数值分析》计算实习报告 第一大题 学院：自动化科学与电气工程学院 专业：控制科学与工程 学生姓名： 学号： 教师： 电话： 完成日期： 2015年11月6日 航空航天大学 Beijing University of Aeronautics and Astronautics 实习题目： 第一题 设有501501?的实对称矩阵A ， ??? ???? ?????????=5011A a b c b c c b c b a 其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1 .0-==???=--=c b i e i i a i i 。矩阵A 的特征值为)501,,2,1(???=i i λ，并且有 ||min ||,501 150121i i s λλλλλ≤≤=≤???≤≤ 1.求1λ，501λ和s λ的值。 2.求A 的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,,2,1(???=k k i λ。 3.求A 的(谱数)条件数2)A (cond 和行列式detA 。 说明： 1.在所用的算法中，凡是要给出精度水平ε的，都取12-10=ε。 2.选择算法时，应使矩阵A 的所有零元素都不储存。 3.打印以下容： （1）全部源程序； （2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。 4.采用e 型输出实型数，并且至少显示12位有效数字。 一、算法设计方案 1、求1λ，501λ和s λ的值。 由于||min ||,501 150121i i s λλλλλ≤≤=≤???≤≤，可知绝对值最大特征值必为1λ和501 λ其中之一，故可用幂法求出绝对值最大的特征值λ，如果λ=0，则1λ=λ，否则 501λ=λ。将矩阵A 进行一下平移： I -A A'λ= （1） 对'A 用幂法求出其绝对值最大的特征值'λ，则A 的另一端点特征值1λ或501λ为'λ+λ。 s λ为按模最小特征值，||min ||501 1i i s λλ≤≤=，可对A 使用反幂法求得。 2、求A 的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,...,2,1(=k k i λ。 计算1)1,2,...,50=(i i λ-k μ，其模值最小的值对应的特征值k λ与k μ最接近。因此对A 进行平移变换： )39,,2,1k -A A k k ==（I μ （2） 对k A 用反幂法求得其模最小的特征值'k λ，则k λ='k λ+k μ。 3、求A 的(谱数)条件数2)(A cond 和行列式detA 。 由矩阵A 为非奇异对称矩阵可得： | | )(min max 2λλ=A cond （3） 其中max λ为按模最大特征值，min λ为按模最小特征值，通过第一问我们求得的λ和s λ可以很容易求得A 的条件数。 在进行反幂法求解时，要对A 进行LU 分解得到。因L 为单位下三角阵，行 列式为1，U 为上三角阵，行列式为主对角线乘积，所以A 的行列式等于U 的行列式，为U 的主对角线的乘积。 支架式水塔水平位移的实用简化计算 资讯类型：技术资料加入时间：2008年6月4日14:18 摘要对水塔进行动力分析时,可简化为单自由度体系。其基本周期可据在水塔水箱重心处单位水平力作用下该处的水平位移δ,按“顶点位移法”来计算。然而δ的计算至今缺少便于设计者应用的手算简化方法,为此本文提出了一种确定支架式水塔δ的简化计算模型及相应公式。公式形式简单,物理意义明确,便于计算。通过具体算例,采用本文方法与用su-persap程序作三维有限元计算对比,两者结果十分接近。 关键词支架式水塔水平位移基本周期简化计算 引言 笔者曾受委托对某建筑工地施工振动对邻近一座支架式水塔的影响进行安全性评估,需要及 时确定该水塔的自振周期,而我国抗震规范[1][2]及有关文献[3]尚未提供类似于筒壁式水塔或烟囱基本周期的计算公式。因此,深感即使在计算机十分普及的情况下,对于一些广泛应用的典型的建、构筑物,提出便捷而可行的简化计算方法,对工程设计仍具有较大的现实意义。水塔属一种高柔构筑物,其质量主要集中于顶部。在动力分析中,通常可以简化为单自由度系统,其基本周期则是一关键特征值。基本周期可由所谓的“顶点位移法”得出:t =2πmeqδ(1) 式中:meq为单自由度体系质量上的等效质量,通常可由下式确定:meq= m0+mt/4(2)m0,mt分别为顶部水箱及塔身的质量;δ为单位力作用下水塔水箱重心处所产生的水平位移。由于支架式塔身为空间格构式结构,且带有一定倾斜度,δ值计算是相当复杂的,至今尚未有便于设计者应用的简化手算方法。为此,本文通过对若干6根支柱水塔进行三维有限元计算,分析了水塔结构内力及位移的本质规律,在此基础上提出一个简便的计算模型,得到了确定δ值精度较高的手算公式,进而解决了基本周期的计算问题。与三维有限元分析结果十分吻合,可供这类水塔结构进行抗震分析,并与现行有关规范配合应用。 一、计算δ的简化模型 1?三维有限元分析主要规律 某典型的6根支柱水塔如图1所示,顶部作用一水平单位力p=1。经过对5个不同尺寸的这类水塔作三维有限元分析,得到如下主要规律: (1)静力分析,x方向作用p=1在各层x方向产生位移,与y方向上作用p=1在各层y方向上产生位移相等。动力分析,前2阶频率相同,分别属x方向及y方向的1阶振型。 (2)各层不同方位柱均存在反弯点,中间各层,反弯点基本上位于中点,而底层反弯点偏上,顶层反弯点偏下,但不及2/3处。 (3)立柱具有一定倾斜度,有效地减小了立柱中的弯矩和剪力,塔身剪力下部小,上部大,沿高度呈线性变化。 (4)各层圈梁中诸横梁受力情况如下: ①p=1沿x方向作用,各横梁主要弯矩位于竖向平面(绕水平轴),侧向弯矩及扭矩均很小,对主要弯矩,各梁都存在反弯点,且基本位于中点。主要弯矩从数值上看,梁中梁端弯矩恰为上述二梁的一半。 ②p=1沿y方向作用,各横梁主要弯矩也是位于竖向平面,侧向弯矩及扭矩较小,约比主要弯矩小一个数量级。梁中反弯点位于中点,而梁,中弯矩为常数,不沿长度变化,即不存在反弯点,且数值仅为上述4根梁梁端弯矩的1/10。 2?δ的简化计算模型 参照上述三维有限元分析所得规律,本文提:eici为i层柱的当量抗弯刚度;ici为各单柱截面绕自身形心轴惯性矩投影到塔身截面计算主轴上的代数和,如p=1沿x方向作用,参照图1,即考虑对y-y轴惯性矩,有对柱1:i1=i-y=bh3/12(b与h 分别为截面的宽度与高度)。对柱2:i2=i-ycos2α+i-xsin2α=bh312cos2α+hb312sin2α对6根柱式支架α=60°,则层间柱当量惯性矩为ic=2i1+4i2=3(i-y+i-x)=bh(h2+b2)/4:ri,r′i分别为第i层圈梁高程处及反弯点高程处相应的半径;hdi,hui分别为第i层反弯点之下、之上到圈梁的距离;ηeib为横梁当量抗弯刚度;ib为实际单根横梁截面惯性矩;η为简化模型中的当量系数。对η本文按如下原则定出:简化模型(图2)中横梁对位移计算的贡献为:δδ′=ri6ηeib(-qi+1hdi+1+-qihui)(hui+hdi+1) 《数值分析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、对矩阵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；最后将所有的实特征值的特征向量输出。 一、算法设计方案 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 ξ???? ? ?=?????? ， 浅谈几种水平位移的方法 【摘要：】本文对常用的几种水平位移的观测方法进行了比较系统的分析和比较，列出了这几种方法的原理，精度分析，优点以及不足，他们适用的场合等内容，对于在生产实践中进行水平位移观测时进行方法的选取具有一定的指导价值。 【关键字：】水平位移，视准线法，测小角法，前方交会法，极坐标法，反演小角法 当要观测某一特定方向（譬如垂直于基坑维护体方向）的位移时，经常采用视准线法、小角度法等观测方法。但当变形体附近难以找到合适的工作基点或需同时观测变形体两个方向位移时，则一般采用前方交会法。水平位移观测观测实践中利用较多的前方交会法主要有两种：测边前方交会法和测角前方交会法。另外还有极坐标法以及一些困难条件下的水平位移观测方法。 视准线法： 当需要测定变形体某一特定方向（譬如垂直于基坑维护体方向）的位移时，常使用视准线法或测小角法。 另外此方法还受到大气折光等因素的影响。 优点： 视准线观测方法因其原理简单、方法实用、实施简便、投资较少的特点, 在水平位移观测中得到了广泛应用，并且派生出了多种多样的观测方法，如分段视准线，终点设站视准线等。 不足： 对较长的视准线而言, 由于视线长, 使照准误差增大, 甚至可能造成照准困难。当即准线太长时，目标模糊，照准精度太差且后视点与测点距离相差太远，望远镜调焦误差较大，无疑对观测成果有较大影响。精度低，不易实现自动观测，受外界条件影响较大，而且变形值（位移标点的位移量）不能超出该系统的最大偏距值，否则无法进行观测。 测小角法： 当需要测定变形体某一特定方向（譬如垂直于基坑维护体方向）的位移时，常使用视准线法或小角度法 原理：如下图所示，如需观测某方向上的水平位移PP′，在监测区域一定距离以外选定工作基点A，水平位移监测点的布设应尽量与工作基点在一条直线上。沿监测点与基准点连线方向在一定远处（100~200m）选定一个控制点B，作为零方向。在B 水平位移几种监测方法的分析和比较 【摘要：】本文对常用的几种水平位移的观测方法进行了比较系统的分析和比较，列出了这几种方法的原理，精度分析，优点以及不足，他们适用的场合等内容，对于在生产实践中进行水平位移观测时进行方法的选取具有一定的指导价值。 【关键字：】水平位移，视准线法，测小角法，前方交会法，极坐标法，反演小角法 当要观测某一特定方向（譬如垂直于基坑维护体方向）的位移时，经常采用视准线法、小角度法等观测方法。但当变形体附近难以找到合适的工作基点或需同时观测变形体两个方向位移时，则一般采用前方交会法。水平位移观测观测实践中利用较多的前方交会法主要有两种：测边前方交会法和测角前方交会法。另外还有极坐标法以及一些困难条件下的水平位移观测方法。 视准线法： 当需要测定变形体某一特定方向（譬如垂直于基坑维护体方向）的位移时，常使用视准线法或测小角法。 原理：如下图所示，点A、B是视准线的两个基准点（端点），1、2、3为水平位移观测点。观测时将经纬仪置于A点，将仪器照准B点，将水平制动装置制动。竖直转动经纬仪，分别转至1、2、3 三个点附近，用钢尺等工具测得水准观测点至A—B这条视准线的距离。根据前后两次的测量距离，得出这段时间内水平位移量。 精度分析： 由基准线的设置过程可知，观测误差主要包括仪器测站点仪器对中误差，视准线照准误差，读数照准误差，其中，影响最大的无疑是读数照准误差。 可知，当即准线太长时，目标模糊，读数照准精度太差；且后视点与测点距离相差太远，望远镜调焦误差较大，无疑对观测成果有较大影响。 另外此方法还受到大气折光等因素的影响。 优点： 视准线观测方法因其原理简单、方法实用、实施简便、投资较少的特点, 在水平位移观测中得到了广泛应用，并且派生出了多种多样的观测方法，如分段视准线，终点设站视准线等。 北航数值分析计算实习 报告一 Company number：【0089WT-8898YT-W8CCB-BUUT-202108】 北京航空航天大学 《数值分析》计算实习报告 第一大题 学 院：自动化科学与电气工程学院 专 业： 控制科学与工程 学 生 姓 名： 学 号： 教 师： 电 话： 完 成 日 期： 2015年11月6日 北京航空航天大学 Beijing University of Aeronautics and Astronautics 实习题目： 第一题 设有501501?的实对称矩阵A ， 其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==???=--=c b i e i i a i i 。矩阵A 的特征值为)501,,2,1(???=i i λ，并且有 1.求1λ，501λ和s λ的值。 2.求A 的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,,2,1(???=k k i λ。 3.求A 的(谱范数)条件数2)A (cond 和行列式detA 。 说明： 1.在所用的算法中，凡是要给出精度水平ε的，都取12-10=ε。 2.选择算法时，应使矩阵A 的所有零元素都不储存。 3.打印以下内容： （1）全部源程序； （2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。 4.采用e 型输出实型数，并且至少显示12位有效数字。 一、算法设计方案 1、求1λ，501λ和s λ的值。 由于||min ||,501 150121i i s λλλλλ≤≤=≤???≤≤，可知绝对值最大特征值必为1λ和501λ其中之 一，故可用幂法求出绝对值最大的特征值λ，如果λ=0，则1λ=λ，否则501λ=λ。将矩阵A 进行一下平移： I -A A'λ= （1） 对'A 用幂法求出其绝对值最大的特征值'λ，则A 的另一端点特征值1λ或501λ为 'λ+λ。 s λ为按模最小特征值，||min ||501 1i i s λλ≤≤=，可对A 使用反幂法求得。 2、求A 的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,...,2,1(=k k i λ。 计算1)1,2,...,50=(i i λ-k μ，其模值最小的值对应的特征值k λ与k μ最接近。因此对A 进行平移变换： )39,,2,1k -A A k k ==（I μ （2） 对k A 用反幂法求得其模最小的特征值'k λ，则k λ='k λ+k μ。水平位移监测方案
北航数值分析大作业第一题幂法与反幂法
水平位移监测方案
北航数值分析报告第三次大作业
北航数值分析大作业第二题精解
北航数值分析大作业第二题
水平位移观测法垂直位移观测法的种类_特点和适用条件(仅供参考版)
北航数值分析第二次大作业--QR分解
北航数值分析报告大作业第八题
北航数值分析计算实习报告一
水塔水平位移的计算
北航数值分析课程第一次大作业讲解
北航数值分析大作业第二次
北航数值分析大作业3
浅谈水平位移的几种方法
水平位移几种监测方法
北航数值分析计算实习报告一
相关主题
文本预览