西南交通大学研究生数值分析作业

数值分析上机报告指导教师：赵海良班级：姓名：学号：电话：2011年12月序随着计算机技术的迅速发展，数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。

要解决工程问题，往往需要处理很多数学模型，不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：解法是否稳定、是否收敛及熟练的速度等。

由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如MATLAB，C++，VB，JA V A的辅助软件来解决，得到一个满足误差限的解。

本文所计算题目，均采用C++编程。

C++是一种静态数据类型检查的、支持多重编程范式的通用程序设计语言。

它支持过程化程序设计、数据抽象、面向对象程序设计、制作图标等等泛型程序设计等多种程序设计风格，在实际工程中得到了广泛应用，对解决一些小型数学迭代问题，C++软件精度已满足相应的精度。

本文使用C++对牛顿法、牛顿-Steffensen法对方程求解，对雅格比法、高斯－赛德尔迭代法求解方程组迭代求解，对Ru n ge-Kutt a 4阶算法进行编程，并通过实例求解验证了其可行性，并使用不同方法对计算进行比较，得出不同方法的收敛性与迭代次数的多少，比较不同方法之间的优缺性，比较各种方法的精确度和解的收敛速度。

目录第一章牛顿法和牛顿-Steffensen法迭代求解的比较 (1)1.1 计算题目 (1)1.2 计算过程和结果 (1)1.3 结果分析 (2)第二章 Jacobi迭代法与Causs-Seidel迭代法迭代求解的比较 (2)2.1 计算题目 (2)2.2 计算过程与结果 (2)2.3 结果分析 (3)第三章 Ru n ge-Kutt a 4阶算法中不同步长对稳定区间的作用 (4)3.1 计算题目 (4)3.2 计算过程与结果 (4)3.3 结果分析 (4)总结 (5)附件 (6)附件 1（1.1第一问牛顿法） (6)附件 2（1.1第一问牛顿-Steffensen法） (6)附件 3（1.1第二问牛顿法） (6)附件 4（1.1第二问牛顿-Steffensen法） (7)附件 5（2.1 Jacobi迭代法） (7)附件 6（2.1Causs-Seidel迭代法） (8)附件 7（3.1 Ru n ge-Kutt a 4阶算法） (9)第一章 牛顿法和牛顿-Steffensen 法迭代求解的比较1.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进行计算。

分析其中遇到的现象与问题。

1.2计算过程和结果1.对方程ln(x +sin x )=0，其导数有些复杂，我们可以对其进行变形，即求解x+sinx=1的解。

