数值分析计算实习一

数值分析大作业

第一题

学号：

姓名：

学院：

授课教师：

2016.10.24

一、算法设计方案

1、求λ1，λ501和λs

a)利用幂法计算出矩阵A 按模最大的特征值，设其为λm 。 b)令矩阵B =A ?λm I （I 为单位矩阵），同样利用幂法计算出矩阵B 按模最大的特征值λm ′。

c)令λm ′′=λm ′+λm 。由计算过程可知λm 和λm ′′

分别为矩阵A 所有特征值按大小排序后，其中两端的值。即，λ1=min?{λm ,λm ′′}，λ501=max?{λm ,λm ′′}。

d) 利用反幂法计算λs 。其中通过使用doolittle 三角分解法求解线性方程组AX =I 的形式求得矩阵A 的逆A ?1。进行doolittle 分解，可以节省内存空间而把矩阵A 压缩，将A 中带上的元素存储为数组C [5][501]。

2、求A 的与数μk =λ1+k

λ501?λ1

最接近的特征值λi k (k =1,2, (39)

e) 令矩阵D k =A ?μk I ，利用反幂法计算出矩阵D k 按模最小的特征值λi k ′,则λi k =λi k ′

μk 。其中，同样使用doolittle 三角分解法求解线性方程组D k X =I 的形式求得矩阵D k 的逆D k ?1,其中k =1,2, (39)

3）求矩阵A 的条件数cond(A )2和行列式det(A )

f) 计算矩阵的条件数cond(A)2=|λ

m λs

|。

g)计算矩阵的行列式det(A )为A 三角分解后得到的上三角矩阵U 的对角线元素之积。

二、程序

#include #include #include using namespace std;

double LU(double C[6][502]);//LU 分解

double solve(double CC[6][502],double array[502],double U[502]);//解线性方程组 double mefa(double arr[][502]);//幂法

double fanmefa(double array[][502]);//反幂法 int max1(int a,int b,int c); int max2(int a ,int b); int min1(int a,int b);

int main() { double a[502],d[6][502],d1[6][502],lam1,lam501,temp,lams,lamki,uk[40],condA,detA=1; //d 为压缩矩阵,d1为进行原点平移时求lam1和lam501的压缩矩阵 int i,j,k; for(i=1;i<=501;i++) { a[i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);

//求解主对角线上的元素

}

for (i=1;i<=5;i++)

for (j=1;j<=501;j++)

