数值计算方法实验报告(含所有)

本科实验报告
课程名称：计算机数值方法实验项目：计算机数值方法实验实验地点：
专业班级：学号：
学生姓名：xxx
指导教师：xxx
太原理工大学学生实验报告
学院名称软件学院专业班级1217班学号201200xxxx 学生姓名xx 实验日期2014.05.21 成绩
课程名称数值计算方法实验题目实验一方程求解
一、实验目的和要求
熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。

选择上述方法中的两种方法求方程：二分法f(x)=x3+4x2-10=0在[1,2]内的一个实根，且要求满足精度|x*-x n|<0.5×10-5
二、主要设备
笔记本 HP ProBook 6470b 一台
编译软件：VC++6.0
三、实验内容和原理
函数f(x)在区间（x，y）上连续，先在区间（x，y）确定a与b，若f(a)，f(b)异号，说明在区间(a，b)内存在零点，然后求f[(a+b)/2]。

假设F(a)<0,F(b)>0,a<b，
①如果f[(a+b)/2]=0，该点即为零点；
②如果f[(a+b)/2]<0，则区间（(a+b)/2，b)内存在零点，(a+b)/2≥a；
③如果f[(a+b)/2]>0，则区间(a,(a+b)/2)内存在零点，(a+b)/2≤b；
返回①重新循环，不断接近零点。

通过每次把f(x)的零点所在区间收缩一半的方法，使区间内的两个端点逐步逼近函数零点，最终求得零点近似值。

四、操作方法与实验步骤
1. 二分法：
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
int main()
{
double a=1.0, b=2.0;
double x,s;
printf(" An\t\tBn\t\tF(Xn)\n");
while(1)
{
x=(a+b)/2;
s=pow(x,3)+4*x*x-10;
if (-0.000005 < s && s < 0.000005)
{
break;
}
else if(s < 0)
{
a=x;
}
else if(s > 0)
{
b=x;
}
printf("%f\t%f\t%f\n",a,b,s);
}
printf("X的值为：%f\n",x);
printf("误差：\t%f\n",s);
return 0;
}
2. 割线法：
#include"stdio.h"
#include"math.h"
int main()
{
float c,a=1.0,b=2.0;
printf("每次得到的X的近似值：\n");
while(1)
{
c=b-(b*b*b+4*b*b-10)*(b-a)/(b*b*b+4*b*b-(a*a*a+4*a*a));
if(fabs(b-c)<0.5*0.00001)break;
b=c;
printf("%f\n",b);
}
printf("X的值为：%f\n",c);
}
五、实验结果与分析
二分法割线法分析：
由程序知，使用二分法和割线法均能计算出方程的根，但利用割线法要比二分法计算的次数少，并且能够较早的达到精度要求。

相比之下，割线法程序代码量较少，精简明了。

六、讨论、心得
本次数值计算方法程序设计实验从习题练习中跳脱出来，直接面对实用性较强的程序代码编写。

效果很好，不仅加深对二分法、割线法的理解，还加强了实际用运能力。

将理论知识成功地转化成实践结果。

