数值分析上机报告前言随着计算机技术的高速发展,越来越多的科技工作者使用计算机进行科学研究和解决工程技术问题。
数值分析(或计算方法)课程的内容是科学工程计算的必备知识,已经成为众多理工科大学生、研究生的必修课程,越来越受到重视。
由于工程实际中所遇到的数学模型求解过程迭代次数很多,计算量很大,所以需要借助很多编程软件来解决,得到一个满足误差限的解。
本文所计算题目,均采用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)(+=',直至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 时,结束迭代。
其迭代结果与迭代次数如下表所示:附:第一问牛顿法#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;}第一问牛顿-Steffensen法#include<iostream.h>#include<math.h>#define EPS 1.0e-5int 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;}第二问牛顿法#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;}第二问牛顿-Steffensen法#include<iostream.h>#include<math.h>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-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"<<n<<"="<<x<<endl;cout<<"迭代次数N="<<n<<endl;}3. 结果分析从结果对比我们可发现牛顿—Steffensen 加速法比牛顿法要收敛的快,牛顿法对于初值的选取特别重要,比如第(1)问中的初值为4的情况,迭代次数算了40次,远大于其余初值的情况;在第(2)问中的初值为1.6的情况,收敛解得31.4159,分析其原因应该是x x f c os )('=,x0=1.62π≈,0)('≈x f ;在牛顿—Steffensen 加速法第(1)问中x=1.5、第二问x=0和3的情况,无收敛解。
二 Jacobi 迭代法与Causs-Seidel 迭代法迭代求解的比较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. 计算过程与结果(1)x0= [0,0, 0]T 为初始值,N 为迭代次数,0.00001为误差精度,X 为收敛解附:雅格比法#include <iostream>#include <cmath>using namespace std ;int main( ){double a[10][10],b[10],sum=0.0,x[10],y[10],c,s=0.0;int i,k=0,n,j;cout<<"请输入维数:";cin>>n;cout<<"请输入数组A:";for(i=0;i<n;i++)for(j=0;j<n;j++)cin>>a[i][j];cout<<"请输入数组B:";for(i=0;i<n;i++)cin>>b[i];cout<<"请输入初始X:";for(i=0;i<n;i++)cin>>x[i];do{for(i=0;i<n;i++){for(j=0;j<i;j++)sum=sum+a[i][j]*x[j];for(j=i+1;j<n;j++)sum=sum+a[i][j]*x[j];y[i]=(b[i]-sum)/a[i][i];sum=0;}c=0;for(i=0;i<n;i++)if(c<fabs(x[i]-y[i]))c=fabs(x[i]-y[i]);if(c>0.00001)c=0;else s=1;for(i=0;i<n;i++)x[i]=y[i];k++;}while(s!=1);for(i=0;i<n;i++)cout<<"解为:"<<x[i]<<endl;cout<<"N="<<k<<endl;return 0;}高斯-赛德尔迭代法#include <iostream>#include <cmath>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<n;i++)for(j=0;j<n;j++)cin>>a[i][j];cout<<"请输入数组B:";for(i=0;i<n;i++)cin>>b[i];cout<<"请输入初始X:";for(i=0;i<n;i++)cin>>x[i];do{for(i=1;i<n;i++)z[i]=x[i];for(i=0;i<n;i++){for(j=0;j<i;j++)sum=sum+a[i][j]*x[j];for(j=i+1;j<n;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<n;i++)if(c<fabs(z[i]-y[i]))c=fabs(z[i]-y[i]);if(c>0.00001)c=0;else s=1;for(i=0;i<n;i++)z[i]=y[i];k++;}while(s!=1);for(i=0;i<n;i++)cout<<"解为:"<<x[i]<<endl;cout<<"N="<<k<<endl;return 0;}3. 结果分析问题(1)中因系数矩阵严格对角占优,所以,无论采用Jacobi迭代法还是Causs-Seidel迭代法,其迭代过程均收敛。