{

if(i==1)

{

if ((j==1)||(j==2))d[i][j]=0;

else d[i][j]=-0.064;

}

else if (i==2)

{

if (j==1)d[i][j]=0;

else d[i][j]=0.16;

}

else if (i==3)

d[i][j]=a[j];

else if(i==4)

{

if (j==501)d[i][j]=0;

else d[i][j]=0.16;

}

else if(i==5)

{

if ((j==500)||(j==501))d[i][j]=0;

else d[i][j]=-0.064;

}

//The purpose of this part is to assign values to matrix d[][],and d[][] is a condensed matrix for A.

lam1=mefa(d);

for (i=1;i<=5;i++)

for (j=1;j<=501;j++)

{

if(i==3)d1[3][j]=d[3][j]-lam1;

else d1[i][j]=d[i][j];

}

//设置d1的值,进行原点平移

lam501=mefa(d1);

lam501=lam501+lam1;

if(lam1>lam501)

{

temp=lam1;

lam1=lam501;

lam501=temp;

}

//将lam1和lam501中小的值赋给lam1

for(i=1;i<=39;i++)

uk[i]=lam1+i*(lam501-lam1)/40;

//求解uk的值

cout<<"近似特征值的结果为:"<

for(i=1;i<=39;i++)

{

for(k=1;k<=5;k++)

{for(j=1;j<=501;j++)

{

if(k==3)d1[k][j]=d[k][j]-uk[i];

else d1[k][j]=d[k][j];

}

//求解uk时的原点平移

LU(d1);

lamki=fanmefa(d1);

cout<<"λk"<

cout<

if(i%2==0)cout<

else cout<<" ";

}

cout<

LU(d); //进行LU分解

lams=fanmefa(d); //进行反幂法

condA=fabs(lam1/lams);//求矩阵的条件数

for(i=1;i<=501;i++)

detA=detA*d[3][i];//求行列式的值

cout<<"矩阵的最小特征值为：λ1="<

cout<<"矩阵的最大特征值为：λ501="<

cout<<"按模最小的特征值为：λs="<

cout<<"矩阵的条件数为：cond(A)="<

cout<<"行列式的值为：det(A)="<

return(0);

}

double mefa(double arr[][502])

//幂法求解按模最大的特征值

{

double u[502],y[502],eita,beta=0,firstbeta=0,sum,sum1;

int i,j,k;

for(i=1;i<=501;i++)

u[i]=1;

for (k=1;k<10000;k++)

{

sum=0;

for(j=1;j<=501;j++)

sum=sum+u[j]*u[j];

eita=sqrt(sum);

for(j=1;j<=501;j++)

y[j]=u[j]/eita;

for(j=1;j<=501;j++)

{

if(j==1)u[j]=arr[3][1]*y[1]+arr[2][2]*y[2]+arr[1][3]*y[3];

if(j==2)u[j]=arr[4][1]*y[1]+arr[3][2]*y[2]+arr[2][3]*y[3]+arr[1][4]*y[4];

if((j>=3)&&(j<=499))

u[j]=arr[5][j-1]*y[j-2]+arr[4][j-1]*y[j-

1]+arr[3][j]*y[j]+arr[2][j+1]*y[j+1]+arr[1][j+2]*y[j+2];

if(j==500)u[j]=arr[5][498]*y[498]+arr[4][499]*y[499]+arr[3][500]*y[500]+arr[2][501] *y[501];

if(j==501)u[j]=arr[5][499]*y[499]+arr[4][500]*y[500]+arr[3][501]*y[501];

}

//幂法中利用矩阵相乘求解uk的值

firstbeta=beta;

sum1=0;

for(j=1;j<=501;j++)

sum1=sum1+y[j]*u[j];

beta=sum1;

if(fabs((beta-firstbeta)/beta)<=1e-12) break;

}

return(beta);

}

double fanmefa(double array[][502])

//利用反幂法求解矩阵按模最小的特征值

{

double u[502],y[502],eita,firstbeta=0,beta=0,sum,sum1;

int i,j,k;

for(i=1;i<=501;i++)

u[i]=1;

for (k=1;k<10000;k++)

{

sum=0;

for(j=1;j<=501;j++)

sum=sum+u[j]*u[j];

eita=sqrt(sum);

for(j=1;j<=501;j++)

y[j]=u[j]/eita;

solve(array,y, u); //调用解线性方程组的程序解出uk

firstbeta=beta;

sum1=0;

for(j=1;j<=501;j++)

sum1=sum1+y[j]*u[j];

beta=sum1;

if(fabs((firstbeta-beta)/firstbeta)<=1e-12) break;//

}

return(1/beta);

}

double LU(double C[6][502])

//此函数对矩阵进行LU分解

{

double sum1,sum2;

int i,j,k,t;

for (k=1;k<=501;k++)

{

for (j=k;j<=min1(k+2,501);j++)

{

sum1=0;

for (t=max1(1,k-2,j-2);t<=k-1;t++)

sum1=sum1+C[k-t+3][t]*C[t-j+3][j];

C[k-j+3][j]=C[k-j+3][j]-sum1;

}

for(i=k+1;i<=min1(k+2,501);i++)

{

sum2=0;

for(t=max1(1,i-2,k-2);t<=k-1;t++)

sum2=sum2+C[i-t+3][t]*C[t-k+3][k];

C[i-k+3][k]=(C[i-k+3][k]-sum2)/C[3][k];

}

return(0);

}

double solve(double CC[6][502],double array[502],double U[502]) //此函数的功能是解线性方程组

{

double sum1,sum2,b[502];

int i,t;

//b[1]=array[1];

for(i=1;i<=501;i++)

b[i]=array[i];

for(i=2;i<=501;i++)

{

sum1=0;

for(t=max2(1,i-2);t<=i-1;t++)

sum1=sum1+CC[i-t+3][t]*b[t];

b[i]=b[i]-sum1;

}

U[501]=b[501]/CC[3][501];

for(i=500;i>=1;i--)

{

sum2=0;

for(t=i+1;t<=min1(i+2,501);t++)

sum2=sum2+CC[i-t+3][t]*U[t];

U[i]=(b[i]-sum2)/CC[3][i];

}

return(0);

}

int max1(int a,int b,int c)

{

int t=a;

if(b>t)t=b;

if(c>t)t=c;

return(t);

}

int max2(int a ,int b)

{

int t=a;

if(b>t)t=b;

return(t);

}

int min1(int a,int b) {

int t=a;

if(b

return(t);

}

三、计算结果

计算结果如下所示：

1、矩阵的最小特征值为：λ1=-1.070011361502e+001

矩阵的最大特征值为：λ501=9.724634098777e+000

按模最小的特征值为：λs=-5.557910794230e-003

2、近似特征值的结果为:

λk

1=-1.018293403315e+001 λk

=-9.585707425068e+000

λk

3=-9.172672423928e+000 λk

=-8.652284007898e+000

λk

5=-8.0934********e+000 λk

=-7.659405407692e+000

λk

7=-7.119684648691e+000 λk

=-6.611764339397e+000

λk

9=-6.0661********e+000 λk

=-5.585101052628e+000

λk

11=-5.114083529812e+000 λk

=-4.578872176865e+000

λk

13=-4.096470926260e+000 λk

=-3.554211215751e+000

λk

15=-3.0410********e+000 λk

=-2.533970311130e+000

λk

17=-2.003230769563e+000 λk

=-1.503557611227e+000

λk

19=-9.935586060075e-001 λk

=-4.870426738850e-001

λk

21=2.231736249575e-002 λk

=5.324174742069e-001

λk

23=1.052898962693e+000 λk

=1.589445881881e+000

λk

25=2.060330460274e+000 λk

=2.558075597073e+000

λk

27=3.080240509307e+000 λk

=3.613620867692e+000

λk

29=4.0913********e+000 λk

=4.603035378279e+000

λk

31=5.132924283898e+000 λk

=5.594906348083e+000

λk

33=6.080933857027e+000 λk

=6.680354092112e+000

λk

35=7.293877448127e+000 λk

=7.717111714236e+000

λk

37=8.225220014050e+000 λk

=8.648666065193e+000

λk

=9.254200344575e+000

3、矩阵的条件数为：cond(A)=1.925204273902e+003

行列式的值为：det(A)=2.772786141752e+118

四、讨论：初始向量的选取对计算结果的影响

当选取初始向量为

(1,1,1,1,...,1),1,2,3,...,501

u i

时，得到以上结果。

现修改初始向量为

[1]1,[]0,2,3,...,501

u u i i

===，得到结果如下：

可以看出，此结果与近似值相差较多。令初始向量为

1,0,(1,2,...,100;101,102, (501)

i j u u i j ====，得出的结果如下所示，

可以看出，此结果也与正确结果相差较多，但是更加靠近正确结果一些。对算法进行分析，如果选取的向量

u 使得

10α=，经管存在舍入误差使其不完全等于

0，但可能由于计算精度较高，使得无法求出按模最大的特征值，而求出次大的特征值。因

此在使用幂法和反幂法时，最好多代入几组初始向量，以确认得到的是正确结果。

数值分析大作业-三、四、五、六、七

大作业三 1. 给定初值 0x 及容许误差，编制牛顿法解方程f (x )=0的通用程序. 解：Matlab 程序如下：函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0; end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +??=-= ?-???解：Matlab 程序为：（1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; while nerrorlim n=n+1; else break ; end x0=x; end disp(['迭代次数: n=',num2str(n)]) disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)]) （2）子函数非线性函数f function y=f(x) y=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918); end

