《计算机数值方法C》实验指导书新.docx

《计算机数值方法C》实验指导书新.docx

《《计算机数值方法C》实验指导书新.docx》由会员分享，可在线阅读，更多相关《《计算机数值方法C》实验指导书新.docx（21页珍藏版）》请在冰点文库上搜索。

《计算机数值方法C》实验指导书新

《计算机数值方法》

实验指导书

（2017版）

太原理工大学

计算机科学与技术学院

二〇一七年三月

目录

实验一方程求根2

实验二线性方程组的直接解法4

实验三线性方程组的迭代解法9

实验四代数插值10

实验五最小二乘法拟合多项式11

实验一方程求根

【目的与要求】

1．了解方程求根的基本方法、基本原理、误差估计；

2．能够按照工程实际要求，选择适当的算法；

3．通过编写程序，进行算法设计和数值求解。

【实验内容】

1.熟悉使用二分法、迭代法、牛顿法、割线法等方程求根的基本方法、基本原理、误差估计的相关理论。

2.选择方程求解方法中的两种方法求方程：

f（x）=x3+4x2-10=0在[1,2]内的一个实根，且要求满足精度|x*-xn|<0.5×10-5。

【示例程序】

增值寻根法又叫做逐步搜索法，是求解方程根的一种方法，也是初步确定方程隔根区间的一种方法。

增值寻根法的基本思想是：

从初值x0开始,按规定的一个初始步长h来增值：

xn+1=xn+hn=0,1,2,3……

同时计算f（xn+1）可能遇到三种情形：

（1）f（xn+1）=0，此时，xn+1即为方程的根x*；

（2）f（xn+1）*f（xn）>0，说明[xn+1,xn]内无根；

（3）f（xn+1）*f（xn）<0，说明[xn+1,xn]内有根；

搜索过程,可从a开始,也可从b开始，这时应取更小的步长h，直到有根区间的长度|xn+1-xn|

（e为所要求的精度），此时f（xn）或f（xn+1）就可近似认为是零，xn或xn+1就是满足精度的方程的近似根。

基于增值寻根法，设计程序如下：

#include

doublefun（doublex）{

doubley=x*x;

returny*x+4*y-10;

}

voidmain（void）{

doublex0=-4,x;

doubleh=1;

doublestep=0;

x=x0;

while（h>0.000001）{

cout<<"-------h="<

while（true）{

step++;

if（fun（x）==0）{

cout<<"x="<

h=0;

break;

}

elseif（fun（x）*fun（x+h）>0）{

x+=h;

continue;

}elseif（fun（x）*fun（x+h）<0）{

cout<<"x1="<

cout<<"x2="<<（x+h）<

break;

}

}

h=h/10;

}

cout<<"step="<

}

实验二线性方程组的直接解法

【目的与要求】

1．了解Gauss消元法、LU分解法等线性方程组直接求解的基本方法、基本原理；

2．能够按照工程实际要求，选择适当的算法；

3．通过编写程序，进行算法设计和数值求解。

【实验内容】

合理利用高斯消元法、高斯主元素消元法、LU分解法求解下列方程组：

①

②

③

【示例程序】

高斯消元法的基本思想是将原有线性方程组化为与之等价的三角形方程组，这个过程也称为消元过程，再通回代过程进行求解。

消元过程使用公式为：

回代过程使用公式如下：

首先，设定数据文件Gaussdata.txt中的数据为要求解的线性方程组数据，如下：

3

2233

4771

-245-7

其中的3表示三阶线性方程组，下面是三行四列的系数矩阵。

编写函数displayA（）以显示系数矩阵及其的运算中的变化如下：

voiddisplayA（）

{

printf（"\n"）;

for（intj=1;j<=n;j++）

{

for（inti=1;i<=n+1;i++）

printf（"%lf",a[j][i]）;

printf（"\n"）;

}

}

定义全局变量如下：