实验地点虎峪校区致远楼B401指导教师xx
太原理工大学学生实验报告
学院名称 软件学院 专业班级 1217班 学号 201200xxxx
学生姓名 xx 实验日期 2014.05.28
成绩
课程名称
数值计算方法
实验题目
实验二 线性方程组的直接解法
一、实验目的和要求
合理利用Gauss 消元法、LU 分解法、追赶法求解下列方程组： ①
⎥⎥
⎥⎦
⎤
⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡1381414
2
21032
1321x x x
② ⎥
⎥⎥⎥⎦⎤⎢⎢⎢
⎢⎣⎡=⎥⎥
⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢
⎣
⎡
--⨯-2178.4617.5911
2
1
259
2.1121130
.6291
.51314.59103.0432115
x x x x
③
⎥
⎥
⎥
⎥⎦⎤
⎢⎢⎢⎢⎣⎡----=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡377220116
12
6384102
785124
4321x x x x ④ ⎥⎥
⎥
⎥
⎥
⎥
⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡----=⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣
⎡-555721
121
121
12
121
n n x x x x （n=5,10,100,…）
二、主要设备
笔记本 HP ProBook 6470b 一台 编译软件：VC++6.0 三、实验内容和原理
高斯消元法：
将原方程组化为三角形方阵的方程组：
l ik =a ik /a kk
a ij = a ij - l ik * a kj
（ k=1,2,…,n-1 i=k+1,k+2, …,n j=k+1,k+2, …,n+1 ）
由回代过程求得原方程组的解：
x n = a nn+1/ a nn
x k =( a kn+1-∑a kj x j )/ a kk
LU 分解法：
将系数矩阵A转化为A=L*U，L为单位下三角矩阵，U为普通上三角矩阵，然后通过解方程组l*y=b，u*x=y,来求解x。

四、操作方法与实验步骤
1. 完全主元素消元法：
#include<stdio.h>
#include<iostream.h>
#include"math.h"
float a[100][101];
float x[10];
int N;
void shuchu()
{
for(int i=1;i<=N;i++)
{
for(int j=1;j<=N+1;j++)
{
cout<<a[i][j]<<" "<<" ";
}
cout<<endl;
}
}
void shuru()
{
cout<<"请输入矩阵阶数:"<<endl;
cin>>N;
cout<<"请输入矩阵各项:"<<endl;
for(int i=1;i<=N;i++)
for(int j=1;j<=N+1;j++)
{
cin>>a[i][j];
}
cout<<endl;
}
void main()
{
int z[10];
int maxi,maxj;
shuru();
for(int i=1;i<=N;i++)
z[i]=i;
for(int k=1;k<N;k++)
{
maxi=k;maxj=k;float maxv=abs(a[k][k]);
for(i=k;i<=N;i++)
for(int j=k;j<=N;j++)
if(abs(a[i][j])>maxv)
{
maxv=abs(a[i][j]);maxi=i;maxj=j;
}
if(maxi!=k)
{
for(int j=1;j<=N+1;j++)
{
float t=a[k][j];a[k][j]=a[maxi][j];a[maxi][j]=t;
}
}
if(maxj!=k)
{
for(i=1;i<=N;i++)
{
float t=a[i][k];a[i][k]=a[i][maxj];a[i][maxj]=t;
}
int t=z[k];z[k]=z[maxj];z[maxj]=t;
}
for(int i=k+1;i<=N;i++)
{
float l=a[i][k]/a[k][k];
for(int j=k;j<=N+1;j++)
{
a[i][j]+=-l*a[k][j];
}
}
}
for(i=N;i>0;i--)
{
float s=0;
for(int j=i+1;j<=N;j++)
{
s+=a[i][j]*x[z[j]];
}
x[z[i]]=(a[i][N+1]-s)/a[i][i];
}
cout<<"完全主元素消去法之后的矩阵为：" <<endl;
shuchu();
for(i=1;i<=N;i++)
cout<<"x["<<i<<"]="<<x[i]<<endl;
}
2. 列主元素消元法：
#include"stdio.h"
int main()
{
float a[3][4]={{1,2,3,14},{0,1,2,8},{2,4,1,13}};
float x[3];
float sum=0;
int k,i,j;
for(k=0;k<2;k++)
for(i=k+1;i<3;i++)
for(j=k+1;j<4;j++)
a[i][j]=a[i][j]-a[i][k]/a[k][k]*a[k][j];
for(i=0;i<3;i++)
{
for(j=0;j<4;j++)
{
printf("a[%d][%d]=%f ",i,j,a[i][j]);
}
printf("\n");
}
x[2]=a[2][3]/a[2][2];
for(k=1;k>=0;k--)
{
sum=0;
for(j=k+1;j<3;j++)
{
sum+=a[k][j]*x[j];
}
x[k]=(a[k][3]-sum)/a[k][k];
}
for(i=0;i<3;i++)
{
printf ("x[%d]=%f\n",i+1,x[i]);
}
printf("\n");
}
3. LU分解法：
#include <stdio.h>
#include <math.h>
#define L 30
double a[L][L],b[L],l[L][L],u[L][L],x[L],y[L]; int main()
{
int n,i,j,k,r;
printf("请输入矩阵元次：\n");
scanf("%d",&n);
printf("请输入矩阵各项：\n");
for(i=1;i<=n;++i)
{
for(j=1;j<=n;++j)
{
scanf("%lf",&a[i][j]);
}
}
printf("请输入方程组的常数项：\n");
for(i=1;i<=n;++i)
{
scanf("%lf",&b[i]);
}
for(i=1;i<=n;++i)
{
for(j=1;j<=n;++j)
{
l[i][j]=0;
u[i][j]=0.0;
}
}
for(k=1;k<=n;++k)
{
for(j=k;j<=n;++j)
{
u[k][j]=a[k][j];
for(r=1;r<k;++r)
{
u[k][j]-=l[k][r]*u[r][j];
}
}
for(i=k+1;i<=n;++i)
{
l[i][k]=a[i][k];
for(r=1;r<k;++r)
{
l[i][k]-=l[i][r]*u[r][k];
}
l[i][k]/= u[k][k];
}
l[k][k]=1.0;
}
for(i=1;i<=n;++i)
{
y[i] = b[i];
for(j=1;j<i;++j)
{
y[i]-=l[i][j]*y[j];
}
}
for(i=n;i>0;--i)
{
x[i] = y[i];
for(j=i+1;j<=n;++j)
{
x[i]-=u[i][j]*x[j];
}
x[i]/= u[i][i];
}
for(i=1;i<=n;++i)
{
printf("%0.2lf\n",x[i]);
}
return 0;
}
五、实验结果与分析
完全主元素消元法：
列主元素消元法：
LU分解法：
分析：
对于两种高斯解方程，完全主元素跟列主元素都是先消元、再回代，由程序段可以发现，始终消去对角线下方的元素。