数值计算方法试题及答案

数值计算方法试题一一、填空题（每空1分，共17分） 1、如果用二分法求方程在区间内的根精确到三位小数，需对分（）次。 2、迭代格式局部收敛的充分条件是取值在（）。 3、已知是三次样条函数，则 =( )，=（），=（）。 4、是以整数点为节点的Lagrange插值基函数，则 ( )，( )，当时( )。 5、设和节点则和。 6、5个节点的牛顿-柯特斯求积公式的代数精度为，5个节点的求积公式最高代数精度为。 7、是区间上权函数的最高项系数为1的正交多项式族，其中，则。 8、给定方程组，为实数，当满足，且时，SOR迭代法收敛。 9、解初值问题的改进欧拉法是阶方法。 10、设，当（）时，必有分解式，其中为下三角阵，当其对角线元素满足（）条件时，这种分解是唯一的。二、二、选择题（每题2分） 1、解方程组的简单迭代格式收敛的充要条件是（）。（1）, (2) , (3) , (4) 2、在牛顿-柯特斯求积公式：中，当系数是负值时，公式的稳定性不能保证，所以实际应用中，当（）时的牛顿-柯特斯求积公式不使用。（1），（2），（3），（4），（1）二次；（2）三次；（3）四次；（4）五次 4、若用二阶中点公式求解初值问题，试问为保证该公式绝对稳定，步长的取值范围为（）。 (1), (2), (3), (4)

