北航研究生 算法设计与分析 Assignment_2

用分支定界算法求以下问题：

某公司于乙城市的销售点急需一批成品，该公司成品生产基地在甲城市。甲城市与乙城市之间共有n 座城市，互

相以公路连通。甲城市、乙城市以及其它各城市之间的公路连通情况及每段公路的长度由矩阵M1 给出。每段公路均

由地方政府收取不同额度的养路费等费用，具体数额由矩阵M2 给出。

请给出在需付养路费总额不超过1500 的情况下，该公司货车运送其产品从甲城市到乙城市的最短运送路线。

具体数据参见文件：

m1.txt: 各城市之间的公路连通情况及每段公路的长度矩阵(有向图); 甲城市为城市Num.1，乙城市为城市Num.50。

m2.txt: 每段公路收取的费用矩阵（非对称）。

思想：利用Floyd算法的基本方法求解。

程序实现流程说明：

1.将m1.txt和m

2.txt的数据读入两个50×50的数组。

2.用Floyd算法求出所有点对之间的最短路径长度和最小费用。

3.建立一个堆栈，初始化该堆栈。

4.取出栈顶的结点，检查它的相邻的所有结点，确定下一个当前最优路径上的结点，被扩展的结点依次

加入堆栈中。在检查的过程中，如果发现超出最短路径长度或者最小费用，则进行”剪枝”，然后回溯。

5.找到一个解后，保存改解，然后重复步骤4。

6.重复步骤4、5，直到堆栈为空，当前保存的解即为最优解。

时间复杂度分析：

Floyd算法的时间复杂度为3

()

O N，N为所有城市的个数。

该算法的时间复杂度等于DFS的时间复杂度，即O(N+E)。其中，E为所有城市构成的有向连通图的边的

总数。但是因为采用了剪枝，会使实际运行情况的比较次数远小于E。

求解结果：

算法所得结果:

甲乙之间最短路线长度是:464

最短路线收取的费用是:1448

最短路径是:

1 3 8 11 15 21 23 26 3

2 37 39 45 47 50

C源代码（注意把m1.txt与m2.txt放到与源代码相同的目录下，下面代码可直接复制运行）：

#include

#include

#include

#include

#define N 50

#define MAX 52

void input(int a[N][N],int b[N][N]);

void Floyd(int d[N][N]);

void fenzhi(int m1[N][N],int m2[N][N],int mindist[N][N],int mincost[N][N]);

int visited[N],bestPath[N];

void main()

