数值分析实践报告
东北大学软件学院
实验一:
1. 实验目的
1、通过上机计算体会迭代法求解线性方程组的特点,并能和消去法比较;
2、运用所学的迭代法算法,解决各类线性方程组,编出算法程序;
3、体会上机计算时,终止步骤<或k>(予给的迭代次数),对迭代法敛散性的意义;
4、体会初始解x,松弛因子的选取,对计算结果的影响。
2. 实验环境
使用平台:Microsoft Visual C++
使用语言:C++
3. 实验关键代码
(1)存放系数矩阵
double **A=new double*[n];
for(i=0;i<n;i++)
A[i]=new double[n];
cout<<"请输入4阶系数矩阵A:"<<endl;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
cin>>A[i][j];
(2)迭代算法
void Diedai(int n,double **A,double b[],double e,double w,double X[]){
int i,j,flag;
double *t=new double[n];
double *Y=new double[n];
double *X_T=new double[n];
cout<<"矩阵为:"<<endl;
cout<<"**************************"<<endl;
for(i=0;i<n;i++){
for(j=0;j<n;j++)
cout<<A[i][j]<<'\t';
cout<<endl;}
cout<<"**************************"<<endl;
for(int k=0;k<2000;k++){
flag=0;
for(i=0;i<n;i++){
Y[i]=0;
for(j=0;j<n;j++){
if(i!=j)
Y[i]=Y[i]+A[i][j]*X[j];}
X_T[i]=(1-w)*X[i]+w*(b[i]-Y[i])/A[i][i];
t[i]=fabs(X[i]-X_T[i]);
X[i]=X_T[i];}
for(i=0;i<n;i++)
if(t[i]<=e)
flag++;
if(flag==n){
cout<<"迭代次数为:"<<endl;
cout<<(k+1)<<endl;
break;}}
if(k==2000)
cout<<"超出最大迭代次数 "<<20##<<"!"<<endl;
delete []t;
delete []X_T;
delete []Y;
}
注:以上仅是关键算法的代码,不是完整代码。
4. 实验结果
依次输入7个超松弛因子,结果如下:
因子为0.15时,需要迭代44次;
因子为0.45时,需要迭代17次;
因子为0.75时,需要迭代10次;
因子为1.0时,需要迭代6次;
因子为1.25时,需要迭代6次;
因子为1.5时,需要迭代32次;
因子为1.75时,需要迭代超过2000次。
超松弛迭代法收敛速度的快慢与松弛因子的选择有密切关系,迭代次数随着超松弛因子增加先减少再增加。
实验二:
1. 实验目的
1、通过上机计算体会迭代法求解线性方程组的特点,并能和消去法比较;
2、运用所学的迭代法算法,解决各类线性方程组,编出算法程序;
3、体会上机计算时,终止步骤<或k>(予给的迭代次数),对迭代法敛散性的意义;
4、体会初始解x,松弛因子的选取,对计算结果的影响。
2. 实验环境
使用平台:Microsoft Visual C++
使用语言:C++
3. 实验关键代码
(1) 系数矩阵的定义
cout<<"<<请输入矩阵的规格(?X?)>>:";
cin>>matrixNum;
matrixA=allocMem(matrixNum*matrixNum);
matrixD=allocMem(matrixNum*matrixNum);
matrixL=allocMem(matrixNum*matrixNum);
matrixU=allocMem(matrixNum*matrixNum);
B=allocMem(matrixNum*matrixNum);
f=allocMem(matrixNum);
x=allocMem(matrixNum);
xk=allocMem(matrixNum);
b=allocMem(matrixNum);
cout<<endl<<endl<<endl<<"<<请输入矩阵中各元素值(为 "<<matrixNum<<"*"<<matrixNum<<",共计 "<<matrixNum*matrixNum<<" 个元素)"<<">>:"<<endl;
for(i=0;i<matrixNum;i++){
cout<<"请输入矩阵中第 "<<i+1<<" 行的 "<<matrixNum<<" 个元素:";
for(j=0;j<matrixNum;j++)
cin>>*(matrixA+i*matrixNum+j);
}
(2)高斯列主元素消去法
void GaussLineMain(double* A,double* x,double* b,int num)
{
int i,j,k;
for(i=0;i<num-1;i++){
double lineMax=fabs(*(A+i*num+i));
int lineI=i;
for(j=i;j<num;j++)
if(fabs(*(A+j*num+i))>fabs(lineMax)) lineI=j;
for(j=i;j<num;j++){
lineMax=*(A+i*num+j);
*(A+i*num+j)=*(A+lineI*num+j);
*(A+lineI*num+j)=lineMax;
//b中对应元素做交换
lineMax=*(b+i);
*(b+i)=*(b+lineI);
*(b+lineI)=lineMax;
}
if(*(A+i*num+i)==0) continue;
for(j=i+1;j<num;j++){
double temp=-*(A+j*num+i)/(*(A+i*num+i));
for(k=i;k<num;k++)
*(A+j*num+k)+=temp*(*(A+i*num+k));
*(b+j)+=temp*(*(b+i));
}
}
double determinantA=1;
for(i=0;i<num;i++)
determinantA*=*(A+i*num+i);
if(determinantA==0){
cout<<endl<<endl<<"通过计算,矩阵A的行列式为|A|=0,即A没有唯一解。\n程序退出..."<<endl<<endl;
exit(0);
}
for(i=num-1;i>=0;i--){
for(j=num-1;j>i;j--)
*(b+i)-=*(A+i*num+j)*(*(x+j));
*(x+i)=*(b+i)/(*(A+i*num+i));
}
}
(3)雅可比迭代算法
void Jacobi(double* x,double* xk,double* B,double* f,int num,int time)
{
int t=1,i,j;
while(t<=time){
for(i=0;i<num;i++){
double temp=0;
for(j=0;j<num;j++)
temp+=*(B+i*num+j)*(*(x+j));
*(xk+i)=temp+*(f+i);
}
for(i=0;i<num;i++)
*(x+i)=*(xk+i);
t++;
}
}
(4)分配内存空间
double* allocMem(int num)
{
double *head;
if((head=new double[num])==NULL){
cout<<"内存空间分配失败,程序终止运行!"<<endl;
exit(0);
}
return head;
}
注:以上仅是关键算法的代码,不是完整代码。
4. 实验结果
结果如图:
则答案如下:
第二问:稳定种群数量为(8481,,2892,1335,601,141)
第三问:x迭代后x5=-259,所以无法达到。
5. 实验总结
1、 超松弛迭代法收敛速度的快慢与松弛因子的选择有密切关系;
2、 迭代次数随着超松弛因子增加先减少再增加;
3、 在实际计算中,可采用试算方法来确定较好的松弛因子。
第二篇:数值分析第一次上机练习实验报告
数值分析第一次上机练习实验报告
——Lagrange插值与三次样条插值
一、问题的描述
设, ,取,.试求出10次Lagrange插值多项式和三次样条插值函数(采用自然边界条件),并用图画出,, .
二、方法描述——Lagrange插值与三次样条插值
我们取,,通过在点的函数值来对原函数进行插值,我们记插值函数为,要求它满足如下条件:
(1)
我们在此处要分别通过Lagrange插值(即多项式插值)与三次样条插值的方法对原函数进行插值,看两种方法的插值结果,并进行结果的比较。
10次的Lagrange插值多项式为:
(2)
其中:
以及
我们根据(2)进行程序的编写,我们可以通过几个循环很容易实现函数的Lagrange插值。
理论上我们根据区间上给出的节点做出的插值多项式近似于,而多项式的次数越高逼近的精度就越好。但实际上并非如此,而是对任意的插值节点,当的时候不一定收敛到;而是有时会在插值区间的两端点附近会出现严重的偏离的现象,即所谓的Runge现象。因此用高次插值多项式近似的效果并不总是好的,因而人们通常在选择插值方式的时候不用高次多项式插值,而用分段低次插值,而这样的插值效果往往是非常好的,能够克服高次多项式插值的弱点,达到令人满意的效果。
分段低次插值包括分段线性插值、分段三次Hermite插值、三次样条插值等。前两种插值函数都具有一致收敛性,但是光滑性较差,而在实际问题中我们往往要求函数具有二阶光滑度,即有二阶连续导数。而对第三种插值方式,我们得到的是一个样条曲线,它是由分段三次曲线拼接而成,在连接点(即样点)上二阶导数连续。
我们记三次样条插值函数为,它在每个小区间上是三次函数,因此在每个区间上需要确定4个参数,总共有10个小区间,因此共需确定40个未知参数。首先我们有插值条件:
(3)
其次在每个节点上满足连续性条件:
(4)
此外在端点处满足自然边界条件:
(5)
我们假设。则在每个小区间上:
(6)
其中:
及
我们利用边界条件(3)(4)(5)可以得到:
(7)
其中:
以及
两端点处的边界条件为:
(8)
将边界条件写成矩阵形式为:
(9)
其中根据自然边界条件(8)有:
我们解方程(9)就可以得到,将他们代入(6)就可以得到各段区间上的的值。
三、方案设计
我们通过编写Matlab程序来进行10次Lagrange插值与三次样条插值的工作。在我们的程序文件中interplotion.m文件是主程序文件;L10.m文件是计算10次Lagrange插值多项式的子程序文件,给它任一个,此程序将返回的值;Mspline.m是根据(9)计算各节点二阶导数值的子程序文件,它将会返回在自然边界条件下的各节点的二阶导数值;然后spline.m是根据以及(6)计算三次样条插值函数的子程序文件。然后运行主程序将给出三幅曲线图,分别是与曲线,与曲线,以及、与三条曲线共同画在一幅图上得到的图象。
解决这个问题的思路很简单,按部就班的来就可以。首先我们计算各节点上的函数值以备后用,然后调用Mspline.m计算。随后我们给出一系列的值,计算,并分别调用L10.m与spline.m分别计算与。然后根据我们得到的数据绘图观察插值结果。具体程序的实现可参见所给程序的相关注释。
四、计算结果及其分析
下面是我们根据程序计算结果得到的数据,其中分别给出了在各典型处的的原函数的值、Lagrange插值结果与样条插值结果;以及绝对误差和,相对误差,。由于在两端点处进行Lagrange插值插值的时候可能出现Runge现象,因此我们在两端点附近多给了几个点的数据。
尽管从数据可以看出一些端倪,但是通过图象我们更能清楚地看到最终插值结果的定性情况。首先我们给出与曲线:
其中蓝色的曲线代表曲线,绿色的曲线代表曲线。可见此时两者之间具有很大的差别,尤其在端点附近会出现严重的偏离的现象,即出现了所谓的Runge现象。而此时曲线与我们用样条插值得到的的曲线为:
其中蓝色的曲线代表曲线,绿色的曲线代表曲线,可见两条曲线几乎完全重合,与符合的很好。
上面我们由曲线定性看到的结论也可以通过表中的数据定量的看出。
五、结论
插值方法中最基本的是多项式插值,而我们可以通过Lagrange多项式来方便的实现这种插值方式。理论上我们根据给定区间上的给定的节点做出的插值多项式近似于,而多项式的次数越高逼近的精度就越好。但实际上对任意的插值节点,当的时候不一定收敛到;而是有时会在插值区间的两端点附近会出现严重的偏离的现象,即所谓的Runge现象。因此用高次插值多项式近似的效果并不总是好的,而我们通过本次试验中的实际计算发现对本次试验中的函数确实出现了Runge现象,插值结果很不令人满意;我们转而采用分段的三次样条插值,得到了非常好的插值效果。