三、1、 2、（15 (1)(1) 试用余项估计其误差。（2）用的复化梯形公式（或复化 Simpson公式）计算出该积分的近似值。四、1、（15分）方程在附近有根，把方程写成三种不同的等价形式（1）对应迭代格式；(2)对应迭代格式；（3）对应迭代格式。判断迭代格式在的收敛性，选一种收敛格式计算附近的根，精确到小数点后第三位。选一种迭代格式建立Steffensen迭代法，并进行计算与前一种结果比较，说明是否有加速效果。 2、（8分）已知方程组，其中，（1）（1）列出Jacobi迭代法和Gauss-Seidel迭代法的分量形式。（2）（2）求出Jacobi迭代矩阵的谱半径，写出SOR 迭代法。五、1、（15分）取步长，求解初值问题用改进的欧拉法求的值；用经典的四阶龙格—库塔法求的值。 2、（8分）求一次数不高于4次的多项式使它满足 ,,,, 六、（下列2题任选一题，4分） 1、1、数值积分公式形如（1）（1）试确定参数使公式代数精度尽量高；（2）设，推导余项公式，并估计误差。 2、2、用二步法求解常微分方程的初值问题时，如何选择参数使方法阶数尽可能高，并求局部截断误差主项，此时该方法是几阶的。数值计算方法试题二一、判断题：（共16分，每小题２分）１、若是阶非奇异阵，则必存在单位下三角阵和上三角阵，使唯一成立。（）

数值分析大作业三四五六七

flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

while flag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1; while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0;

end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +?? =-= ?-???解： Matlab 程序为：（1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1;

数值分析(计算方法)总结

第一章绪论误差来源：模型误差、观测误差、截断误差（方法误差）、舍入误差是的绝对误差，是的误差，为的绝对误差限（或误差限）为的相对误差，当较小时，令相对误差绝对值得上限称为相对误差限记为：即：绝对误差有量纲，而相对误差无量纲若近似值的绝对误差限为某一位上的半个单位，且该位直到的第一位非零数字共有n位，则称近似值有n位有效数字，或说精确到该位。例：设x==3.1415926…那么,则有效数字为1位，即个位上的3，或说精确到个位。科学计数法：记有n位有效数字，精确到。由有效数字求相对误差限：设近似值有n位有效数字，则其相对误差限为由相对误差限求有效数字：设近似值的相对误差限为为则它有n位有效数字令 1.x+y近似值为和的误差（限）等于误差（限）的和 2.x-y近似值为 3.xy近似值为 4. 1．避免两相近数相减 2．避免用绝对值很小的数作除数 3．避免大数吃小数

