数值分析

数值分析上机报告

前言

随着计算机技术的高速发展，越来越多的科技工作者使用计算机进行科学研究和解决工程技术问题。数值分析（或计算方法）课程的内容是科学工程计算的必备知识，已经成为众多理工科大学生、研究生的必修课程，越来越受到重视。

由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助很多编程软件来解决，得到一个满足误差限的解。本文所计算题目，均采用C++编程。在本文中使用C++编写了牛顿法、牛顿-Steffensen法方程求解的程序和雅格比法、高斯－赛德尔迭代法求解方程组的程序及Ru n ge-Kutt a4阶算法，并通过实例求解验证了其可行性，比较了求解同一种问题时不同方法之间的优缺性，其中包含解的精确度和解的收敛速度两个重要指标。

一 牛顿法和牛顿-Steffensen 法迭代求解的比较

1. 计算题目

分别用牛顿法，及基于牛顿算法下的Steffensen 加速法

(1) 求ln(x +sin x )=0的根。初值x0分别取0.1, 1,1.5, 2, 4进行计算。 (2) 求sin x =0的根。初值x0分别取1，1.4，1.6, 1.8，3进行计算。 分析其中遇到的现象与问题。

2. 计算过程和结果

1.对方程ln(x +sin x )=0，其导数有些复杂，我们可以对其进行变形，即求解x+sinx=1的解。使用牛顿法，令1sin )(-+=x x x f ，则x x f cos 1)(+='，直至

5

110

1||-+?<-k k x x 时，结束迭代；然后再使用基于牛顿法的Steffensen 加速法进

行计算，直至51101||-+?<-k k x x 时，结束迭代。其迭代结果与迭代次数如下表所示（注N1为牛顿法迭代次数，N2为基于牛顿法Steffensen 加速法迭代次数）：

2.对方程sin x =0，使用牛顿法时，令x x f sin )(=，使用牛顿法计算，直至

5

110

1||-+?<-k k x x 时，结束迭代；然后依据Steffensen 加速法进行编程计算，直

至51101||-+?<-k k x x 时，结束迭代。其迭代结果与迭代次数如下表所示：

附：

第一问牛顿法

#include

#include

#define EPS 1.0e-5

int main()

{

double x=0.0;

double y=0.0;

int n=0;

cout<<"请输入初值X0"<

cin>>x;

do

{

y=x;

x=x-(x+sin(x)-1)/(1+cos(x));

n++;

}

while(fabs(x-y)>=EPS);

cout<<"X"<

cout<<"迭代次数为N="<

}

第一问牛顿-Steffensen法

#include

#include

#define EPS 1.0e-5

int main()

{

double x=0.0;

double y=0.0;

double z=0.0;

double a=0.0;

int n=0;

cout<<"请输入初值x0"<

cin>>x;

do

{

a=x;

y=x-(x+sin(x)-1)/(1+cos(x));

z=y-(y+sin(y)-1)/(1+cos(y));

x=x-(y-x)*(y-x)/(z-2*y+x);

n++;

}

while(fabs(x-a)>=EPS);

cout<<"x"<

cout<<"迭代次数N="<

}

第二问牛顿法

#include

#include

#define EPS 1.0e-5

int main()

{

double x=0.0;

double y=0.0;

int n=0;

cout<<"请输入初值X0"<

cin>>x;

do

{

y=x;

x=x-tan(x);

n++;

}

while(fabs(x-y)>=EPS);

cout<<"X"<

cout<<"迭代次数为N="<

}

第二问牛顿-Steffensen法

#include

#include

int main()

{

double x=0.0;

double y=0.0;

double z=0.0;

double a=0.0;

int n=0;

cout<<"请输入初值x0"<

cin>>x;

do

{

a=x;

y=x-tan(x);

z=y-tan(y);

x=x-(y-x)*(y-x)/(z-2*y+x);

n++;

}

while(fabs(x-a)>=1.0E-05);

cout<<"x"<

cin>>n;

cout<<"请输入数组A:";

for(i=0;i

for(j=0;j

cin>>a[i][j];

cout<<"请输入数组B:";

for(i=0;i

cin>>b[i];

cout<<"请输入初始X:";

for(i=0;i

cin>>x[i];

do

{for(i=0;i

{for(j=0;j

sum=sum+a[i][j]*x[j];

for(j=i+1;j

sum=sum+a[i][j]*x[j];

y[i]=(b[i]-sum)/a[i][i];

sum=0;

}

c=0;

for(i=0;i

if(c

c=fabs(x[i]-y[i]);

if(c>0.00001)

c=0;

else s=1;

for(i=0;i

x[i]=y[i];

k++;

}while(s!=1);

for(i=0;i

cout<<"解为："<

cout<<"N="<

return 0;

}

高斯－赛德尔迭代法

#include

#include

using namespace std ;

int main( )

{

double a[10][10],b[10],sum=0.0,x[10],y[10],c,s=0.0,z[10];

int i,k=0,n,j;

cout<<"请输入维数:";

cin>>n;

cout<<"请输入数组A:";

for(i=0;i

for(j=0;j

cin>>a[i][j];

cout<<"请输入数组B:";

for(i=0;i

cin>>b[i];

cout<<"请输入初始X:";

for(i=0;i

cin>>x[i];

do

{for(i=1;i

z[i]=x[i];

for(i=0;i

{for(j=0;j

sum=sum+a[i][j]*x[j];

for(j=i+1;j

sum=sum+a[i][j]*x[j];

y[i]=(b[i]-sum)/a[i][i];

sum=0;

x[i]=y[i];

}

c=0;

for(i=0;i

if(c

c=fabs(z[i]-y[i]);

if(c>0.00001)

c=0;

else s=1;

for(i=0;i

z[i]=y[i];

k++;

}while(s!=1);

for(i=0;i

cout<<"解为："<

cout<<"N="<

return 0;

}

3. 结果分析

问题（1）中因系数矩阵严格对角占优，所以，无论采用Jacobi迭代法还是Causs-Seidel迭代法，其迭代过程均收敛。从迭代次数来看，方程组Ax=b1中Jacobi 迭代法使用迭代次数20次，而Causs-Seidel迭代法使用迭代次数12次；方程组Ax=b2中Jacobi迭代法使用迭代次数26次，而Causs-Seidel迭代法使用迭代次数17次。可明显看出，在方程组绝对收敛的条件下，Causs-Seidel迭代法比Jacobi 迭代法收敛速度快。

问题（2）中，使用Jacobi迭代法方法时，其方程组不收敛，而采用Causs-Seidel 迭代法迭代时，方程Ax=b1迭代37次，得到准确解；方程Ax=b2迭代44次，

得到准确解。这是因为，Jacobi 迭代矩阵为：0

8

.08

.08.008

.08

.08.00------ ，经计算其特征值λ为-1.6，0.8，0.8，其谱半径为=)(B ρ 1.6>1，所以其迭代矩阵发散。而

Causs-Seidel 迭代矩阵为：768

.0128

.00

16.064.00

8

.08.00

---，经计算其特征值λ为0，0.7040 + 0.1280i ，0.704280i ，其谱半径为=)(B ρ0.7155<1，所以采用Causs-Seidel 迭代法时，方程收敛。

问题（3）中，采用Jacobi 迭代法和Causs-Seidel 迭代法其方程均不收敛，这是因为，Jacobi 迭代矩阵为：

7

30-，经计算其特征值λ为-4.5826i ，4.5826i

其谱半径为=)(B ρ 4.5826>1,迭代方法发散。Causs-Seidel 迭代矩阵为：21

30--，

经计算其特征值λ为0,21，其谱半径为=)(B ρ,21>1,迭代方法发散。如果将方程系数矩阵两行进行交换，则为严格对角优势阵，无论采用Jacobi 迭代法还是Causs-Seidel 迭代法，其矩阵均收敛。

三 Ru n ge-Kutt a 4阶算法中不同步长对稳定区间的作用

1. 计算题目

用Ru n ge-Kutt a 4阶算法对初值问题y /=-20*y ，y (0)=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。 注：此方程的精确解为：y =e -20x

2. 计算过程与结果

使用C++中函数调用功能，对Ru n ge-Kutt a 4阶算法，分别采用步长h=0.1,0.2计算，得：

附：

Ru n ge-Kutt a 4阶算法

#include

void main(){

double x0,y0,h,x;

double x1,y1,k1,k2,k3,k4;

cout<<"请输入初值x0,y0:";

cin>>x0>>y0;

cout<<"请输入最大X:";

cin>>x;

cout<<"请输入步长的大小:"; cin>>h;

double f(double x,double y);

for(int n=1;n<=x/h;n++,x0=x1,y0=y1)

{

x1=x0+h;

k1=f(x0,y0);

k2=f(x0+h/2,y0+(h/2)*k1);

k3=f(x0+h/2,y0+(h/2)*k2);

k4=f(x0+h,y0+h*k3);

y1=y0+h*(k1+2*k2+2*k3+k4)/6;

}

cout<<"y1="<

}

double f(double x,double y){

return -20*y;

}

3. 结果分析

λ,h=0.2时，λh=-4,而Ru n ge-Kutt a 4阶算法的绝对稳定区间是=

-

20

[-2.78,0],故h=0.2时计算不稳定；而h=0.1时，λh=-2，在绝对稳定区间内，计算稳定，结果可靠。

总结

通过这次上机练习，让我对数值分析所介绍的迭代求解方法及其理论有了更深层次的理解，了解了各种方法之间的优缺点，并且认识到了自己在以前的学习中所存在的问题，及时的修补了自己知识上的漏洞。同时也提高了我在编写程序上的熟练程度，复习了以前所学的知识。本次上机实习是非常有意义的，给予我在编程与数值分析学习很大的帮助。