使用牛顿法，令1sin )(-+=x x x f ，则x x f cos 1)(+='，直至51101||-+⨯<-k k x x 时，结束迭代；然后再使用基于牛顿法的Steffensen 加速法进行计算，直至51101||-+⨯<-k k x x 时，结束迭代。

其迭代结果与迭代次数如下表所示（注N1为牛顿法迭代次数，N2为基于牛顿法Steffensen 加速法迭代次数）：2.对方程sin x =0，使用牛顿法时，令x x f sin )(=，使用牛顿法计算，直至51101||-+⨯<-k k x x 时，结束迭代；然后依据Steffensen 加速法进行编程计算，直至51101||-+⨯<-k k x x 时，结束迭代。

其迭代结果与迭代次数如下表所示：1.3结果分析（1）从牛顿与Steffensen 加速法可以看出，牛顿—Steffensen 加速法迭代速度明显快于牛顿加速法，从计算结果可以看出，x=0.510973处取得精确值，而在x0=0.1,1处收敛速度较快，说明牛顿法与基于牛顿法Steffensen 加速法在单根附近有较快的收敛速度，迭代法是否收敛，与初始近似值x0的好坏很有关系。

（2）从sinx=0使用两种方法可以看出，Steffensen 加速法迭代速度明显快于牛顿加速法，且牛顿法和基于牛顿法Steffensen 加速法在单根附近有较快的收敛速度，迭代法是否收敛，与初始近似值x0的好坏很有关系，如在x0=1.6处，得到的收敛解很大，这是因为0)/1806.1cos()6.1(≈⨯='πf ，很难说明f(x)在此处是否发散或者是否收敛。

第二章 Jacobi 迭代法与Causs-Seidel 迭代法迭代求解的比较2.1计算题目用雅格比法与高斯－赛德尔迭代法解下列方程组Ax =b ，研究其收敛性，上 机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。

(1)A 行分别为A 1=[6,2,-1],A 2=[1,4,-2],A 3=[-3,1,4]； b 1=[-3,2,4]T , b 2=[100,-200,345]T ,(2) A 行分别为A 1=[1,0,8,0.8],A 2=[0.8,1,0.8],A 3=[0.8,0.8,1]；b 1=[3,2,1] T , b 2=[5,0,-10]T ,(3)A 行分别为A 1=[1,3],A 2=[-7,1]；b =[4,6]T ,2.2计算过程与结果1.采用C++进行编程，分别采用用雅格比法与高斯－赛德尔迭代法对问题（1）进行迭代，x0= [0,0,0]T 为初始值。

2.采用C++进行编程，分别采用用雅格比法与高斯－赛德尔迭代法对问题（2）进行迭代，x0= [0,0,0]T 为初始值。

3.采用C++进行编程，分别采用用雅格比法与高斯－赛德尔迭代法对问题（3）进行迭代，x0= [0,0]T 为初始值。

2.3 结果分析问题（1）中因系数矩阵严格对角占优，所以，无论采用Jacobi 迭代法还是Causs-Seidel 迭代法，其迭代过程均收敛。

从迭代次数来看，方程组Ax =b 1中Jacobi 迭代法使用迭代次数20次，而Causs-Seidel 迭代法使用迭代次数12次；方程组Ax =b 2中Jacobi 迭代法使用迭代次数26次，而Causs-Seidel 迭代法使用迭代次数17次。

可明显看出，在方程组绝对收敛的条件下，Causs-Seidel 迭代法比Jacobi 迭代法收敛速度快。

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

这是因为，Jacobi 迭代矩阵为：08.08.08.008.08.08.00------ ，经计算其特征值λ为-1.6，0.8，0.8，其谱半径为=)(B ρ 1.6>1，所以其迭代矩阵发散。

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

问题（3）中，采用Jacobi 迭代法和Causs-Seidel 迭代法其方程均不收敛，这是因为，Jacobi 迭代矩阵为：0730-，经计算其特征值λ为-4.5826i ，4.5826i其谱半径为=)(B ρ 4.5826>1,迭代方法发散。

Causs-Seidel 迭代矩阵为：21030--，经计算其特征值λ为0,21，其谱半径为=)(B ρ,21>1,迭代方法发散。

如果将方程系数矩阵两行进行交换，则为严格对角优势阵，无论采用Jacobi 迭代法还是Causs-Seidel 迭代法，其矩阵均收敛。

第三章 Ru n ge-Kutt a 4阶算法中不同步长对稳定区间的作用3.1计算题目用Ru n ge-Kutt a 4阶算法对初值问题y /=-20*y ，y (0)=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。

注：此方程的精确解为：y =e -20x3.2 计算过程与结果使用C++中函数调用功能，对Ru n ge-Kutt a 4阶算法，分别采用步长h=0.1,0.2计算，得：3.3 结果分析从上述计算结果可以看出，当步长h=0.1时，其结果精确，在x=1.0处，误差为1.693293885E-05而步长为h=0.2时，其迭代过程不收敛，说明用Ru n ge-Kutt a4阶算法是否收敛与步长h有关。

因当λ=-20，h=0.2时，λh=-4，不在Ru n ge-Kutt a 4阶算法的绝对稳定区间[-2.785,0]之间，计算不稳定。

而当h=0.1时，λh=-2，在绝对稳定区间内，计算稳定，结果可靠。

总结通过这次上机练习，让我对数值分析所介绍的迭代求解方法及其理论有了更深层次的理解，通过对比，使我认识到各种方法的优缺点，并使我认识到学习数值分析中的不足，时的修补了自己知识上的漏洞，以便在复习与以后的工作中得以加强。

同时，通过本次上机实习，复习了本科阶段所学习的C++编程，提高了编程的熟练度，并将编程能应用到实际的工作中去，解决实际遇到的工程问题以及数学问题，加强了两门学科之间的联系。

总之，本次上机实习使我获益匪浅，给予我编程与数值分析学习很大的帮助。

附件附件 1（1.1第一问牛顿法）#include<iostream.h>#include<math.h>#define EPS 1.0e-5int main(){double x=0.0;double y=0.0;int n=0;cout<<"请输入初值X0"<<endl; cin>>x;do{y=x;x=x-(x+sin(x)-1)/(1+cos(x)); n++;}while(fabs(x-y)>=EPS);cout<<"X"<<n<<"="<<x<<endl; cout<<"迭代次数为N="<<n<<endl; }附件 2（1.1第一问牛顿-Steffensen法）#include<iostream.h>#include<math.h>#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"<<endl; 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"<<n<<"="<<x<<endl; cout<<"迭代次数N="<<n<<endl; }附件 3（1.1第二问牛顿法）#include<iostream.h>#include<math.h>#define EPS 1.0e-5int main(){double x=0.0;double y=0.0;int n=0;cout<<"请输入初值X0"<<endl; cin>>x;do{y=x;x=x-tan(x);n++;}while(fabs(x-y)>=EPS);cout<<"X"<<n<<"="<<x<<endl; cout<<"迭代次数为N="<<n<<endl; }附件 4（1.1第二问牛顿-Steffensen法）#include<iostream.h>#include<math.h>#define EPS 1.0e-6int main(){double x=0.0;double y=0.0;double z=0.0;double a=0.0;int n=0;cout<<"请输入初值x0"<<endl; 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)>=EPS);cout<<"x"<<n<<"="<<x<<endl; cout<<"迭代次数N="<<n<<endl; }附件 5（2.1 Jacobi迭代法）#include<iostream.h>#include<math.h>#define EPS 1.0e-5int main(){int N,bk=0;int i,j,l=0;float fnum;cout<<"请输入N值";cin>>N;double a[12][12],b[10],c[12],x[12], y[12], t=0,temp=0;cout<<"请输入矩阵A"<<endl; for(i=0;i<N;i++)for(j=0;j<N;j++){cout<<"a"<<i<<j<<"=";cin>>a[i][j];}{for(j=0;j<N;j++)cout<<a[i][j]<<" ";cout<<endl;}cout<<"请输入b"<<endl;for(i=0;i<N;i++){cout<<"b"<<i<<"="<<" ";cin>>b[i];}cout<<"请输入x0"<<endl;for(i=0;i<N;i++)cin>>x[i];for(i=0;i<N;i++){c[i]=1/a[i][i];a[i][i]=0;}while(bk!=1){for(i=0;i<N;i++){for(j=0;j<N;j++)temp=temp-a[i][j]*x[j];y[i]=c[i]*(temp+b[i]);temp=0;}for(i=0;i<N;i++)fnum=float(fabs(x[i]-y[i]));if(t<fnum) t=fnum;}if(t<=0.00001) bk=1;t=0;for(i=0;i<N;i++)x[i]=y[i];l++;}for(i=0;i<N;i++)cout<<x[i]<<" "<<endl;cout<<"迭代次数N为"<<l;}附件 6（2.1Causs-Seidel迭代法）#include<iostream.h>#include<math.h>#define EPS 1.0e-5int main(){int N,bk=0;int i,j,l=0;float fnum;cout<<"请输入N值";cin>>N;double a[12][12],b[10],c[12],x[12],y[12],z[12],t=0,temp=0;cout<<"请输入矩阵A"<<endl;for(j=0;j<N;j++){cout<<"a"<<i<<j<<"="; cin>>a[i][j];}for(i=0;i<N;i++){for(j=0;j<N;j++)cout<<a[i][j]<<" ";cout<<endl;}cout<<"请输入b"<<endl;for(i=0;i<N;i++){cout<<"b"<<i<<"="<<" ";cin>>b[i];}cout<<"请输入x0"<<endl;for(i=0;i<N;i++)cin>>x[i];for(i=0;i<N;i++){c[i]=1/a[i][i];a[i][i]=0;}while(bk!=1){for(i=0;i<N;i++)z[i]=x[i];for(i=0;i<N;i++)for(j=0;j<N;j++)temp=temp-a[i][j]*z[j];y[i]=c[i]*(temp+b[i]);z[i]=y[i];temp=0;}for(i=0;i<N;i++){fnum=float(fabs(x[i]-y[i]));if(t<fnum) t=fnum;}if(t<=0.00001) bk=1;t=0;for(i=0;i<N;i++)x[i]=y[i];l++;}for(i=0;i<N;i++)cout<<x[i]<<" "<<endl;cout<<"迭代次数N为"<<l;}附件 7（3.1 Ru n ge-Kutt a 4阶算法）#include <iostream.h>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="<<y1;}double f(double x,double y){return -20*y;}。