4．尽量减少计算工作量第二章非线性方程求根 1.逐步搜索法设f (a) <0, f (b)> 0，有根区间为 (a, b)，从x0=a出发，按某个预定步长(例如h=(b-a)/N)一步一步向右跨，每跨一步进行一次根的搜索，即判别f(x k)=f(a+kh)的符号，若f(x k)>0(而f(x k-1)<0),则有根区间缩小为[x k-1,x k] (若f(x k)=0，x k即为所求根), 然后从 x k-1出发，把搜索步长再缩小，重复上面步骤，直到满足精度：|x k-x k-1|< 为止，此时取 x*≈(x k+x k-1)/2作为近似根。 2.二分法设f(x)的有根区间为[a,b]= [a0,b0], f(a)<0, f(b)>0.将[a0,b0]对分，中点x0= ((a0+b0)/2),计算f(x0)。 3.比例法一般地，设 [a k,b k]为有根区间，过(a k, f(a k))、 (b k, f(b k))作直线，与x轴交于一点x k,则： 1.试位法每次迭代比二分法多算一次乘法，而且不保证收敛。 2.比例法不是通过使求根区间缩小到0来求根，而是在一定条件下直接构造出一个点列（递推公式），使该点列收敛到方程的根。——这正是迭代法的基本思想。事先估计: 事后估计局部收敛性判定定理：局部收敛性定理对迭代函数的要求较弱，但对初始点要求较高，即初始点必须选在精确解的附近 Steffensen迭代格式： Newton法： Newton下山法：是下山因子弦割法：

数值计算方法大作业

目录第一章非线性方程求根 (3) 1.1迭代法 (3) 1.2牛顿法 (4) 1.3弦截法 (5) 1.4二分法 (6) 第二章插值 (7) 2.1线性插值 (7) 2.2二次插值 (8) 2.3拉格朗日插值 (9) 2.4分段线性插值 (10) 2.5分段二次插值 (11) 第三章数值积分 (13) 3.1复化矩形积分法 (13) 3.2复化梯形积分法 (14) 3.3辛普森积分法 (15) 3.4变步长梯形积分法 (16) 第四章线性方程组数值法 (17) 4.1约当消去法 (17) 4.2高斯消去法 (18) 4.3三角分解法 (20)

4.4雅可比迭代法 (21) 4.5高斯—赛德尔迭代法 (23) 第五章常积分方程数值法 (25) 5.1显示欧拉公式法 (25) 5.2欧拉公式预测校正法 (26) 5.3改进欧拉公式法 (27) 5.4四阶龙格—库塔法 (28)