{

clock_t start,finish;

double duration;

int i,j,mindist[N][N],mincost[N][N],m1[N][N],m2[N][N]; /* m1[N][N]和m2[N][N]分别代表题目所给的距离矩阵和代价矩阵*/

// int visited[N],bestPath[N];

FILE *fp,*fw;

// system("cls");

time_t ttime;

time(&ttime);

printf("%s",ctime(&ttime));

start=clock();

for(i=0;i

{

visited[i]=0;

bestPath[i]=0;

}

fp=fopen("m1.txt","r"); /* 把文件中的距离矩阵m1读入数组mindist[N][N] */

if(fp==NULL)

{

printf("can not open file\n");

return;

}

for(i=0;i

for(j=0;j

fscanf(fp,"%d",&mindist[i][j]);

fclose(fp); /* 距离矩阵m1读入完毕*/

fp=fopen("m2.txt","r"); /* 把文件中的代价矩阵m2读入数组mincost[N][N] */

if(fp==NULL)

{

printf("can not open file\n");

return;

}

for(i=0;i

for(j=0;j

fscanf(fp,"%d",&mincost[i][j]);

fclose(fp); /* 代价矩阵m2读入完毕*/

input(m1,mindist); /* mindist[N][N]赋值给m1[N][N]，m1[N][N]代表题目中的距离矩阵*/

input(m2,mincost); /* mincost[N][N]赋值给m2[N][N]，m2[N][N]代表题目中的代价矩阵*/

for(i=0;i

{

mindist[i][i]=9999;

mincost[i][i]=0;

}

Floyd(mindist); /* 用弗洛伊德算法求任意两城市之间的最短距离，结果存储在数组mindist[N][N]中*/

/*

fw=fopen("1.txt","w");

for(i=0;i

for(j=0;j

fprintf(fw,"%4d ",mindist[i][j]);

fprintf(fw,"\n");

}

fclose(fw);

// getchar();

//*/

Floyd(mincost); /* 用弗洛伊德算法求任意两城市之间的最小代价，结果存储在数组mincost[N][N]中*/

/*

fw=fopen("2.txt","w");

for(i=0;i

for(j=0;j

fprintf(fw,"%4d ",mincost[i][j]);

fprintf(fw,"\n");

}

fclose(fw);

// getchar();

//*/

fenzhi(m1,m2,mindist,mincost); /* 调用分支定界的实现函数，寻找出所有的可行路径并依次输出*/

finish=clock();

duration = (double)(finish - start) / CLOCKS_PER_SEC;

printf( "%f seconds\n", duration );

//*/

}

void Floyd(int d[N][N]) /* 弗洛伊德算法的实现函数*/

{

int v,w,u,i;

for(u=0;u

{

for(v=0;v

for(w=0;w

if(d[v][u]+d[u][w]

{

//printf("v,w,u,d[v][u],d[u][w],d[v][w] %d %d %d %d %d %d",v+1,w+1,u+1,d[v][u],d[u][w],d[v][w]);getchar();

d[v][w]=d[v][u]+d[u][w];

}

}

}

}

void input(int a[N][N],int b[N][N]) /* 把矩阵b赋值给矩阵a */

{

int i,j;

for(i=0;i

for(j=0;j

a[i][j]=b[i][j];

}

void fenzhi(int m1[N][N],int m2[N][N],int mindist[N][N],int mincost[N][N])

{

int stack[MAX],depth=0,next,i,j; /* 定义栈，depth表示栈顶指针；next指向每次遍历时当前所处城市的上一个已经遍历的城市*/

int bestLength,shortestDist,minimumCost,distBound=9999,costBound=9999;

int cur,currentDist=0,currentCost=0; /* cur指向当前所处城市，currentDist和currentCost分别表示从甲城市到当前所处城市的最短距离和最小代价，

currentDist和currentCost初值为0表示从甲城市出发开始深度优先搜索*/

stack[depth]=0; /* 对栈进行初始化*/

stack[depth+1]=0;

visited[0]=1; /* visited[0]=1用来标识从甲城市开始出发进行遍历，甲城市已被访问*/ while(depth>=0) /* 表示遍历开始和结束条件，开始时从甲城市出发，栈空，depth=0；结束时遍历完毕，所有节点均被出栈，故栈也为空，depth=0 */

/* 整个while（）循环体用来实现从当前的城市中寻找一个邻近的城市*/

{

cur=stack[depth]; /* 取栈顶节点赋值给cur，表示当前访问到第cur号城市*/

next=stack[depth+1]; /* next指向当前所处城市的上一个已经遍历的城市*/

for(i=next+1;i

{

if((currentCost+mincost[cur][N-1]>costBound)||(currentDist+mindist[cur][N-1]>=distBound)){ /* 所试探的城市满足剪枝条件，进行剪枝*/

//printf("here1 %d %d %d %d %d %d %d\n",cur,currentCost,mincost[cur][49],costBound,currentDist,mindist[cur][49],dist Bound); getchar();

// printf("%d %d %d %d %d %d",cur,i,m1[cur][i],currentCost,mincost[cur][49],costBound); getchar();

continue;

}

if(m1[cur][i]==9999) continue; /* 所试探的城市不连通*/

if(visited[i]==1) continue; /* 所试探的城市已被访问*/

if(i

if(i==N) /* 判断for循环是否是由于搜索完所有城市而终止的，如果是（i==N），进行回溯*/ {

// printf("here");getchar();

depth--;

currentDist-=m1[stack[depth]][stack[depth+1]];

currentCost-=m2[stack[depth]][stack[depth+1]];

visited[stack[depth+1]]=0;

}

else /* i!=N,表示for循环的终止是由于寻找到了当前城市的一个可行的邻近城市*/

{

//

printf("%d %d %d %d %d %d\n",cur,i,m1[stack[depth]][i],m2[stack[depth]][i],currentCost,currentDist);//getchar();

currentDist+=m1[stack[depth]][i]; /* 把从当前所处城市到所找到的可行城市的距离加入currentDist */

currentCost+=m2[stack[depth]][i]; /* 把从当前所处城市到所找到的可行城市的代价加入currentCost */

depth++; /* 所找到的可行城市进栈*/

stack[depth]=i; /* 更新栈顶指针，指向所找到的可行城市*/

stack[depth+1]=0;

visited[i]=1; /* 修改所找到的城市的访问标志*/

if(i==N-1) /* i==N-1表示访问到了乙城市，完成了所有城市的一次搜索，找到一条通路*/

{

// printf("here\n");

for(j=0;j<=depth;j++) /* 保存当前找到的通路所经过的所有节点*/

bestPath[j]=stack[j];

bestLength=depth; /* 保存当前找到的通路所经过的所有节点的节点数*/

shortestDist=currentDist; /* 保存当前找到的通路的距离之和*/

minimumCost=currentCost; /* 保存当前找到的通路的代价之和*/

//costBound=currentCost;

distBound=currentDist; /* 更新剪枝的路径边界，如果以后所找到的通路路径之和大于目前通路的路径之和，就剪枝*/

if(minimumCost>1500) continue; /* 如果当前找到的通路的代价之和大于1500，则放弃这条通路*/

printf("最短路径:%3d,路径代价:%3d,所经历的节点数目:%3d,所经历的节点如下:\n",shortestDist,minimumCost,bestLength+1); /* 输出找到的通路的结果*/

bestPath[bestLength]=49;

for(i=0;i<=bestLength;i++) /* 输出所找到的通路所经过的具体的节点*/

printf("%3d ",bestPath[i]+1);

printf("\n");

depth--; /* 连续弹出栈顶的两个值，进行回溯，开始寻找新的可行的通路*/

currentDist-=m1[stack[depth]][stack[depth+1]];

currentCost-=m2[stack[depth]][stack[depth+1]];

visited[stack[depth+1]]=0;

depth--;

currentDist-=m1[stack[depth]][stack[depth+1]];

currentCost-=m2[stack[depth]][stack[depth+1]];

visited[stack[depth+1]]=0;

// getchar();

}

}

}

}

北航2010-2011年研究生数值分析期末模拟试卷1-3

数值分析模拟试卷1 一、填空（共30分，每空3分） 1 设??? ? ??-=1511A ，则A 的谱半径=)(a ρ______，A 的条件数)(1A cond =________. 2 设 ,2,1,0,,53)(2==+=k kh x x x f k ，则],,[21++n n n x x x f ＝________, ],,[321+++n n n n x x x x f ，＝________. 3 设?????≤≤-++≤≤+=2 1,121 0,)(2 323x cx bx x x x x x S ，是以0，1，2为节点的三次样条函数，则b=________,c=________. 4 设∞=0)]([k k x q 是区间［0，1］上权函数为x x =)(ρ的最高项系数为1的正交多项式族，其中1)(0=x q ，则 ?=1 )(dx x xq k ________，=)(2 x q ________. 5 设???? ??????=11001a a a a A ，当∈a ________时，必有分解式，其中L 为下三角阵，当 其对角线元素)3,2,1(=i L ii 满足条件________时，这种分解是唯一的. 二、（14分）设4 9,1,41,)(2102 3 === =x x x x x f , （1）试求)(x f 在]4 9,41[上的三次Hermite 插值多项式)(x H 使满足 2,1,0),()(==i x f x H i i ，)()(11x f x H '='. （2）写出余项)()()(x H x f x R -=的表达式. 三、（14分）设有解方程0cos 2312=+-x x 的迭代公式为n n x x cos 3 2 41+ =+， （1） 证明R x ∈?0均有? ∞ →=x x n x lim （? x 为方程的根）； （2） 取40=x ，用此迭代法求方程根的近似值，误差不超过，列出各次迭代值； （3）此迭代的收敛阶是多少？证明你的结论. 四、(16分) 试确定常数A ，B ，C 和，使得数值积分公式 有尽可能高的代数精度. 试问所得的数值积分公式代数精度是多少？它是否为Gauss 型的？

北航数值分析大作业一

《数值分析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 #include #include #include #include #include #include #define E 1.0e-12 /*定义全局变量相对误差限*/ int max2(int a,int b) /*求两个整型数最大值的子程序*/ { if(a>b) return a; else return b; } int min2(int a,int b) /*求两个整型数最小值的子程序*/ { if(a>b) return b; else return a; } int max3(int a,int b,int c) /*求三整型数最大值的子程序*/ { int t; if(a>b) t=a; else t=b; if(t

北航数值分析大作业第二题

数值分析第二次大作业 史立峰 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。

北航数值分析报告大作业第八题

北京航空航天大学 数值分析大作业八 学院名称自动化 专业方向控制工程 学号 学生姓名许阳 教师孙玉泉 日期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**+==。

北航数值分析大作业第一题幂法与反幂法

《数值分析》计算实习题目 第一题： 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 #include void init_a();//初始化A double get_an_element(int,int);//取A 中的元素函数 double powermethod(double);//原点平移的幂法 double inversepowermethod(double);//原点平移的反幂法 int presolve(double);//三角LU 分解 int solve(double [],double []);//解方程组 int max(int,int); int min(int,int); double (*u)[502]=new double[502][502];//上三角U 数组 double (*l)[502]=new double[502][502];//单位下三角L 数组 double a[6][502];//矩阵A int main() { int i,k; double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;

北航数值分析计算实习报告一

航空航天大学 《数值分析》计算实习报告 第一大题 学院：自动化科学与电气工程学院 专业：控制科学与工程 学生姓名： 学号： 教师： 电话： 完成日期： 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 的主对角线的乘积。

北航数值分析报告第三次大作业

数值分析第三次大作业 一、算法的设计方案： （一）、总体方案设计： 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 =-，

北航2010-2015年研究生数值分析报告期末模拟试卷与真题

数值分析模拟卷A 一、填空（共30分，每空3分） 1 设???? ??-=1511A ，则A 的谱半径=)(a ρ______，A 的条件数)(1A cond =________. 2 设 ,2,1,0,,53)(2==+=k kh x x x f k ，则],,[21++n n n x x x f ＝________, ],,[321+++n n n n x x x x f ，＝________. 3 设?????≤≤-++≤≤+=2 1,1210,)(2323x cx bx x x x x x S ，是以0，1，2为节点的三次样条函数，则b=________,c=________. 4 设∞=0)]([k k x q 是区间［0，1］上权函数为x x =)(ρ的最高项系数为1的正交多项式族，其中1)(0=x q ，则?=10)(dx x xq k ________，=)(2x q ________. 5 设??????????=11001a a a a A ，当∈a ________时，必有分解式，其中L 为下三角阵，当 其对角线元素)3,2,1(=i L ii 满足条件________时，这种分解是唯一的. 二、（14分）设4 9,1,41,)(21023 ====x x x x x f , （1）试求)(x f 在]4 9 ,41[上的三次Hermite 插值多项式)(x H 使满足2,1,0),()(==i x f x H i i ，)()(11x f x H '='. （2）写出余项)()()(x H x f x R -=的表达式.

三、（14分）设有解方程0cos 2312=+-x x 的迭代公式为n n x x cos 3 241+ =+， （1） 证明R x ∈?0均有?∞→=x x n x lim （?x 为方程的根）； （2） 取40=x ，用此迭代法求方程根的近似值，误差不超过 ，列出各次迭代值； （3）此迭代的收敛阶是多少？证明你的结论. 四、(16分) 试确定常数A ，B ，C 和，使得数值积分公式 有尽可能高的代数精度. 试问所得的数值积分公式代数精度是多少？它是否为Gauss 型的？ 五、（15分） 设有常微分方程的初值问题???=='00 )(),(y x y y x f y ，试用Taylor 展开原理构造形如)()(11011--++++=n n n n n f f h y y y ββα的方法，使其具有二阶精度，并推导其局部截断误差主项.

北航数值分析第一次大作业(高斯gauss lu分解)

一、问题分析及算法描述 编写程序，分别用列主元的Gauss 消去法和LU 分解法求解下面线型代数方程组AX=b 的解，其中A 为N ×N 矩阵，N=50，其中第i(i ≥1)行、第j(i ≥1)列元素 a ij =1 i+j ?1， 右端向量b 的第i(i ≥1)个分量为 b i = 10 i+j ?1N j=1. 列主元素Gauss 消去过程中，要用到两种初等行变换。第一种，交换两行的位置；第二种，用一个数乘某一行加到另一行上。在第k 次消元之前，先对增广矩阵 A (k),b (k) 作第一种行变换，使得a ik (k) 中绝对值最大的元素交换到第k 行的主对角线位置上，然后再使用第二种行变换进行消元。如此往复，最后得到一个上三角系数矩阵，并回代求解解向量。由于每次消元前选取了列主元素，因此与顺序Guass 消元法相比，可提高数值计算的稳定性，且其计算量与顺序Guass 消元法相同。列主元的Gauss 消去法要求系数矩阵A 非奇异。 LU 分解法，即通过一系列初等行变换将系数矩阵A 分解成一个下三角矩阵L 与一个上三角矩阵U 的乘积，进一步通过求解两个三角矩阵得出解向量。若L 为单位下三角矩阵，U 是上三角矩阵，则称为Doolittle 分解；若L 为下三角矩阵，U 是单位上三角矩阵，则称为Crout 分解。若系数矩阵A 的前n-1阶顺序主子式不为零，则Doolittle\Crout 分解具有唯一性。若在每步行变换中选取主元，可提高数值计算稳定性。本算例中采用选主元的Doolittle 分解。 通过分析可知，本算例中待求解线型方程组系数矩阵为非奇异矩阵，且其前n-1阶顺序主子式不为零。方程组的解向量为x = 10,10,?,10 T 。满足列主元高斯消去法以及LU 分解法的基本使用条件。为了验证上述两种方法对本算例的适用性，笔者利用Microsoft Visual C++6.0编写了该算例的列主元高斯消去法以及LU 分解法的程序代码，并进行了运算求解。

北航数值分析大作业第二题精解

目标：使用带双步位移的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 #include #include #define N 10 #define E 1.0e-12 #define MAX 10000 //以下是符号函数 double sgn(double a) { double z; if(a>E) z=1; else z=-1; return z; } //以下是矩阵的拟三角分解 void nishangsanjiaodiv(double A[N][N]) { int i,j,k; int m=0; double d,c,h,t; double u[N],p[N],q[N],w[N]; for(i=0;i

北航数值分析课程第一次大作业讲解

《数值分析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 #include #include #include #include #include #include #define E 1.0e-12 /*定义全局变量相对误差限*/ int max2(int a,int b) /*求两个整型数最大值的子程序*/ { if(a>b) return a; else return b; } int min2(int a,int b) /*求两个整型数最小值的子程序*/ { if(a>b) return b; else return a; } int max3(int a,int b,int c) /*求三整型数最大值的子程序*/ { int t; if(a>b) t=a; else t=b; if(t

北航研究生数值分析A作业三

《数值分析（A）》计算实习题目三

一、题目 关于x，y，t，u，v，w的下列方程组 0.5cost+u+v+w-x=2.67 t+0.5sinu+v+w-y=1.07 0.5t+u+cosv+w-x=3.74 t+0.5u+v+sinw-y=0.79 以及关于z，t，u的下列二维数表 确定了一个二元函数z=f(x，y)。 1．试用数值方法求出f(x，y)在区域D={(x,y)︱0≤x≤0.8,0.5≤y≤1.5}上的一个

近似表达式 ,0 (,)k r s rs r s p x y c x y == ∑ 要求p(x,y)一最小的k 值达到以下的精度 10 20 27 ((,)(,))10i j i j i j f x y p x y σ-=== -≤∑∑ 其中x i =0.08i ，y j =0.5+0.05j 。 2．计算f(x i *,y j *),p(x i *,y j *)(i=1,2,…,8;j=1,2,…,5)的值，以观察p(x,y)逼近f(x,y) 的效果，其中x i *=0.1i,y j * =0.5+0.2j 。

二、算法设计方案 1.将0.08(0,1,,10) i x i i * == 和0.50.05(0,1,,20)j y j j *=+= 代入非线性方程 组中，用牛顿法解出i t 和j u ； 2.以采取分片二次插值，选择（m ，n ）满足 ,232 2 m i m h h t t t m -<≤+ ≤≤ ,23 2 2 n j n u u u n τ τ - <≤+ ≤≤ 如果12 i h t t ≤+ 或42 i h t t >- ，则m=1或4；如果12 j u u τ ≤+ 或42 j u u τ >- ，则n=1 或n=4。选择(,)(1,,1;1,,1) k r t u k m m m r n n n =-+=-+为插值节点，相应的Lagrange 形式的插值多项式为 ),()(~ )(),(1 11 1 22r k r m m k n n r k u t f u l t l u t p ∑∑ +-=+-== 其中 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 r w n w w r w r y y y y u l (r=n-1, n, n+1) 并将i t 和j u 代入22(,)p t u ，便得到了数表,, (,) i j i j x y f x y 。 3.进行曲面拟和系数矩阵[]rs c * =C ，1 1 () () T T T --=C B B B U G G G 其中 001110 1011[()]1k k r i k x x x x x x x ??? ?? ??==?? ????? ?B ，001110 1011[()]1 k k s j k y y y y G y y y ψ?? ????== ?????? ?? [(,)]i j f x y =U k 从0逐渐增大，直到7 10 σ -≤，便得到了要求精度的系数rs c 。 4.由前面得到的函数关系，根据重新取值的x ，y 可以分别得到新的数表，比 较两组数据观察逼近效果

BUAA数值分析大作业三

北京航空航天大学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. 指出有效数49×102，0.0490，490.00的绝对误差限、相对误差限和有效数字位数． 2. 将 3.142作为π的近似值，它有几位有效数字，相对误差限和绝对误差限各为多少？ 3. 要使101的近似值x * 的相对误差限不超过4102 1?×，问查开方表时x * 需要保留几位有效数字？ 4. 已知近似数x * 有两位有效数字，试估计其相对误差限． 5. 设x * 为x 的近似数, 证明n x * 的相对误差大约为x * 相对误差的n 1倍． 6. 某矩形的长和宽大约为100cm 和50cm, 应该选用最小刻度为多少cm 的测量工具, 才能保证计算出的面积误差（绝对值）不超过0.15cm 2． 7. 已知三角形面积c ab S sin 2 1=，测量a , b , c 时产生的相对误差为)(*a e r ，)(*b e r ，)(*c e r ，其中2 ,0*π<>2时的情形．用所设计的算法以及二次方程求根公式计算05.240=p ，00.1=q 时方程根的近似值(计算过程保留2位小数)，并给出它们的相对误差限（根的准确值为L 0916683.4801?=x , L 002082935.02?=x ）．

北航数值分析大作业第二次

《数值分析》计算实习作业 （第二题）

算法设计方案： 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；最后将所有的实特征值的特征向量输出。

北航数值分析计算实习报告一

北航数值分析计算实习报告一

————————————————————————————————作者：————————————————————————————————日期:

北京航空航天大学 《数值分析》计算实习报告 第一大题 学院：自动化科学与电气工程学院 专业: 控制科学与工程 学生姓名： 学号: 教师: 电话: 完成日期: 2015年11月6日 北京航空航天大学 Ｂeｉjing Univeｒsiｔy of Ａeｒｏnaｕtｉcs and Ａstroｎautｉcs

实习题目： 第一题 设有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.求Ａ的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,,2,1(???=k k i λ。 ３.求Ａ的(谱范数)条件数2)A (cond 和行列式d ｅtA 。 说明： 1.在所用的算法中,凡是要给出精度水平ε的，都取12-10=ε。 2.选择算法时,应使矩阵A 的所有零元素都不储存。 3.打印以下内容： （1）全部源程序； （2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。 4.采用e 型输出实型数，并且至少显示１２位有效数字。

北航数值分析大作业3

一、算法设计方案 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 ξ???? ? ?=?????? ，

北航数值分析计算实习报告一

北航数值分析计算实习 报告一 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 μ。