即，为了节约内存及时效，可以不必计算出主元素下方数据。

列主元素消元法的算法设计上优于完全主元素消元法，它只需依次按列选主元素然后换行使之变到主元素位置，再进行消元即可。

列主元素消元法的耗时比完全主元素法少很多，常采用之。

对于LU分解法，分解矩阵为单位下三角阵L与上三角阵U的乘积，然后解方程组Ly=b，回代，解方程组Ux=y。

其中的L为n阶单位下三角阵、U为上三角阵.
六、讨论、心得
本次试验中，感觉是最难的一次，完全主元素消元法程序编写过程相对来说花了好长时间。

纠正各种语法、算法、思路错误。

最后勉强成功，但还是有几处警告，不得解决之法。

感到程序学习的不足，再加之对高斯的不甚了解。

编写过程很是痛苦。

查阅各种内外部资料，这点有利有弊。

突然觉得，应该再把数据结构之类的重新学习一下才行。

以后多花时间在编程吧，重在理解。

必须反省一下自己的C、C++学习了，还是得多加练习，平时必须养成一种好的算法思维习惯。

实验地点虎峪校区致远楼B401指导教师xx
太原理工大学学生实验报告
学院名称软件学院专业班级1217班学号201200xxxx 学生姓名xx 实验日期2014.06.04 成绩
课程名称数值计算方法实验题目实验三线性方程组的迭代解法一、实验目的和要求
使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。

二、主要设备
笔记本 HP ProBook 6470b 一台
编译软件：VC++6.0
三、实验内容和原理
设线性方程组 A x=b
的系数矩阵A可逆，且主对角元素a11,a22,…,a nn均不为零，令
D=diag(a11,a22,…,a nn)
并将A分解成 A=(A-D)+D
从而线性方程组可写成 D x=(D-A)x+b
则有迭代公式
x(k+1)=B1x(k)+f1
其中，B1=I-D-1A,f1=D-1b。