doublea[15][16],a0[15][16];//a和a0用以记录方程组系数，

//其中：

a为本次运算的结果，a0为上次运算的结果

doublel[15],tmp;//l用做消元过程使用系数l的计算和存储，tmp为零时变量

intn;//n为线性方程组阶数

完整的程序如下：

#include"stdafx.h"

#include"stdio.h"

doublea[15][16],a0[15][16];

doublel[15],tmp;

intn;

voiddisplayA（）

{

printf（"\n"）;

for（intj=1;j<=n;j++）

{

for（inti=1;i<=n+1;i++）

printf（"%lf",a[j][i]）;

printf（"\n"）;

}

}

voidmain（）

{

FILE*f;

inti,j,k;

f=fopen（"Gaussdata.txt","r"）;

fscanf（f,"%d",&n）;

for（i=1;i<=n;i++）

{

for（j=1;j<=n+1;j++）

{

fscanf（f,"%lf",&a[i][j]）;

a0[i][j]=a[i][j];

}

}

displayA（）;

fclose（f）;

//消元过程

k=1;

do

{

for（i=k+1;i<=n;i++）

{

l[i]=a0[i][k]/a0[k][k];

printf（"l[%i][%i]=%lf",i,k,l[i]）;

for（j=k+1;j<=n+1;j++）

{

a[i][j]=a0[i][j]-l[i]*a0[k][j];

}

displayA（）;

}

k++;

if（k==n）break;

for（j=1;j<=n;j++）

{

for（i=1;i<=n+1;i++）

a0[j][i]=a[j][i];

}

}while

（1）;

//回代过程

l[n]=a[n][n+1]/a[n][n];

for（k=n-1;k>=1;k--）

{

tmp=0;

for（j=k+1;j<=n;j++）tmp+=a[k][j]*l[j];

l[k]=（a[k][n+1]-tmp）/a[k][k];

}

for（i=1;i<=n;i++）printf（"x[%i]=%lf\n",i,l[i]）;

}

运行结果如下：

2.0000002.0000003.0000003.000000

4.0000007.0000007.0000001.000000

-2.0000004.0000005.000000-7.000000

l[2][1]=2.000000

2.0000002.0000003.0000003.000000

0.0000003.0000001.000000-5.000000

-2.0000004.0000005.000000-7.000000

l[3][1]=-1.000000

2.0000002.0000003.0000003.000000

0.0000003.0000001.000000-5.000000

0.0000006.0000008.000000-4.000000

l[3][2]=2.000000

2.0000002.0000003.0000003.000000

0.0000003.0000001.000000-5.000000

0.0000000.0000006.0000006.000000

x[1]=2.000000

x[2]=-2.000000

x[3]=1.000000

实验三线性方程组的迭代解法

【目的与要求】

1．了解雅可比迭代法、高斯-赛德尔迭代法等线性方程组迭代求解的基本方法、基本原理；

2．能够按照工程实际要求，选择适当的算法；

3．通过编写程序，进行算法设计和数值求解。

【实验内容】

使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。

【内容提示】

雅可比迭代法的公式如下：

高斯-赛德尔迭代法的公式如下：

实验五代数插值

【目的与要求】

1．了解拉格朗日插值法或牛顿插值法的基本方法、基本原理；

2．通过编写程序，进行算法设计和数值求解。

【实验内容】

使用拉格朗日插值法或牛顿插值法求解：

已知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

【内容提示】

拉格朗日基函数为：

拉格朗日插值多项式为：

实验六最小二乘法拟合多项式

【目的与要求】

1．了解最小二乘法拟合多项式的基本方法、基本原理；

2．通过编写程序，进行算法设计和数值求解。

【实验内容】

给定数据点（xi,yi），用最小二乘法拟合数据的多项式，并求平方误差。

xi

0

0.5

0.6

0.7

0.8

0.9

1.0

yi

1

1.75

1.96

2.19