数值计算方法第一章非线性方程求根 1.1迭代法程序代码： Private Sub Command1_Click() x0 = Val(InputBox("请输入初始值x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = (Exp(2 * x0) - x0) / 5 If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求f(x)=e2x-6x=0在x=0.5附近的根（ep=10-10）

1.2牛顿法程序代码： Private Sub Command1_Click() b = Val(InputBox("请输入被开方数x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = x0 - (x0 ^ 2 - b) / (2 * b) If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求56的值。（ep=10-10）

北航数值分析大作业一

《数值分析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

数值分析大作业

数值分析报大作业班级：铁道2班专业：道路与铁道工程姓名：蔡敦锦学号：13011260

一、序言该数值分析大作业是通过C语言程序编程在Microsoft Visual C++ 6.0编程软件上运行实现的。本来是打算用Matlab软间来计算非线性方程的根的。学习Matlab也差不多有一个多月了，感觉自己编程做题应该没什么问题了；但是当自己真心的去编程、运行时才发现有很多错误，花了一天时间修改、调试程序都没能得到自己满意的结果。所以，我选择了自己比较熟悉的C程序语言来编程解决非线性的求值问题，由于本作业是为了比较几种方法求值问题的收敛速度和精度的差异，选择了一个相对常见的非线性函数来反映其差异，程序运行所得结果我个人比较满意。编写C语言，感觉比较上手，程序出现问题也能比较熟练的解决。最终就决定上交一份C程序语言编程的求值程序了！

二、选题本作业的目的是为了加深对非线性方程求根方法的二分法、简单迭代法、、牛顿迭代法弦截法等的构造过程的理解；能将各种方法的算法描述正确并且能够改编为程序并在计算机上实现程序的正确合理的运行，能得到自己满意的结果，并且能调试修改程序中可能出现的问题和程序功能的增减修改。本次程序是为了比较各种方法在求解同一非线性方程根时，在收敛情况上的差异。为了达到上面的条件我选择自己比较熟悉的语言—C语言来编程，所选题目为计算方程f(x)=x3-2x-5=0在区间[2,3]内其最后两近似值的差的绝对值小于等于5 ?的根的几种方法的比较。 110- 本文将二分法、牛顿法、简单迭代法、弦截法及加速收敛法这五种方法在同一个程序中以函数调用的方式来实现，比较简洁明了，所得结果能很好的比较，便于分析；发现问题和得出结论。

数值分析计算方法试题集及答案

数值分析复习试题第一章绪论一. 填空题 1.* x 为精确值 x 的近似值；() **x f y =为一元函数 ()x f y =1的近似值； ()**,*y x f y =为二元函数()y x f y ,2=的近似值，请写出下面的公式：**e x x =-： *** r x x e x -= ()()()*'1**y f x x εε≈? ()() () ()'***1**r r x f x y x f x εε≈ ? ()()()() ()* *,**,*2**f x y f x y y x y x y εεε??≈?+??? ()()()()() ** * *,***,**222r f x y e x f x y e y y x y y y ε??≈ ?+??? 2、计算方法实际计算时，对数据只能取有限位表示，这时所产生的误差叫舍入误差。 3、分别用2.718281，2.718282作数e 的近似值，则其有效数字分别有 6 位和 7 位；又取 1.73≈-21 1.73 10 2 ≤?。 4、设121.216, 3.654x x ==均具有3位有效数字，则12x x 的相对误差限为 0.0055 。 5、设121.216, 3.654x x ==均具有3位有效数字，则12x x +的误差限为 0.01 。 6、已知近似值 2.4560A x =是由真值T x 经四舍五入得到,则相对误差限为 0.0000204 . 7、递推公式,??? ? ?0n n-1y =y =10y -1,n =1,2, 如果取0 1.41y ≈作计算,则计算到10y 时,误差为 81 10 2 ?;这个计算公式数值稳定不稳定不稳定 . 8、精确值 14159265.3* =π，则近似值141.3*1=π和1415.3*2=π分别有 3

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

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

《数值计算方法》试题及答案

数值计算方法考试试题一、选择题（每小题4分，共20分） 1. 误差根据来源可以分为四类，分别是（ A ） A. 模型误差、观测误差、方法误差、舍入误差； B. 模型误差、测量误差、方法误差、截断误差； C. 模型误差、实验误差、方法误差、截断误差； D. 模型误差、建模误差、截断误差、舍入误差。 2. 若132)(3 56++-=x x x x f ，则其六阶差商 =]3,,3,3,3[6210 f （ C ） A. 0； B. 1； C. 2； D. 3 。 3. 数值求积公式中的Simpson 公式的代数精度为（ D ） A. 0； B. 1； C. 2； D. 3 。 4. 若线性方程组Ax = b 的系数矩阵A 为严格对角占优矩阵，则解方程组的Jacobi 迭代法和Gauss-Seidel 迭代法（ B ） A. 都发散； B. 都收敛 C. Jacobi 迭代法收敛，Gauss-Seidel 迭代法发散； D. Jacobi 迭代法发散，Gauss-Seidel 迭代法收敛。 5. 对于试验方程y y λ='，Euler 方法的绝对稳定区间为（ C ） A. 02≤≤-h ； B. 0785.2≤≤-h ； C. 02≤≤-h λ； D. 0785.2≤≤-h λ ；二、填空题（每空3分，共18分） 1. 已知 ? ??? ??--='-=4321,)2,1(A x ，则 =2 x 5，= 1Ax 16 ，=2A 22115+ 2. 已知 3)9(,2)4(==f f ，则 f (x )的线性插值多项式为)6(2.0)(1+=x x L ,且用线性插值可得f (7)= 2.6 。 3. 要使 20的近似值的相对误差界小于0.1%，应至少取 4 位有效数字。三、利用下面数据表， 1. 用复化梯形公式计算积分 dx x f I )(6 .28 .1? =的近似值；解：1.用复化梯形公式计算取 2.048 .16.2,4=-= =h n 1分分分分7058337 .55))6.2()2.08.1(2)8.1((22.04)) ()(2)((231 1 1 4=+++=++=∑∑=-=f k f f b f x f a f h T k n k k 10.46675 8.03014 6.04241 4.42569 3.12014 f (x ) 2.6 2.4 2.2 2.0 1.8 x

计算流体力学常用数值方法简介[1]

计算流体力学常用数值方法简介李志印　熊小辉　吴家鸣 (华南理工大学交通学院) 关键词　计算流体力学　数值计算一　前　言任何流体运动的动力学特征都是由质量守恒、动量守恒和能量守恒定律所确定的,这些基本定律可以由流体流动的控制方程组来描述。利用数值方法通过计算机求解描述流体运动的控制方程,揭示流体运动的物理规律,研究流体运动的时一空物理特征,这样的学科称为计算流体力学。计算流体力学是一门由多领域交叉而形成的一门应用基础学科,它涉及流体力学理论、计算机技术、偏微分方程的数学理论、数值方法等学科。一般认为计算流体力学是从20世纪60年代中后期逐步发展起来的,大致经历了四个发展阶段:无粘性线性、无粘性非线性、雷诺平均的N-S方程以及完全的N-S方程。随着计算机技术、网络技术、计算方法和后处理技术的迅速发展,利用计算流体力学解决流动问题的能力越来越高,现在许多复杂的流动问题可以通过数值计算手段进行分析并给出相应的结果。经过40年来的发展,计算流体力学己经成为一种有力的数值实验与设计手段,在许多工业领域如航天航空、汽车、船舶等部门解决了大量的工程设计实际问题,其中在航天航空领域所取得的成绩尤为显著。现在人们已经可以利用计算流体力学方法来设计飞机的外形,确定其气动载荷,从而有效地提高了设计效率,减少了风洞试验次数,大大地降低了设计成本。此外,计算流体力学也己经大量应用于大气、生态环境、车辆工程、船舶工程、传热以及工业中的化学反应等各个领域,显示了计算流体力学强大的生命力。随着计算机技术的发展和所需要解决的工程问题的复杂性的增加,计算流体力学也己经发展成为以数值手段求解流体力学物理模型、分析其流动机理为主线,包括计算机技术、计算方法、网格技术和可视化后处理技术等多种技术的综合体。目前计算流体力学主要向二个方向发展:一方面是研究流动非定常稳定性以及湍流流动机理,开展高精度、高分辩率的计算方法和并行算法等的流动机理与算法研究;另一方面是将计算流体力学直接应用于模拟各种实际流动,解决工业生产中的各种问题。二　计算流体力学常用数值方法流体力学数值方法有很多种,其数学原理各不相同,但有二点是所有方法都具备的,即离散化和代数化。总的来说其基本思想是:将原来连续的求解区域划分成网格或单元子区

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

目标：使用带双步位移的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

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

《数值分析》计算实习题目第一题： 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;

上海大学_王培康_数值分析大作业

数值分析大作业（2013年5月）金洋洋（12721512），机自系 1.下列各数都是经过四舍五入得到的近似值，试分别指出它们的绝对误差限, 相对误差限和有效数字的位数。 X1 =5.420, x 2 =0.5420, x 3=0.00542, x 4 =6000, x 5=50.610? 解：根据定义：如果*x 的绝对误差限不超过x 的某个数位的半个单位，则从*x 的首位非零数字到该位都是有效数字。显然根据四舍五入原则得到的近视值，全部都是有效数字。因而在这里有：n1=4, n2=4, n3=3, n4=4, n5=1 （n 表示x 有效数字的位数）对x1：有a1=5, m1=1 （其中a1表示x 的首位非零数字，m1表示x1的整数位数）所以有绝对误差限 143 11 (1)101022 x ε--≤ ?=? 相对误差限 31() 0.510(1)0.00923%5.4201 r x x x εε-?= == 对x2：有a2=5, m2=0 所以有绝对误差限 044 11 (2)101022 x ε--≤ ?=? 相对误差限 42() 0.510(2)0.00923%0.54202 r x x x εε-?= == 对x3：有a3=5, m3=-2 所以有绝对误差限 235 11 (3)101022 x ε---≤ ?=? 相对误差限 53() 0.510(3)0.0923%0.005423 r x x x εε-?= == 对x4：有a4=0, m4=4 所以有绝对误差限 4411(4)1022 x ε-≤?= 相对误差限 4() 0.5 (4)0.0083%6000 4 r x x x εε= = = 对x5：有a5=6, m5=5 所以有绝对误差限 514 11(5)101022 x ε-≤ ?=? 相对误差限 45() 0.510(5)8.3%600005 r x x x εε?= ==

数值计算方法试题及答案

数值计算方法试题一一、填空题（每空1分，共17分） 1、如果用二分法求方程043=-+x x 在区间]2,1[内的根精确到三位小数，需对分（）次。 2、迭代格式 ) 2(2 1-+=+k k k x x x α局部收敛的充分条件是α取值在（）。 3、已知?????≤≤+-+-+-≤≤=31)1()1()1(2110)(2 33x c x b x a x x x x S 是三次样条函数，则 a =( )， b =（）， c =（）。 4、)(,),(),(10x l x l x l n Λ是以整数点n x x x ,,,10Λ为节点的Lagrange 插值基函数，则 ∑== n k k x l 0)(( )， ∑== n k k j k x l x 0 )(( )，当2≥n 时 = ++∑=)()3(20 4 x l x x k k n k k ( )。 5、设 1326)(247+++=x x x x f 和节点,,2,1,0,2/Λ==k k x k 则=],,,[10n x x x f Λ 和=?07f 。 6、5个节点的牛顿-柯特斯求积公式的代数精度为，5个节点的求积公式最高代数精度为。 7、{}∞=0)(k k x ?是区间]1,0[上权函数x x =)(ρ的最高项系数为1的正交多项式族，其中1)(0=x ?，则?= 1 4)(dx x x ? 。

8、给定方程组?? ?=+-=-221121b x ax b ax x ，a 为实数，当a 满足，且 20<<ω时，SOR 迭代法收敛。 9、解初值问题00 (,)()y f x y y x y '=?? =?的改进欧拉法 ??? ??++=+=++++)],(),([2),(] 0[111] 0[1n n n n n n n n n n y x f y x f h y y y x hf y y 是阶方法。 10、设 ?? ??? ?????=11001a a a a A ，当∈a （）时，必有分解式T LL A =，其中L 为下三角阵，当其对角线元素)3,2,1(=i l ii 满足（）条件时，这种分解是唯一的。二、二、选择题（每题2分） 1、解方程组b Ax =的简单迭代格式g Bx x k k +=+)() 1(收敛的充要条件是（）。（1）1)(A ρ, (4) 1)(>B ρ 2、在牛顿-柯特斯求积公式： ?∑=-≈b a n i i n i x f C a b dx x f 0 )() ()()(中，当系数 ) (n i C 是负值时，公式的稳定性不能保证，所以实际应用中，当（）时的牛顿-柯特斯求积公式不使用。（1）8≥n ，（2）7≥n ，（3）10≥n ，（4）6≥n ， 3、有下列数表