下载文档原格式
(单选题)1: 当输入信号的复频率等于系统函数的零点时，系统的强迫响应分量为（）。
A: 无穷大
B: 不为零的常数
C: 0
D: 随输入信号而定
正确答案: C
(单选题)2: 满足傅氏级数收敛条件时，周期信号f(t)的平均功率（）。
A: 大于各谐波分量平均功率之和
B: 不等于各谐波分量平均功率之和
C: 小于各谐波分量平均功率之和
D: 等于各谐波分量平均功率之和
正确答案: D
(单选题)3: 卷积δ(t)*f(t)*δ(t)的结果为（）。
A: δ(t)
B: δ(2t)
C: f(t)
D: f(2t)
正确答案: C
(单选题)4: 信号的时宽与信号的频宽之间呈（）。
A: 正比关系
B: 反比关系
C: 平方关系
D: 没有关系
正确答案: B
(单选题)5: 设一个矩形脉冲的面积为S，则矩形脉冲的傅氏变换在原点处的函数值等于（）。
A: S/2
B: S/3
C: S/4
D: S
正确答案: D
(单选题)6: 线性系统具有（）。
A: 分解特性
B: 零状态线性
C: 零输入线性
D: 以上全对
正确答案: D
(单选题)7: 如果一连续时间二阶系统的系统函数H(s)的共轭极点在虚轴上，则它的h(t)应是（）。
基于内聚力模型(CZM)的单筋拉拔数值分析方法研究景剑;强峰;施凯【摘要】目前化学植筋粘结性能数值模拟中界面单元均以双弹簧单元为主,但是模拟结果与相应的试验结果有较大偏差.为了提高数值模拟的精确度,本文基于双线性内聚力模型(Cohesive Zone Model)进行了单筋拉拔试验的有限元模拟,采用双线性应力-张开位移模型定义内聚力单元本构关系,进行了参数分析,得到了内聚力参数对计算结果的影响规律,并对一些试验的荷载-位移曲线进行参数拟合以确定合理参数,从而验证了该植筋模拟方法的有效性.【期刊名称】《低温建筑技术》【年(卷),期】2018(040)007【总页数】5页(P57-60,64)【关键词】化学植筋;单筋拉拔试验;内聚力模型;参数分析【作者】景剑;强峰;施凯【作者单位】江苏省建筑工程质量检测中心有限公司,南京210008;河海大学土木与交通学院,南京210098;江苏省建筑工程质量检测中心有限公司,南京210008【正文语种】中文【中图分类】TU5020 引言化学植筋是目前加固改造领域应用相当广泛的后锚固连接技术，现有植筋承载力和力学性能的相关研究大多限于单筋拉拔试验研究，由于拉拔试验试件制作及试验装置比较简单，试验结果便于分析，长期以来一直作为研究化学植筋性能的有效方法，但是由于试验中存在诸多不确定性因素，如果通过大量的拉拔试验研究化学植筋性能，不仅耗费过多的试验材料，而且需要很长的试验周期，给研究带来诸多不便。
近些年来，应用有限元分析方法进行化学植筋锚固性能研究已成为一种方便有效的数值模拟方法。
在早期的植筋锚固系统研究中，国内外同行已发表了一些有关粘结锚固的研究成果。
Cook等人[1]通过单筋拉拔试验总结出了在混凝土构件中，植筋的破坏模式，研究了单个钢筋锚固的破坏过程和机理，给出了单筋的粘结锚固建议；郭晓飞[2]提出了采用双弹簧单元模拟混凝土与植筋胶界面单元和钢筋与植筋胶界面单元，并考虑了植筋胶的厚度，采用四边形滑移单元对植筋胶进行模拟。
特征值法对元素为实数或复数的n×n矩阵A，求数λ和n维非零向量x使A x=λx，这样的问题称为代数特征值问题，也称矩阵特征值问题，λ和x分别称为矩阵A的特征值和特征向量。
代数特征值问题的数值解法是计算数学的主要研究课题之一，它常出现于动力系统和结构系统的振动问题中。
在常微分方程和偏微分方程的数值分析中确定连续问题的近似特征系，若用有限元方法或有限差分方法求解，最终也化成代数特征值问题。
此外，其他数值方法的理论分析，例如确定某些迭代法的收敛性条件和初值问题差分法的稳定性条件，以及讨论计算过程对舍入误差的稳定性问题等都与特征值问题有密切联系。
求解矩阵特征值问题已有不少有效而可靠的方法。
矩阵A的特征值是它的特征多项式P n(λ)det(λI-A)的根，其中I为单位矩阵。
但阶数超过4的多项式一般不能用有限次运算求出根，因而特征值问题的计算方法本质上是迭代性质的，基本上可分为向量迭代法和变换方法两类。
向量迭代法是不破坏原矩阵A，而利用A对某些向量作运算产生迭代向量的求解方法，多用来求矩阵的部分极端特征值和相应的特征向量，特别适用于高阶稀疏矩阵。
乘幂法、反幂法都属此类，隆措什方法也常作为迭代法使用。
变换方法是利用一系列特殊的变换矩阵(初等下三角阵、豪斯霍尔德矩阵、平面旋转矩阵等)，从矩阵A出发逐次进行相似变换，使变换后的矩阵序列趋于容易求得特征值的特殊形式的矩阵(对角阵、三角阵、拟三角阵等)；多用于求解全部特征值问题，其优点是收敛速度快，计算结果可靠，但由于原矩阵A被破坏，当A是稀疏矩阵时，在计算过程中很难保持它的稀疏性，因而大多数变换方法只适于求解中小规模稠密矩阵的全部特征值问题。
雅可比方法、吉文斯－豪斯霍尔德方法以及LR方法、QR方法等都属此类。
乘幂法计算矩阵的按模最大的特征值及对应特征向量的一种向量迭代法。
设A为具有线性初等因子的矩阵，它的n个线性无关的特征向量是u i(i=1，2，…，n)，特征值排列次序满足是一个n维非零向量，于是若λ1>λ2，则当α1≠0，且k足够大时，A k z0除相差一个纯量因子外趋于λ1所对应的特征向量，这就是乘幂法的基本思想。
北航《信号与系统》在线作业一一、单选题：1.理想低通滤波器是( )( )。
(满分:3)A. 因果系统B. 物理可实现系统C. 非因果系统D. 响应不超前于激励发生的系统正确答案:C2.理想低通滤波器一定是( )( )。
(满分:3)A. 稳定的物理可实现系统B. 稳定的物理不可实现系统C. 不稳定的物理可实现系统D. 不稳定的物理不可实现系统正确答案:B3.因果系统是物理上( )( )( )系统。
(满分:3)A. 不可实现的B. 可实现的C. 未定义的D. 以上都不对正确答案:B4.哪种滤波器功能是只允许信号中的低频成分通过( )( )。
(满分:3)A. 理想低通滤波器B. 带通滤波器C. 高通滤波器D. 以上全对正确答案:A5.信号〔ε(t)－ε(t－2)〕的拉氏变换的收敛域为( )( )。
(满分:3)A. Re[s]>0B. Re[s]>2C. 全S平面D. 不存在正确答案:C6.激励为零，仅由系统的( )( )引起的响应叫做系统的零输入响应。
(满分:3)A. 初始状态B. 中间状态C. 最终状态D. 以上全不对正确答案:A7.如果一连续时间二阶系统的系统函数H(s)的共轭极点在虚轴上，则它的h(t)应是( )( )。
(满分:3)A. 指数增长信号B. 指数衰减振荡信号C. 常数D. 等幅振荡信号正确答案:D8.欲使信号通过系统后只产生相位变化，则该系统一定是( )( )。
(满分:3)A. 高通滤波网络B. 带通滤波网络C. 全通网络D. 最小相移网络正确答案:C9.已知某连续时间系统的系统函数H（s）= 1/(s+1)，该系统属于什么类型( )( )。
(满分:3)A. 高通滤波器B. 低通滤波器C. 带通滤波器D. 带阻滤波器正确答案:B10.一信号x(t)的最高频率为500Hz，则利用冲激串采样得到的采样信号x(nT)能唯一表示出原信号的最大采样周期为( )( )。
(满分:3)A. 500B. 1000C. 0.05D. 0.001正确答案:D二、多选题：1.偶周期信号的傅氏级数中只有( )( )。
传热学（一）第一部分选择题•单项选择题（本大题共 10 小题，每小题 2 分，共 20 分）在每小题列出的四个选项中只有一个选项是符合题目要求的，请将正确选项前的字母填在题后的括号内。
1. 在稳态导热中 , 决定物体内温度分布的是 ( )A. 导温系数B. 导热系数C. 传热系数D. 密度2. 下列哪个准则数反映了流体物性对对流换热的影响 ?( )A. 雷诺数B. 雷利数C. 普朗特数D. 努谢尔特数3. 单位面积的导热热阻单位为 ( )A. B. C. D.4. 绝大多数情况下强制对流时的对流换热系数 ( ) 自然对流。
A. 小于B. 等于C. 大于D. 无法比较5. 对流换热系数为 100 、温度为 20 ℃的空气流经 50 ℃的壁面，其对流换热的热流密度为（） A. B. C. D.6. 流体分别在较长的粗管和细管内作强制紊流对流换热，如果流速等条件相同，则（）A. 粗管和细管的相同B. 粗管内的大C. 细管内的大D. 无法比较7. 在相同的进出口温度条件下，逆流和顺流的平均温差的关系为（）A. 逆流大于顺流B. 顺流大于逆流C. 两者相等D. 无法比较8. 单位时间内离开单位表面积的总辐射能为该表面的（）A. 有效辐射B. 辐射力C. 反射辐射D. 黑度9. （）是在相同温度条件下辐射能力最强的物体。
A. 灰体B. 磨光玻璃C. 涂料D. 黑体10. 削弱辐射换热的有效方法是加遮热板，而遮热板表面的黑度应（）A. 大一点好B. 小一点好C. 大、小都一样D. 无法判断第二部分非选择题•填空题（本大题共 10 小题，每小题 2 分，共 20 分）11. 如果温度场随时间变化，则为。
12. 一般来说，紊流时的对流换热强度要比层流时。
13. 导热微分方程式的主要作用是确定。
14. 当 d 50 时，要考虑入口段对整个管道平均对流换热系数的影响。
15. 一般来说，顺排管束的平均对流换热系数要比叉排时。
目录1等参单元及其应用 (1)1.1概述 (1)1.1.1等参单元的概念、原理 (1)1.1.2等参单元对有限元法工程应用的意义 (1)1.2等参单元的数值积分方法 (1)1.2.1等参单元刚度矩阵的数值积分方法 (1)1.2.2确定积分阶的原理 (2)1.2.3全积分单元与减缩积分单元讨论和评价 (3)1.3线性等参单元 (3)1.3.1全积分、减缩积分线性等参单元有关问题的分析讨论 (3)1.4等参单元的应用 (5)2分析与计算 (6)2.1四节点平面等参单元的收敛协调性 (6)2.2八节点平面等参单元 (8)2.33节点平面三角形单元 (9)2.420节点六面体等参单元 (10)2.520节点六面体等参单元 (11)3上机实验 (15)3.1实验一 (15)3.1.1实验题目 (15)3.1.2实验目的 (15)3.1.3建模概述 (15)3.1.4计算结果分析与结论 (16)3.1.5实验体会与总结 (32)3.2实验二 (33)3.2.1实验题目 (33)3.2.2实验目的 (33)3.2.3建模概述 (33)3.2.4计算结果分析与讨论 (34)3.2.5实验体会与总结 (36)3.3实验三 (36)3.3.1实验题目 (36)3.3.2实验目的 (36)3.3.3建模概述 (37)3.3.4计算结果分析与结论 (37)3.3.5实验体会与总结 (44)1 等参单元及其应用1.1 概述1.1.1 等参单元的概念、原理普通单元受到两个方面的限制：（1）单元的精度。
单元的节点数越多，单元精度越高；（2）单元几何上的限制。
普通矩形和六面体单元都不能模拟任意形状几何体，所有几种普通单元都是直线边界，处理曲边界几何体误差较大。
为了解决上述矛盾，方法就是突破矩形单元和六面体单元几何方面的限制，使其成为任意四边形和任意六面体单元，这类单元位移模式和形函数的构造和单元列式的导出不能沿用构造简单单元的方法，必须引入等参变换，采用相同的插值函数对单元的节点坐标和节点位移在单元上进行插值。
北航《运筹学》在线作业3一、单选题（共 10 道试题，共 30 分。
）1. 基本可行解是满足非负条件的基本解。
（）A. 正确B. 错误C. 不一定D. 无法判断正确答案：A2. 不满足匈牙利法的条件是A. 问题求最小值B. 效率矩阵的元素非负C. 人数与工作数相等D. 问题求最大值正确答案：D3. 连通图G有n个点，其部分树是T，则有A. T有n个点n条边B. T的长度等于G的每条边的长度之和C. T有n个点n－1条边D. T有n－1个点n条边正确答案：C4. 线性规划标准型中，决策变量（）是非负的。
A. 一定B. 一定不C. 不一定D. 无法判断正确答案：A5. 运输问题可以用( )法求解。
A. 定量预测B. 单纯形C. 求解线性规划的图解D. 关键线路正确答案：B6. 下列说法正确的是A. 割集是子图B. 割量等于割集中弧的流量之和C. 割量大于等于最大流量D. 割量小于等于最大流量正确答案：C7. 动态规划求解的一般方法是什么？（）A. 图解法B. 单纯形法C. 逆序求解D. 标号法正确答案：C8. 运输问题的数学模型属于A. 0-1规划模型B. 整数规划模型C. 网络模型D. 以上模型都是正确答案：C9.通过什么方法或者技巧可以把产销不平衡运输问题转化为产销平衡运输问题( )A. 非线性问题的线性化技巧B. 静态问题的动态处理C. 引入虚拟产地或者销地D. 引入人工变量正确答案：C10. 下列结论正确的有A.运输问题的运价表第r行的每个Cij同时加上一个非零常数k，其最优调运方案不变B.运输问题的运价表第p列的每个Cij同时乘以一个非零常数k，其最优调运方案不变C. 运输问题的运价表的所有Cij同时乘以一个非零常数k,其最优调运方案变化D. 不平衡运输问题不一定存在最优解正确答案：A北航《运筹学》在线作业3二、多选题（共 10 道试题，共 40 分。
）1. 下列结论不正确的有A.运输问题的运价表第r行的每个Cij同时加上一个非零常数k，其最优调运方案不变B.运输问题的运价表第p列的每个Cij同时乘以一个非零常数k，其最优调运方案不变C. 运输问题的运价表的所有Cij同时乘以一个非零常数k,其最优调运方案变化D. 不平衡运输问题不一定存在最优解正确答案：BCD2. 下面命题不正确的是（）。
一、算法设计方案：按题目要求，本程序运用带双步位移的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.txtpReadFile = 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;elseB[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){//计算drsum = 0;for (i = r+1; i < N; i++){sum += a[i][r] * a[i][r];}dr = sqrt(sum);//计算crif (a[r+1][r] > 0)cr = -dr;elsecr = dr;//计算hrhr = cr * cr - cr * a[r+1][r];//取向量ur[]for (i = 0; i < N; i++){if (i < r+1)ur[i] = 0;elseif (i == r+1)ur[i] = a[i][r] - cr;elseur[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;}//计算trsum = 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];elsewr[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];//取矩阵R1[][]为A[][]for (i = 0; i < N; i++)for (j = 0; j < N; j++)R[i][j] = A[i][j];//取矩阵Q1[][]为单位矩阵for (i = 0; i < N; i++)for (j = 0; j < N; j++){if (i == j)Q[i][j] = 1;elseQ[i][j] = 0;}//for (r = 0; r < N-1; r++){//判断R[i][r](i=r+1,r+2,...,N-1)是否全为零sum = 0;//变量sum使用前清零for (i = r+1; i < N; i++){sum = sum || R[i][r];}if (sum){//计算drsum = 0;for (i = r; i < N; i++){sum += R[i][r] * R[i][r];}dr = sqrt(sum);//计算crif (R[r][r] > 0)cr = -dr;elsecr = dr;//计算hrhr = cr * cr - cr * R[r][r];//取向量ur[]for (i = 0; i < N; i++){if (i < r)ur[i] = 0;elseif (i == r)ur[i] = R[i][r] - cr;elseur[i] = R[i][r];}//计算向量wr[]for (i = 0; i < N; i++){sum = 0;for (j = r; j < N; j++)sum += Q[i][j] * ur[j];wr[i] = sum;}//计算新的矩阵Qr+1[][],存储在Q[][]里面for (i = 0; i < N; i++)for (j = 0; j < N; j++)Q[i][j] = Q[i][j] - wr[i] * ur[j] / hr;//计算向量pr[]for (i = 0; i < N; i++){sum = 0;for (j = r; j < N; j++)sum += R[j][i] * ur[j];pr[i] = sum / hr;}//计算新产生的R[][]for (i = 0; i < N; i++)for (j = 0; j < N; j++)R[i][j] = R[i][j] - ur[i] * pr[j];}}}/*************************************Mk[][]矩阵的QR分解及求Ak+1[][]的计算形参m为每次降阶之后的值无返回值**************************************/void QR_decompositionMk(double mk[][N],int m, double a[][N]){ intr,i,j;double tr,hr,cr,dr,sum;double ur[N],pr[N],qr[N],wr[N],vr[N],br[N][N],Cr[N][N];//取矩阵br[][]for (i = 0; i <= m; i++)for (j = 0; j <= m; j++)br[i][j] = mk[i][j];//取矩阵Cr[][]for (i = 0; i <= m; i++)for (j = 0; j <= m; j++)Cr[i][j] = a[i][j];//for (r = 0; r <= m-1; r++){//判断br[i][r](i=r+1,r+2,...,m)是否全为零sum = 0;//变量sum使用前清零for (i = r+1; i <= m; i++){sum = sum || br[i][r];}if (sum){//计算drsum = 0;for (i = r; i <= m; i++){sum += br[i][r] * br[i][r];}dr = sqrt(sum);//计算crif (br[r][r] > 0)cr = -dr;elsecr = dr;//计算hrhr = cr * cr - cr * br[r][r];//取向量ur[]for (i = 0; i <= m; i++){if (i < r)ur[i] = 0;elseif (i == r)ur[i] = br[i][r] - cr;elseur[i] = br[i][r];}//计算向量vr[]for (i = 0; i <= m; i++){sum = 0;for (j = r; j <= m; j++)sum += br[j][i] * ur[j];vr[i] = sum / hr;}//计算新的矩阵Br+1[][],存储在Br[][]里面for (i = 0; i <= m; i++)for (j = 0; j <= m; j++)br[i][j] = br[i][j] - ur[i] * vr[j];//计算向量pr[]for (i = 0; i <= m; i++){sum = 0;for (j = r; j <= m; j++)sum += Cr[j][i] * ur[j];pr[i] = sum / hr;}//计算向量qr[]for (i = 0; i <= m; i++){sum = 0;for (j = r; j <= m; j++)sum += Cr[i][j] * ur[j];qr[i] = sum / hr;}//计算trsum = 0;for (i = r; i <= m; i++)sum += pr[i] * ur[i];tr = sum / hr;//计算向量wr[]for (i = 0; i <= m; i++){if (i < r)wr[i] = qr[i];elsewr[i] = qr[i] - tr * ur[i];}//计算新产生的矩阵Cr[][]for (i = 0; i <= m; i++)for (j = 0; j <= m; j++)Cr[i][j] = Cr[i][j] - wr[i] * ur[j] - ur[i] * pr[j]; }}//将计算出的Cr[][]矩阵中的值赋给矩阵A[][]，得到新的矩阵。