2.44

2.71

3.00

【示例程序】

定义要拟合的数据文件数据如下：

9

110

35

44

52

61

71

82

93

104

2

一共9个点，最后的2是要拟合的次数。

程序使用完全主元素消去法对方程组进行了求解，得到2次拟合函数。

源程序如下：

#include"stdio.h"

#include"math.h"

#definenum10

doublex[num],y[num],sumX[num],sumY[num];

doublea[num][num],a0[num][num];

doubleneiji（doubleb[],doublec[]）/*内积函数*/

{

intp;

doublenj=0;

for（p=1;p

nj+=c[p]*b[p];

returnnj;

}

doublepower（double&a,intn）

{

doubleb=1;

for（inti=0;i

{

b*=a;

}

returnb;

}

//完全主元素消去法相关程序

voiddisplayA（doublea[][num],intn）

{

printf（"\n"）;

for（intj=1;j<=n+1;j++）

{

for（inti=1;i<=n+1;i++）

printf（"%lf",a[j][i]）;

printf（"\n"）;

}

}

voidmain（）

{

inti,j,n,k,index;

doublel[15],x[15],sum,max;

intmaxL,maxC;

FILE*f;

f=fopen（"zxec.txt","r"）;

fscanf（f,"%d",&n）;

printf（"已知点的个数n=%i\n",n）;

sumX[1]=sumY[1]=0;

printf（"No.xy\n"）;

for（i=1;i<=n;i++）

{

fscanf（f,"%lf",&x[i]）;

printf（"%i%lf",i,x[i]）;

sumX[1]+=x[i];

fscanf（f,"%lf",&y[i]）;

printf（"%lf\n",y[i]）;

sumY[1]+=y[i];

}

printf（"\nsumX[1]=%lf\nsumY[1]=%lf",sumX[1],sumY[1]）;

fscanf（f,"%d",&index）;

printf（"\n拟和次数=%i\n",index）;

fclose（f）;

i=n;

sumX[0]=i;

for（i=2;i<=2*index;i++）

{

for（j=1;j<=n;j++）

{

sumX[i]+=power（x[j],i）;

}

printf（"sumX[%d]=%lf\n",i,sumX[i]）;

}

for（i=2;i<=index+1;i++）

{

for（j=1;j<=n;j++）

{

sumY[i]+=power（x[j],i-1）*y[j];

}

printf（"sumY[%d]=%lf\n",i,sumY[i]）;

}

f=fopen（"Gaussdata.txt","w"）;

fprintf（f,"%i\n",index+1）;

printf（"\n\n%i\n",index+1）;

for（i=0;i<=index;i++）

{

for（j=i;j<=index+i;j++）

{

fprintf（f,"%lf",sumX[j]）;

printf（"%lf",sumX[j]）;

}

fprintf（f,"%lf\n",sumY[i+1]）;

printf（"%lf\n",sumY[i+1]）;

}

fclose（f）;

//完全主元素消去法相关程序

f=fopen（"Gaussdata.txt","r"）;

fscanf（f,"%d",&n）;

for（i=1;i<=n;i++）

{

for（j=1;j<=n+1;j++）

{

fscanf（f,"%lf",&a[i][j]）;

a0[i][j]=a[i][j];

}

}

for（j=1;j<=n+1;j++）

{

a[n+1][j]=j;

}

displayA（a,n）;

fclose（f）;

//消元过程

k=1;

do

{

//选择主元素

max=0;maxL=k;maxC=k;

for（i=k;i<=n;i++）

{

for（j=k;j<=n;j++）

if（fabs（max）

{

max=a[i][j];

maxL=i;maxC=j;

}

}

printf（"（%i,%i）<-->（%i,%i）",k,k,maxL,maxC）;

//进行行变换

if（maxL!

=k）

{

for（i=1;i<=n+1;i++）

{

sum=a[maxL][i];

a[maxL][i]=a[k][i];

a[k][i]=sum;

a0[maxL][i]=a0[k][i];

a0[k][i]=sum;

}

}

//进行列变换

if（maxC!

=k）

{

for（i=1;i<=n+1;i++）

{

sum=a[i][k];

a[i][k]=a[i][maxC];

a[i][maxC]=sum;

a0[i][k]=a0[i][maxC];

a0[i][maxC]=sum;

}

}

displayA（a,n）;

//进行消元

for（i=k+1;i<=n;i++）

{

l[i]=a0[i][k]/a0[k][k];

printf（"l[%i][%i]=%lf",i,k,l[i]）;

//for（j=k+1;j<=n+1;j++）//书上的公式

for（j=1;j<=n+1;j++）//出题用

{

a[i][j]=a0[i][j]-l[i]*a0[k][j];

}

displayA（a,n）;

}

k++;

if（k==n）break;

for（j=1;j<=n;j++）

{

for（i=1;i<=n+1;i++）

a0[j][i]=a[j][i];

}

}while

（1）;

//回代过程

l[n]=a[n][n+1]/a[n][n];

for（k=n-1;k>=1;k--）

{

sum=0;

for（j=k+1;j<=n;j++）sum+=a[k][j]*l[j];

l[k]=（a[k][n+1]-sum）/a[k][k];

}

for（j=1;j<=n;j++）

for（i=1;i<=n;i++）

{

if（j!

=a[n+1][i]）continue;

x[j]=l[i];

printf（"x[%i]=%lf\n",j,x[j]）;

break;

}

printf（"y=%lf",x[1]）;

for（i=2;i<=n;i++）

printf（x[i]>=0?

"+%lfx^%i":

"%lfx^%i",x[i],i-1）;

printf（"\n"）;

}