四、操作方法与实验步骤
高斯—赛德尔迭代法
#include "iostream"
#include "iomanip"
using namespace std;
int main()
{
int i,j,k=0,m,n;
double t1,t2,e1,e2=0.0;
cout<<"请输入精度e：";
cin>>e1;
cout<<"请输入系数矩阵行数：";
cin>>m;
cout<<"请输入系数矩阵列数：";
cin>>n;
cout<<endl;
double (**a)=new double *[m];
for(i=0;i<=m;i++)
{
a[i]=new double[n];
}
double (*b)=new double [m];
double (*x)=new double [n];
cout<<"请输入系数矩阵："<<endl;
cout<<"-------------------------------------------------------------"<<endl;
for(int num1=0;num1<m;num1++)
{
for(int num2=0;num2<n;num2++)
{
cin>>a[num1][num2];
}
cout<<endl;
}
cout<<"输入的系数矩阵为："<<endl;
for (int num3=0;num3<m;num3++)
{
for(int num4=0;num4<n;num4++)
{
cout<<a[num3][num4]<<" ";
}
cout<<endl;
}
cout<<"请输入矩阵b："<<endl;
for(int num5=0;num5<m;num5++)
{
cin>>b[num5];
}
cout<<"输入的矩阵b为："<<endl;
for(int num6=0;num6<m;num6++)
{
cout<<b[num6]<<" ";
}
cout<<endl;
for(int num7=0;num7<n;num7++)
{
x[num7]=0.0000;
}
do
{
cout<<"第"<<k<<"次迭代值：";
e2=0.0;
for(i=0;i<m;i++)
{
double sum=0.0;
for(j=0;j<n;j++)
{
if(j!=i)
sum+=a[i][j]*x[j];
}
t1=x[i];
t2=e2;
x[i]=(b[i]-sum)/a[i][i];
e2=(x[i])-t1>=0?(x[i])-t1:t1-(x[i]);
e2=(e2>=t2?e2:t2);
cout<<setprecision(8)<<x[i]<<" ";
}
cout<<endl;
k++;
}while(e2>=e1&&k<30);
cout<<"共迭代了"<<k<<"次";
delete[]a;
delete[]b;
delete[]x;
return 0 ;
}
雅克比迭代法：
#include <stdio.h>
#include <math.h>
int main()
{
float a[3][3]={{10,-1,-2},{-1,10,-2},{-1,-1,5}},b[3]={7.2,8.3,4.2};
float x[3]={0,0,0},sum;
int i,j,k,n=3;
printf("\t\t X[1]\t\t X[2]\t\t X[3]\n");
for(k=0;k<8;k++)
{
for(i=0;i<3;i++)
{
sum=0;
for(j=0;j<n;j++)
{
if(i==j)continue;
sum=sum+a[i][j]*x[j];
}
x[i]=(b[i]-sum)/a[i][i];
}
printf("第%d次迭代：\t",k+1);
for(i=0;i<n;i++)
{
printf("%f\t",x[i]);
}
printf("\n");
}
}
五、实验结果与分析
高斯赛德尔迭代法：
雅克比迭代：
分析：
使用高斯-赛德尔和雅克比迭代都可以求出方程组的解，但是利用高斯-赛德尔迭代法所需的迭代次数比雅克比迭代少，能够更早的达到精度要求。

从程序中可以看出，雅克比定义的sum只有一个，而高斯赛德尔需要两个。

时效性上后者要好些。

六、讨论、心得
这次试验算是比较成功，要归功于授课时候的认真听讲。

程序编写之前，对书本的理论知识进行了进一步的探索。

预习准备工作很彻底，自然随后的一切也都很顺利。

实验地点虎峪校区致远楼B401指导教师xx
太原理工大学学生实验报告
学院名称 软件学院 专业班级 1217班 学号 201200xxxx
学生姓名 xx 实验日期 2014.06.11
成绩
课程名称
数值计算方法
实验题目
实验四 代数插值和最小二乘法拟合
一、实验目的和要求
(1)使用拉格朗日插值法或牛顿插值法求解：已知f(x)在6个点的函数值如下表所示，运用插值方法，求f(0.596)的近似值。

X 0.40 0.55 0.65 0.80 0.90 1.05 f(x)
0.41075
0.57815
0.69675
0.88811
1.02652
1.25386
(2)给定数据点（x i ,y i ），用最小二乘法拟合数据的多项式，并求平方误差。

二、 主要设备
笔记本 HP ProBook 6470b 一台 编译软件：VC++6.0
三、实验内容和原理
(1) 设函数在区间[a,b]上n+1互异节点
x 0,x 1,…,x n
上的函数值分别为
y 0,y 1,…,y n ,求n 次插值多项式P n
(x),满足条件
P n (x j )=y j , j=0,1,…,n
令
L n (x)=y 0l 0(x)+y 1l 1(x)+…+y n l n (x)= ∑y i l i (x)
其中l 0(x),l 1(x),…, l n (x) 为以x 0,x 1,…,x n 为节点的n 次插值基函数，
则Ln(x)是一次数不超过n 的多项式，且满足
L n (x j )=y j , L=0,1,…,n
再由插值多项式的唯一性，得
P n (x)≡L n (x)
(2 ) 建立正规方程组：
x i
0 0.5 0.6 0.7 0.8 0.9 1.0 y i
1
1.75
1.96
2.19
2.44
2.71
3.00
∑（∑x i j+k）a k=∑x i j y i ,j=0,1,…,n
平方误差：
I=∑（∑a k x i k-y i）2
对给定数据点{(Xi，Yi)}(i=0,1，…，m），在取定的函数类Φ 中，求p(x）∈Φ，使误差的平方和E2最小，E2=∑[p(Xi)-Yi]2。

从几何意义上讲，就是寻求与给定点{(Xi，Yi)}(i=0,1，…，m）的距离平方和为最小的曲线y=p(x）。

函数p(x）称为拟合函数或最小二乘解，求拟合函数p(x）的方法称为曲线拟合的最小二乘法。

得到的两个关于a0、a1为未知数的两个方程组，解这两个方程组得出：
a0 = （∑Yi) / m - a1（∑Xi) / m
a1 = [m∑Xi Yi - （∑Xi ∑Yi)] / [m∑Xi2 - （∑Xi)2 )]
即最终的拟合多项式各项系数。

四、操作方法与实验步骤
(1) 代数插值
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <alloc.h>
void difference(float *x,float *y,int n)
{
float *f;
int k,i;
f=(float *)malloc(n*sizeof(float));
for(k=1;k<=n;k++)
{
f[0]=y[k];
for(i=0;i<k;i++)f[i+1]=(f[i]-y[i])/(x[k]-x[i]);
y[k]=f[k];
}
return;
}
int main()
{
int i,n;
float x[20],y[20],xx,yy;
printf("请输入数据个数n：");
scanf("%d",&n);
printf("\n");
for(i=0;i<=n-1;i++)
{
printf("x[%d]=",i);
scanf("%f",&x[i]);
printf("y[%d]=",i);
scanf("%f",&y[i]);
printf("\n");
}
printf("\n");
difference(x,(float *)y,n);
printf("请输入插值X：");
scanf("%f",&xx);
yy=y[20];
for(i=n-1;i>=0;i--)
yy=yy*(xx-x[i])+y[i];
printf("\n近似值为：F(%f)=%f\n",xx,yy);
}
(2)最小二乘法拟合
#include<iostream.h>
#include<fstream.h>
#define N 15
double power(double &a,int n)
{
double b=1;
for(int i=0;i<n;i++)b*=a;
return b;
}
void Gauss();
double X[N],Y[N],sumX[N],sumY[N],a[N][N],b[N],l[N][N],x[N];
int main()
{
ofstream outdata;
ifstream indata;
double s;
int i,j,k,n,index;
cout<<"请输入已知点的个数n=";
cin>>n;
cout<<endl;
cout<<"请输入X和Y:"<<endl;
for(i=0;i<n;i++)
{
cout<<"X["<<i<<"]=";cin>>X[i];sumX[1]+=X[i];
cout<<"Y["<<i<<"]=";cin>>Y[i];sumY[1]+=Y[i];
cout<<endl;
}
cout<<"sumX[1]="<<sumX[1]<<"\t"<<"sumY[1]="<<sumY[1]<<endl;
cout<<"请输入拟合次数index=";cin>>index;
cout<<endl;
i=n;
sumX[0]=i;
for(i=2;i<=2*index;i++)
{
sumX[i]=0;
for(j=0;j<n;j++)sumX[i]+=power(X[j],i);
cout<<"sumX["<<i<<"]="<<sumX[i]<<endl;
}
for(i=2;i<=index+1;i++)
{
sumY[i]=0;
for(j=0;j<n;j++)sumY[i]+=power(X[j],i-1)*Y[j];
cout<<"sumY["<<i<<"]="<<sumY[i]<<endl;
}
for(i=1;i<=index+1;i++)
{
for(j=1;j<=index+1;j++)a[i][j]=sumX[i+j-2];b[i]=sumY[i];
}
k=1;
do
{
for(j=k+1;j<=index+1;j++) l[j][k]=a[j][k]/a[k][k];
for(i=k+1;i<=index+1;i++)
{
for(j=k+1;j<=index+1;j++)
a[i][j]=a[i][j]-l[i][k]*a[k][j];
b[i]=b[i]-l[i][k]*b[k];
}
if(k==index+1)break;
k++;
}while(1);
x[index+1]=b[index+1]/a[index+1][index+1];
for(i=index;i>=1;i--)
{
s=0;
for(j=i+1;j<=index+1;j++)s=s+a[i][j]*x[j]; x[i]=(b[i]-s)/a[i][i]; }
cout<<"拟合系数为:";
for(i=1;i<=index+1;i++)cout<<x[i]<<"\t";
double m=0;
cout<<endl<<"平方误差为:";
for(i=0;i<n;i++)
{
double t=x[1]+x[2]*X[i]-Y[i];
m=m+power(t,2);
}
cout<<m<<endl;
}
五、实验结果与分析
(1) 代数插值
分析：
拉格朗日插值的优点是插值多项式特别容易建立，缺点是增加节点是原有多项式不能利用，必须重新建立，即所有基函数都要重新计算，这就造成计算量的增加。

牛顿插值法则很好地避免了上述问题。

(2)最小二乘法拟合
分析：
数据拟合的具体作法是：对给定的数据（x i ,y i）（i=0,1,…,m），在取定的函
数类中，求p(x)属于此函数类，使误差r
i =p(x
i
)- y
i
(i=0,1,…,m)的平方和最小，
即∑r i2=∑（∑p(x i)-y i）2=min
从几何意义上讲，就是寻求与给定点（x i,y i）（i=0,1,…,m）的距离平方和为最小的曲线y=p(x)。

六、讨论、心得
本实验中，插值法有两种插值方法可以选用，最终选用了牛顿插值法。

最小二乘法拟合中，很好地实现了最小二乘法的程序模拟。

对程序编写又更进一步，分析算法、尝试编写。

检查并修改其中错误。

感觉收获很大。

特别是对平方误差的计算模块，一次成功。

最后还对程序进行优化，删去冗余以及标准化模块。

实验地点虎峪校区致远楼B401指导教师xx
word文档可自由复制编辑。