程序运行结果如下：

已知点的个数n=9

No.xy

11.00000010.000000

23.0000005.000000

34.0000004.000000

45.0000002.000000

56.0000001.000000

67.0000001.000000

78.0000002.000000

89.0000003.000000

910.0000004.000000

sumX[1]=53.000000

sumY[1]=32.000000

拟和次数=2

sumX[2]=381.000000

sumX[3]=3017.000000

sumX[4]=25317.000000

sumY[2]=147.000000

sumY[3]=1025.000000

3

9.00000053.000000381.00000032.000000

53.000000381.0000003017.000000147.000000

381.0000003017.00000025317.0000001025.000000

9.00000053.000000381.00000032.000000

53.000000381.0000003017.000000147.000000

381.0000003017.00000025317.0000001025.000000

1.0000002.0000003.0000004.000000

（1,1）<-->（3,3）

25317.0000003017.000000381.0000001025.000000

3017.000000381.00000053.000000147.000000

381.00000053.0000009.00000032.000000

3.0000002.0000001.0000004.000000

l[2][1]=0.119169

25317.0000003017.000000381.0000001025.000000

0.00000021.4673147.59663524.851839

381.00000053.0000009.00000032.000000

3.0000002.0000001.0000004.000000

l[3][1]=0.015049

25317.0000003017.000000381.0000001025.000000

0.00000021.4673147.59663524.851839

0.0000007.5966353.26626416.574594

3.0000002.0000001.0000004.000000

（2,2）<-->（2,2）

25317.0000003017.000000381.0000001025.000000

0.00000021.4673147.59663524.851839

0.0000007.5966353.26626416.574594

3.0000002.0000001.0000004.000000

l[3][2]=0.353870

25317.0000003017.000000381.0000001025.000000

0.00000021.4673147.59663524.851839

0.0000000.0000000.5780447.780278

3.0000002.0000001.0000004.000000

x[1]=13.459664

x[2]=-3.605309

x[3]=0.267571

y=13.459664-3.605309x^1+0.267571x^2

提醒：

每个实验结束后，应整理出实验报告，实验报告模版到学校教务处网站下载。