计算方法实验报告汇总

中 南 林 业 科 技 大 学

计算机与信息工程学院

《数值分析》

实 验 报 告

      姓名：          尹一君                  __

学号：          20134651                  

指导老师：      赵红敏                    

专业班级：      计算机科学与技术2班      

指导教师评语：                              成绩：

签名：           

实验二 非线性方程求根    - 6 -

实验三 线性方程组的直接解法    - 13 -

实验四 解线性方程组的迭代法    - 21 -

实验五 函数插值方法    - 27 -

实验六 函数逼近与曲线拟合    - 31 -

实验七 数值积分与数值微分    - 36 -

参考文献    - 42 -

实验一 算法设计基础

1．问题提出

  编程实现以下矩阵相乘：

2．要求

   1.编织一个程序进行运算，最后打印出结果；

   2.将编程结果与手工演算结果进行比较；

   3.原始数据使用ASCII码数据文件的形式访问；

   4.最好能将程序结果输出到ASCII码数据文件。

3．实验目的和意义

   1.通过实验进一步熟悉矩阵运算和多重循环编程；

   2.熟悉数据文件访问操作。

四．计算公式

    在计算机中，一个矩阵实际上就是一个二维数组。一个m行n列的矩阵与一个n行p列的矩阵可以相乘，得到的结果是一个m行p列的矩阵，其中的第i行第j列位置上的数为第一个矩阵第i行上的n个数与第二个矩阵第j列上的n个数对应相乘后所得的n个乘积之和。

五．结构程序设计

1.

#define M 4

#define N 3

#define L 3

#include"stdio.h"

void main()

{

    FILE *f;

    float A[M][L],B[L][N],C[M][N];

    int i,j,k;

    if((f=fopen("Array1.txt

        return;

for(i=0;i        for(j=0;j            fscanf(f,"%f",&A[i][j]);

for(i=0;i        for(j=0;j            fscanf(f,"%f",&B[i][j]);

for(i=0;i    {

        for(j=0;j        {

            C[i][j]=0;

            for(k=0;k                C[i][j]=C[i][j]+A[i][k]*B[k][j];

            printf("%10.2f",C[i][j]);

        }

        printf("\\n");

    }

}

2.

#define M 2

#define N 3

#define L 4

#include"stdio.h"

void main()

{

    FILE *f1;

    int A[M][L],B[L][N],C[M][N];

    int i,j,k;

    if((f1=fopen("Array2.txt

        return;

for(i=0;i        for(j=0;j            fscanf(f1,"%d",&A[i][j]);

for(i=0;i        for(j=0;j            fscanf(f1,"%d",&B[i][j]);

for(i=0;i    {

        for(j=0;j        {

            C[i][j]=0;

            for(k=0;k                C[i][j]=C[i][j]+A[i][k]*B[k][j];

            printf("%4d",C[i][j]);

        }

        printf("\\n");

    }

}

六．结果展示

  程序1和2的运行结果分别如图1-1和1-2所示：

          图1-1                                 图1-2

七．实验讨论与总结

1.本实验考察了对数据文件的建立、读入和试用，是一个突破；

  2.第一个矩阵相乘中出现小数，故应用float数据类型；第二个矩阵中全是整数，故用int型即可。

  3.对于输出结果的排版需要格外注意，应灵活使用输出格式符使其美观完整。

  4.文件操作的基本函数（打开、关闭、输入和输出）

    ①fopen（打开文件）

    相关函数 open，fclose

    表头文件 #include

    定义函数 FILE * fopen(const char * path,const char * mode);

    函数说明 参数path字符串包含欲打开的文件路径及文件名，参数mode字符串则代表着流形态。

    ②fclose（关闭文件）

    相关函数 close，fflush，fopen，setbuf

    表头文件 #include

    定义函数 int fclose(FILE * stream);

    函数说明 fclose()用来关闭先前fopen()打开的文件。此动作会让缓冲区内的数据写入文件中，并释放系统所提供的文件资源。

    返回值 若关文件动作成功则返回0，有错误发生时则返回EOF并把错误代码存到errno。

错误代码 EBADF表示参数stream非已打开的文件。

    ③fscanf（输入函数）

    功能：从一个流中执行格式化输入

    表头文件：#include

    函数原型：int fscanf(FILE *stream, char *format[,argument...]);

    FILE* 一个FILE型的指针

    char* 格式化输出函数，和scanf里的格式一样

    返回值：成功时返回转换的字节数，失败时返回一个负数

    fp = fopen("/local/test.c

fscanf(fp,"%s",str);

    ④fprintf

    功能：传送格式化输出到一个文件中

    表头文件：#include

    函数原型：int fprintf(FILE *stream, char *format[, argument,...]);

    FILE* 一个FILE型的指针

    char* 格式化输入函数，和printf里的格式一样

    返回值：成功时返回转换的字节数，失败时返回一个负数

    fp = fopen("/local/test.c

    fprintf(fp,"%s\\n",str)

实验二 非线性方程求根

1．问题提出

二．要求

  1.编制一个程序进行运算，最后打印出每种迭代格式的敛散情况；

  2.用事后误差估计|xk+1-xk|<ε来控制迭代次数，并且打印出迭代的次数；

  3.初始值的选取对迭代收敛有何影响；

  4.分析迭代收敛和发散的原因。

三．实验目的和意义

  1.通过实验进一步了解方程求根的算法；

  2.认识选择计算格式的重要性；

  3.掌握迭代算法和精度控制；

  4.明确迭代收敛性与初值选取的关系。

四．计算公式

    对于给定的线性方程组x=Bx+f（这里的xBf同为矩阵，任意线性方程组都可以变换成此形式），用公式x(k+1）=Bx(k)+f（括号中为上标，代表迭代k次得到的x，初始时k=0)逐步带入求近似解的方法称为迭代法（或称一阶定常迭代法）。如果k趋向无穷大时limx（k）存在，记为x*，称此迭代法收敛。显然x*就是此方程组的解，否则称为迭代法发散。

五．结构程序设计

各函数及其导数建立如图2-1和图2-2所示：

  图2-1                      图2-2

主函数如下：

void main()

{

    double i,j;

    double eps=1e-10;

    int k;

if(fabs(f11(1))>1 || fabs(f11(2))>1)

            printf("不满足收敛条件\\n");

            else

            {

               k=0;

               i=f1(1);j=i; i=f1(j);

while(fabs(i-j)>eps)

               {

                  j=i;

                  i=f1(j);k++;

               }

               printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

            }

if(fabs(f22(1))>1 || fabs(f22(2))>1)

            printf("不满足收敛条件\\n");

            else

            {

               k=0;

               i=f2(1);j=i; i=f2(j);

while(fabs(i-j)>eps)

               {

                  j=i;

                  i=f2(j);k++;

               }

               printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

            }

if(fabs(f33(1))>1 || fabs(f33(2))>1)

            printf("不满足收敛条件\\n");

            else

            {

               k=0;

               i=f3(1);j=i; i=f3(j);

while(fabs(i-j)>eps)

               {

                  j=i;

                  i=f3(j);k++;

               }

               printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

            }

if(fabs(f44(1))>1 || fabs(f44(2))>1)

            printf("不满足收敛条件\\n");

            else

            {

               k=0;

               i=f4(1);j=i; i=f4(j);

while(fabs(i-j)>eps)

               {

                  j=i;

                  i=f4(j);k++;

               }

               printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

            }

if(fabs(f55(1))>1 || fabs(f55(2))>1)

            printf("不满足收敛条件\\n");

            else

            {

               k=0;

               i=f5(1);j=i; i=f5(j);

while(fabs(i-j)>eps)

               {

                  j=i;

                  i=f5(j);k++;

               }

               printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

            }

if(fabs(f66(1))>1 || fabs(f66(2))>1)

            printf("不满足收敛条件\\n");

            else

            {

               k=0;

               i=f6(1);j=i; i=f6(j);

while(fabs(i-j)>eps)

               {

                  j=i;

                  i=f6(j);k++;

               }

               printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

            }

}

五．结果展示

程序运行结果如图2-3所示：

图2-3

六．实验讨论与总结

1.完成本实验相当于实现了将迭代算法转化成程序语言这一突破；

2.本实验考察众多容易忽视的细节和学生的仔细程度，比如六个函数转化成C语言时，语法会出现众多不适，括号较多，很容易出错；另一方面，在用pow函数求x的y次方时，容易将x的1/3次方写成pow（x，1/3），这样写是错误的。因为1/3为两int型整数相除，系统执行结果为0，故应写成pow（x，1.0/3）。这是一个极容易忽视的地方；

☆  3.方法改进

由于求导法判定迭代收敛时，必定在一个区间内，而这个区间内只能存在一个根。因此用该方法判断出来的迭代结果如果是不收敛，并不能说明它对于其余的根也不收敛。因此，求导法只能求出f(x)的一个根，不能满足题目要求，而且求导需人工执行、过程繁琐。故对该算法进行改进：

#include"math.h"

#include"stdio.h"

double f1(double x)

{ return ((3*x+1)/(x*x));}

    double f2(double x)

{ return((pow(x,3)-1)/3);}

double f3(double x)

{ return(pow(3*x+1,1.0/3));}

    double f4(double x)

{ return(1/(x*x-3));}

    double f5(double x)

{ return(pow((3+1/x),1.0/2));}

    double f6(double x)

{ return(x-(x*x*x-3*x-1)/(3*(x*x-1)));}       //无需建立导函数

void main()

{

    double i,j;

    double eps=1e-5;

    int k;

    k=0;

    i=f1(1);j=i; i=f1(j);

while(fabs(i-j)>eps )

    {

        j=i;

        i=f1(j);k++;

    }

    if(i<-99999) printf("不满足收敛条件\\n");

    else printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

    k=0;

    i=f2(1);j=i; i=f2(j);

while(fabs(i-j)>eps)

    {

        j=i;

        i=f2(j);k++;

    }

    if(i<-99999) printf("不满足收敛条件\\n");

    else printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

    k=0;

    i=f3(1);j=i; i=f3(j);

while(fabs(i-j)>eps)

    {

        j=i;

        i=f3(j);k++;

    }

    if(i<-99999) printf("不满足收敛条件\\n");

    else printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

    k=0;

    i=f4(1);j=i; i=f4(j);

while(fabs(i-j)>eps)

    {

        j=i;

        i=f4(j);k++;

    }

    if(i<-99999) printf("不满足收敛条件\\n");

    else printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

    k=0;

    i=f5(1);j=i; i=f5(j);

while(fabs(i-j)>eps)

    {

        j=i;

        i=f5(j);k++;

    }

    if(i<-99999) printf("不满足收敛条件\\n");

    else printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

    k=0;

    i=f6(1);j=i; i=f6(j);

while(fabs(i-j)>eps)

    {

        j=i;

        i=f6(j);k++;

    }

    if(i<-99999) printf("不满足收敛条件\\n");

    else printf("满足收敛条件，x≈%f,迭代次数为%d次\\n",i,k);

}

改进算法的程序运行结果如图2-4所示：

                   图2-4

    该方法的思想是：无区间迭代，如果迭代结束得到的结果为无穷大或无穷小，则说明该迭代不收敛。这种做法更加实用，而且可以看出，求得了两个根。但就科学性而言，教材上定理2.1的内容显示了迭代法收敛的条件。定义明确指出用求导判断收敛性，根据定义判断才是最科学的方法。但考虑到计算机语言的便利性就在于全过程无人工计算，而手动求导的方法与之相违背，从这个角度考虑，改进法更加具有说服力。但若只需求一个根而又要讲究科学性，我们仍然可以选择求导法。因此，我们应做到根据情况的不同来选择不同的适当方法来解决问题。 

    4.初始值的选取对迭代收敛有何影响？

初始值的选取虽然不会改变迭代收敛的性质，但是会影响其收敛速度。为了提高收敛速度，可设法提高初值的精度以减少迭代的次数。

5.分析迭代收敛和发散的原因。

迭代公式收敛与否，完全取决于迭代公式ψ(x)的性态。这里考察到迭代法的几何意义：方程x=ψ(x)的求根问题在几何上就是确定曲线y=ψ(x)与直线y=x的交点P的横坐标。按教材图2.3演示的路径走下去，在曲线y=ψ(x)上得到点列P0,P1,…,其横坐标分别按迭代公式所确定的迭代值x1,x2,…,如果迭代收敛，则点P1,P2,…将越来越逼近所求的交点P。否则迭代法发散。这是迭代收敛和发散的本质原因。

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

1．问题提出

    给出下列几个不同类型的线性方程组，请用适当算法计算其解。

二．要求

  1.对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改方根法；追赶法求解（选择其一）；

  2.应用结构程序设计编出通用程序；

  3.比较计算结果，分析数值解误差的原因；

  4.尽可能利用相应模块输出系数矩阵的三角分解式。

三．实验目的和意义

  1.通过该课题的实验，体会模块化结构程序设计方法的优点；

  2.运用所学的计算方法，解决各类线性方程组的直接算法；

  3.提高分析和解决问题的能力，做到学以致用；

  4.通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。

四．计算公式

    本次实验的三个方程组，在此分别使用Gauss顺序消去法、平方根法和追赶法求解。

下面给出这三种方法的计算公式。

    1.Gauss顺序消去法计算公式

   2.平方根法计算公式

   3.追赶法计算公式

五．结构程序设计

1.用Gauss顺序消去法求解线性方程组的程序设计如下（问题1）：

//Gauss顺序消去法--求解线性方程组

//设计人：尹一君  20134651

#include"stdio.h"

#include"math.h"

void main()

{

    FILE *f;

    double a[20][20],b[20],x[20];

    double m,s;

    int i,j,k;

    if((f=fopen("Array3.txt

        return;

for(i=1;i<=10;i++)

        for(j=1;j<=10;j++)

            fscanf(f,"%lf",&a[i][j]);

for(i=1;i<=10;i++)

        fscanf(f,"%lf",&b[i]);

    fclose(f);

for(k=1;k<10;k++)

        for(i=k+1;i<=10;i++)

        {

            for(j=k+1;j<=10;j++)

            {

                m=a[i][k]/a[k][k];

                a[i][j]=a[i][j]-m*a[k][j];

            }

            b[i]=b[i]-m*b[k];

        }

    x[10]=b[10]/a[10][10];

for(i=9;i>0;i--)

    {

        s=0;

        for(j=i+1;j<=10;j++)

            s=s+a[i][j]*x[j];

        x[i]=(b[i]-s)/a[i][i];

    }

for(i=1;i<=10;i++)

    {

        if(fabs(x[i])<0.5)

            x[i]=fabs(x[i]);

        printf("x[%d]=%f\\n",i,x[i]);

    }

}

2.用平方根法求解对称正定阵系数阵线方程组的程序设计如下（问题2）：

//平方根法--求解对称正定阵系数阵线方程组

//设计人：尹一君  20134651

#include"stdio.h"

#include"math.h"

void main()

{

    FILE *f;

    double a[10][10],b[10],l[10][10],x[10],y[10];

    double s1,s2,s3,s4;

    int i,j,k;

    if((f=fopen("Array4.txt

        return;

for(i=1;i<=8;i++)

        for(j=1;j<=8;j++)

            fscanf(f,"%lf",&a[i][j]);

for(i=1;i<=8;i++)

        fscanf(f,"%lf",&b[i]);

    fclose(f);

    l[1][1]=pow(a[1][1],1.0/2);  

for(i=1,j=2;j<=8;j++)

        l[j][i]=a[j][i]/l[1][1];

for(i=2;i<=8;i++)

for(j=i;j<=8;j++)

           if(i==j)

           {

               s1=0;

               for(k=1;k<=i-1;k++)

                   s1=s1+l[i][k]*l[i][k];

               l[i][i]=pow((a[i][i]-s1),1.0/2);

           }  

           else

           {

               s2=0;

               for(k=1;k<=i-1;k++)

                   s2=s2+l[j][k]*l[i][k];

               l[j][i]=(a[j][i]-s2)/l[i][i];

           }

     y[1]=b[1]/l[1][1];

for(i=2;i<=8;i++)

     {

         s3=0;

         for(k=1;k<=i-1;k++)

             s3=s3+l[i][k]*y[k];

         y[i]=(b[i]-s3)/l[i][i];

     }

         x[8]=y[8]/l[8][8];

for(i=7;i>=1;i--)

     {

         s4=0;

         for(k=i+1;k<=8;k++)

             s4=s4+l[k][i]*x[k];

         x[i]=(y[i]-s4)/l[i][i];

     }

for(i=1;i<=8;i++)

    {

        if(fabs(x[i])<0.5)

            x[i]=fabs(x[i]);

        printf("x[%d]=%f\\n",i,x[i]);

    }

}

3.用追赶法求解三对角形线性方程组的程序设计如下（问题3）：

//追赶法--求解三对角形线性方程组

//设计人：尹一君  20134651

#include"stdio.h"

#include"math.h"

void main()

{

    FILE *p;

    double a[15],b[15],c[15],f[15],l[15],u[15],x[15],y[15];

    int i;

    if((p=fopen("Array5.txt

        return;

for(i=1;i<=10;i++)

        fscanf(p,"%lf",&a[i]);

for(i=1;i<=10;i++)

        fscanf(p,"%lf",&b[i]);

for(i=1;i<=10;i++)

        fscanf(p,"%lf",&c[i]);

for(i=1;i<=10;i++)

        fscanf(p,"%lf",&f[i]);

    fclose(p);

for(i=1;i<=9;i++)

    {

        if(i==1)  l[i]=b[i];

        u[i]=c[i]/l[i];

        l[i+1]=b[i+1]-a[i+1]*u[i];

    }

for(i=1;i<=10;i++)

        if(i==1)  

            y[i]=f[i]/l[i];

        else

            y[i]=(f[i]-a[i]*y[i-1])/l[i];

for(i=10;i>=1;i--)

    {

        if(i==10)

            x[i]=y[i];

        else

            x[i]=y[i]-u[i]*x[i+1];

    }

for(i=1;i<=10;i++)

    {

        if(fabs(x[i])<0.5)

            x[i]=fabs(x[i]);

        printf("x[%d]=%f\\n",i,x[i]);

    }

}

六．结果展示

   问题1、问题2和问题3的程序运行结果分别如图3-1、3-2和3-3所示：

图3-1                      图3-2                      图3-3

七．实验讨论与总结

1.一般线性方程组使用高斯消去法求解时，在消元过程中可能会出现akk=0的情况，这时消去法将无法进行；即使akk≠0，但它的绝对值很小，用其作除数，会导致其他元素数量级的严重增长和舍入误差的扩散，将严重影响计算结果的精度。实际计算时必须避免这类情况的发生。主元素消去法可以弥补这一缺陷，这里不再详细说明。

2.平方根法的源代码程序较长，空间复杂度高；但是需要的乘除运算次数为n3/6数量级，比LU分解节省近一半的工作量。而且，平方根法所求得的中间量lik是完全可以控制的，故舍入误差的增长也是可控的，因而计算过程是稳定的。

3.追赶法顺序求y的消元过程称为“追”的过程，逆序求x的回代过程称为“赶”的过程，追赶法因此而得名。可以看到，虽然追赶法的原理与高斯消去法相同，但由于三对角方程组的特殊性，计算时将系数中的大量零元素撇开，从而大大节省了计算量与存储量。因此，在实际编制程序时，只需要3个一维数组分别存储三对角方程组系数矩阵A的元素{ai}、{bi}、{ci}即可。

4.稀疏线性方程组解法的特点。

稀疏线性方程组指的是对应矩阵的大多数分量都是零元素, 非零元素只占其中少部分。 我们知道线性方程组可写成 [A] {x} = {b} 的形式. 在很多数值方法 (比如有限元法) 得到的矩阵 A 往往是一个稀疏矩阵, 如果对其中每一个元素都要储存和计算的话, 将会浪费很大的内存和CPU计算时间, 合理利用线性方程组的稀疏性可以节省存储空间和计算时间。

5.模块化结构程序设计方法的优点。

每个模块都可以分配给不同的程序员完成，从而缩短开发周期；各个模块高聚合、模块之间低耦合，只要模块之间确定了参数传递的接口，不管哪个模块内部的改动，均不会影响其它模块， 从而使软件产品的生产更加灵活； 系统细化到模块，条理清楚，系统更加容易理解和实现；容易维护、系统可靠。特点：

·各模块相对、功能单一、结构清晰、接口简单； 

·避免程序开发的重复劳动； 

·易于维护和功能扩充； 

·程序设计的复杂性得到了有效控制等。

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

一．问题提出

对于下面的方程组，试分别选用Jacobi迭代法，Gauss-Seidol迭代法和SOR方法计算

其解。

          10x1  -  x2   -2x3 = 7.2

           -x1  +10x2   -2x3  = 8.3

           -x1  -  x2   +5x3 = 4.2

2．要求

  1.体会迭代法求解线性方程组，并能与消去法做以比较；

  2.分别对不同精度要求，如ε=10-3,10-4,10-5由迭代次数体会该迭代法的收敛快慢；

  3.对方程组2,3使用SOR方法时，选取松弛因子w=0.8,0.9,1,1.1,1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者；

  4.给出各种算法的设计程序和计算结果。

三．实验目的和意义

1.通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较；

2.运用所学的迭代算法，解决各类线性方程组，编出算法程序；

3.体会上机计算时，终止步骤||x(k+1)-x(k)||∞ ＜ ε或k＞(给予的迭代次数),对迭代法敛散性的意义；

4.体会初始解x(0)，松弛因子的选取，对计算结果的影响。

四．计算公式

  1.Jacobi迭代法(aii≠0；i=1,2,…,n；j=1,2,…,n；k=0,1,2,…)

  2.Gauss-Seidol迭代法(aii≠0；i=1,2,…,n；j=1,2,…,n；k=0,1,2,…)

  3.SOR方法(aii≠0；i=1,2,…,n；j=1,2,…,n；k=0,1,2,…)

五．结构程序设计

1.Jacobi迭代法

//Jacobi迭代法解线性方程组

//编制人：尹一君 20134651

#include"stdio.h"

#include"math.h"

void main()

{

    FILE *f;

    double a[5][5],b[5],z[5],x[5]={0,0,0,0,0};

    double s,y,eps=1e-10;

    int i,j,k=0;                                 //k用于统计迭代次数

    if((f=fopen("Array6.txt

        return;

for(i=1;i<=3;i++)

        for(j=1;j<=3;j++)

            fscanf(f,"%lf",&a[i][j]);

for(i=1;i<=3;i++)

        fscanf(f,"%lf",&b[i]);

    while(1)                                

    {

        for(i=1;i<=3;i++)

        {

            s=0;

            for(j=1;j<=3;j++)

                if(j!=i)

                    s=s+a[i][j]*x[j];

            y=x[i];   

            z[i]=(b[i]-s)/a[i][i];   

if(fabs(y-z[i])<=eps) goto loop;

        }

        for(i=1;i<=3;i++)

            x[i]=z[i];

        k++;

    }

loop:

for(i=1;i<=3;i++)

        printf("x[%d]=%f\\n",i,x[i]);

        printf("\\nJacobi迭代法----迭代次数为：%d次\\n\\n",k);

}

2.Gauss-Seidol迭代法

//Gauss-Seidol迭代法解线性方程组

//编制人：尹一君 20134651

#include"stdio.h"

#include"math.h"

void main()

{

    FILE *f;

    double a[5][5],b[5],x[5]={0,0,0,0,0};

    double s,y,eps=1e-10;

    int i,j,k=0;                                 //k用于统计迭代次数

    if((f=fopen("Array6.txt

        return;

for(i=1;i<=3;i++)

        for(j=1;j<=3;j++)

            fscanf(f,"%lf",&a[i][j]);

for(i=1;i<=3;i++)

        fscanf(f,"%lf",&b[i]);

    while(1)                                

    {

        for(i=1;i<=3;i++)

        {

            s=0;

            for(j=1;j<=3;j++)

                if(j!=i)   

                    s=s+a[i][j]*x[j];

            y=x[i];   

            x[i]=(b[i]-s)/a[i][i];   

            if(fabs(y-x[i])<=eps)  goto loop;

        }

        k++;

    }

loop:

for(i=1;i<=3;i++)

        printf("x[%d]=%f\\n",i,x[i]);

        printf("\\nGauss-Seidol迭代法----迭代次数为：%d次\\n\\n",k);

}

3.SOR方法

//SOR方法解线性方程组

//编制人：尹一君 20134651

#define w 1.05                        //松弛因子，已通过多次试验得出其最佳值为1.05

#include"stdio.h"

#include"math.h"

void main()

{

    FILE *f;

    double a[5][5],b[5],x[5]={0,0,0,0,0};

    double s,y,eps=1e-10;

    int i,j,k=0;                        //k用于统计迭代次数

    if((f=fopen("Array6.txt

        return;

for(i=1;i<=3;i++)

        for(j=1;j<=3;j++)

            fscanf(f,"%lf",&a[i][j]);

for(i=1;i<=3;i++)

        fscanf(f,"%lf",&b[i]);

    while(1)                                

    {

        for(i=1;i<=3;i++)

        {

            s=0;

            for(j=1;j<=3;j++)

                if(j!=i)   

                    s=s+a[i][j]*x[j];

            y=x[i];   

            x[i]=(1-w)*x[i]+w*(b[i]-s)/a[i][i];   

            if(fabs(y-x[i])<=eps)  goto loop;

        }

        k++;

    }

loop:

for(i=1;i<=3;i++)

        printf("x[%d]=%f\\n",i,x[i]);

        printf("\\nSOR方法----迭代次数为：%d次\\n\\n",k);

}

六．结果展示

   问题1、问题2和问题3的程序运行结果分别如图4-1、4-2和4-3所示：

                                图4-1

                                图4-2

                                图4-3

七．实验讨论与总结

1.求解线性方程组的总结。

求解线性方程组的方法主要有两类，即直接解法和迭代解法。直接法比较适合于方程组的系数矩阵是满秩矩阵，且系数矩阵阶数不太高的情况。在许多实际问题中引出的线性方程组，往往是稀疏的，即系数矩阵中的大多数元素为零。对于这类方程组，如果它有不具有带状性，若采用直接法，消元过程会使大量零元素转化为非零元素，并破坏系数矩阵的形状，导致计算量的增加和存储空间的浪费。因此，对于此类稀疏问题，往往采用另一种方法：迭代法。

2.迭代法与消去法的区别。

迭代法与消去法的区别在于它不是通过有限次的算术运算求得方程组的精确解，而是通过迭代逐步逼近方程组的精确解。也就是从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法。

3.松弛因子

    在使用SOR方法求解方程组时，松弛因子w的选取对算法的迭代次数产生了很大的影响。为了保证迭代过程收敛，常要求0＜w＜2。当0＜w＜1时，称为低松弛法；当1＜w＜2时，称为超松弛法。一般统称超松弛法。迭代次数与w值的关系，类似一个开口向上的抛物线，一定存在一个w值，使迭代次数达到最少，而不论w在此基础上增大或是减小，都会使迭代次数同程度的增加。

4.终止步骤对迭代敛散性的影响。

终止步骤||x(k+1)-x(k)||∞ ＜ ε，对迭代法敛散性的影响在于迭代次数和迭代精度。若选取的误差值达到足够小，则必然伴随着迭代次数的增加，但另一方面，迭代精度也会达到更高的要求。若选取迭代次数k作为终止条件也是一样，指定的迭代次数越大，越可能达到所需精度要求。因此，两种方法可描述为：控精度，变次数；控次数，变精度。

实验五 函数插值方法

1．问题提出

2．要求

  1.利用Lagrange插值公式编写出插值多项式程序；

  2.给出插值多项式或分段三次插值多项式的表达式；

  3.根据节点选取原则，对问题2用三点插值或二点插值，其结果如何；

  4.对此插值问题用Newton差值多项式其结果如何。

3．实验目的和意义 

  1.学会常用的插值方法，求函数的近似表达式，以解决其它实际问题；

  2.明确插值多项式和分段插值多项式各自的优缺点；

  3.熟悉插值方法的程序编制；

  4.如果绘出插值函数的曲线，观察其光滑性。 

4．计算公式

  在第二目要求中Lagrange公式已给出，此处只给出分段线性插值的计算公式：

  设f(x)在n+1个节点a=x05．结构程序设计

  对于问题（1），首先采用Lagrange插值公式为算法原理，并直接套用公式编制程序：

#include"stdio.h"

#include"math.h"

double x[6]={0.4,0.55,0.65,0.80,0.95,1.05};

double y[6]={0.41075,0.57815,0.69675,0.90,1.00,1.25382};

double l0(double X)

{

    return ((X-x[1])*(X-x[2])*(X-x[3])*(X-x[4])*(X-x[5]))/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4])*(x[0]-x[5]));

}

double l1(double X)

{

    return ((X-x[0])*(X-x[2])*(X-x[3])*(X-x[4])*(X-x[5]))/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4])*(x[1]-x[5]));

}

double l2(double X)

{

    return ((X-x[0])*(X-x[1])*(X-x[3])*(X-x[4])*(X-x[5]))/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4])*(x[2]-x[5]));

}

double l3(double X)

{

    return ((X-x[0])*(X-x[1])*(X-x[2])*(X-x[4])*(X-x[5]))/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4])*(x[3]-x[5]));

}

double l4(double X)

{

    return ((X-x[0])*(X-x[1])*(X-x[2])*(X-x[3])*(X-x[5]))/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3])*(x[4]-x[5]));

}

double l5(double X)

{

    return ((X-x[0])*(X-x[1])*(X-x[2])*(X-x[3])*(X-x[4]))/((x[5]-x[0])*(x[5]-x[1])*(x[5]-x[2])*(x[5]-x[3])*(x[5]-x[4]));

}

double L5(double X)

{

    return y[0]*l0(X)+y[1]*l1(X)+y[2]*l2(X)+y[3]*l3(X)+y[4]*l4(X)+y[5]*l5(X);

}

void main()

{

    printf("f(0.596)=%lf\\n",L5(0.596));

    printf("f(0.99)=%lf\\n",L5(0.99));

}

  由此看出，此程序算法虽然直白，但未免过于繁琐复杂，因此考虑到将所有插值基函数写到一个for语句里面。改进算法在问题（2）中展现：

#include"stdio.h"

#include"math.h"

double x[7]={1,2,3,4,5,6,7};

double y[7]={0.368,0.135,0.050,0.018,0.007,0.002,0.001};

double L6(double X)

{

    int i,j;

    double l6=0,p,q;

for(i=0;i<=6;i++)

    {

        p=q=1;

        for(j=0;j<=6;j++)

            if(j!=i)

            {

                p=p*(X-x[j]);

                q=q*(x[i]-x[j]);

            }

        l6+=y[i]*p/q;

    }

    return l6;

}

void main()

{

    printf("f(1.8)=%lf\\n",L6(1.8));

    printf("f(6.15)=%lf\\n",L6(6.15));

}

6．结果展示

  问题（1）、（2）的程序运行结果如图6-1和图6-2所示：

                图6-1                                图6-2

7．实验讨论与总结

1.用三点插值或两点插值的讨论。

三点插值或两点插值，无非是选取三个节点或两个节点进行插值多项式构造的过程，那么得到的插值多项式为二次多项式（或一次多项式），利用这样的插值多项式计算函数值，其误差必然会更大。

2.对此插值问题用牛顿插值法其结果如何？

牛顿插值多项式与拉格朗日相比，不仅克服了“增加一个节点时整个计算工作重新开始

”的缺点，而且每增加一个节点时只需增加相应的一项即可，节省了乘除法运算次数。然而，插值多项式具有唯一性，因此不论用牛顿插值法还是拉格朗日插值法，得到的插值多项式是完全一样的，因此计算结果也是相同的。

3.分段线性插值法的讨论

由于高次插值多项式可能会产生龙格现象，而且从舍入误差看，高次插值误差的传播也较为严重。为克服在区间上进行高次插值所造成的龙格现象，采用分段插值的方法，将插值区间分成若干个小的区间，在每个小区间进行线性插值，然后相互连接，用连接相邻节点的折线逼近被插函数，这种把插值区间分段的方法就是分段线性插值法。

实验六 函数逼近与曲线拟合

1．问题提出

2．要求

  1.用最小二乘法进行曲线拟合；

  2.近似解析表达式为φ(t)=a1t+a2t2+a3t3；

  3.打印出拟合函数φ(t)，并打印出φ(tj)与y(tj)的误差，j=1,2,…,12；

  4.另外选取一个近似表达式，尝试拟合效果的比较；

  5.*绘制出曲线拟合图。

3．实验目的和意义

  1.掌握曲线拟合的最小二乘法；

  2.最小二乘法亦可用于解超定线代数方程组；

  3.探索拟合函数的选择与拟合精度间的关系。

四．计算公式

5．结构程序设计

#include "stdio.h"

#include "math.h"

#define num 10

float neiji(float b[num],float c[num]) /*内积函数*/

{  

    int p;

    float nj=0;

for (p=1;p       nj+=c[p]*b[p];

    return nj;

}

float s[num],x[num],y[num],fai[num][num],afa[num];

float beida[num],a[num],xfai[num],yd[num],max,pcpfh;

void main()

{ 

    int i,j,k,n,index,flag;

    char conti;

    conti=' ';

    printf("请输入已知点的个数n=");

    scanf("%d",&n);

    printf("请输入x 和y\\n:");

for(i=1;i<=n;i++)

{

     printf("x[%d]=",i);

         scanf("%f",&x[i]);

         printf("y[%d]=",i);

         scanf("%f",&y[i]);

    }

    while(conti==' ')

{ 

    printf("请输入拟和次数=");

        scanf("%d",&index);

        pcpfh=0;

        afa[1]=0;

        a[0]=0;

for(i=1;i<=n;i++)

        {

            afa[1]+=x[i];

            a[0]+=y[i];

            fai[0][i]=1;

        }

        afa[1]=afa[1]/n;

        a[0]=a[0]/n;

for (i=1;i<=n;i++)

        {

            fai[1][i]=x[i]-afa[1];

        }

        a[1]=neiji(fai[1],y)/neiji(fai[1],fai[1]);

for(k=1;k        {

for(i=1;i<=n;i++)

            xfai[i]=x[i]*fai[k][i];

            afa[k+1]=neiji(fai[k],xfai)/neiji(fai[k],fai[k]);

            beida[k]=neiji(fai[k],fai[k])/neiji(fai[k-1],fai[k-1]);

for(j=1;j<=n;j++)

                fai[k+1][j]=(x[j]-afa[k+1])*fai[k][j]-beida[k]*fai[k-1][j];

            a[k+1]=neiji(fai[k+1],y)/neiji(fai[k+1],fai[k+1]);

        }

        printf("%d 次拟和结果为\\n",index);

for(i=0;i<=index;i++)

        printf("a[%d]=%f\\n",i,a[i]);

for(i=1;i<=index;i++)

        printf("afa[%d]=%f\\n",i,afa[i]);

for(i=1;i            printf("beida[%d]=%f\\n",i,beida[i]);

for(i=1;i<=n;i++)

        { 

for(k=0;k<=index;k++)

            s[i]+=a[k]*fai[k][i];

            yd[i]=float(fabs(y[i]-s[i]));

            pcpfh+=yd[i]*yd[i];

            s[i]=0;

        }

        max=0;

for(i=1;i<=n;i++)

if(yd[i]>max)

        {

            max=yd[i];

            flag=i;

        }

        printf("当x=%f 时,偏差最大=%f,偏差平方和为%f\\n",x[flag],max,pcpfh);

        printf("继续拟和请按space,按其他键退出");

        conti=getchar();

        conti=getchar();

    }

}

6．结果展示

  该程序的运行结果如图6-1所示：

                                图6-1

7．实验讨论与总结

1.曲线拟合的最小二乘法是处理实验数据的常用方法。多项式拟合是线性最小二乘法拟合问题的一种特殊情况，其特点是拟合多项式形式简单，但当n较大时，法方程组往往是病态的。用离散正交多项式进行曲线拟合，不用解线性方程组，只需按递推式进行计算，避免了法方程组病态所造成的麻烦，并且当逼近次数增加一次时，只要在原基础上增加一项，计算程序十分简单。

2.关于非线性最小二乘曲线拟合问题，一般求解比较困难，但对一些特殊情形，可以转换为线性最小二乘拟合问题。在实际计算时，要选择合理的拟合多项式的次数有时十分困难的。一般可对数据作分析，例如在方格纸上作草图，从草图中观察应作几次多项式精度较好，以选择最佳的拟合多项式的次数。

3.近似表达式的拟合效果。

由于本实验所选取的近似解析表达式为φ(t)=a1t+a2t2+a3t3，这样通过执行得出的拟合曲线为y=-0.0052t2+0.2634t+0.0178，拟合效果可如下图6-2表示： 

                             图6-2

    由此得出拟合效果在后期出现病态。因此根据所求，另外选取一个二次的近似表达式，即φ(t)=a0+a1t+a2t2，这样通过执行得出的拟合曲线为y=-0.0024t2+0.2037t+0.2305，拟合效果可如下图6-3表示：

                             图6-3

从图中看出，该表达式的拟合效果非常好。由此可见，在利用数据的最小二乘法求拟合曲线时，选取合适的近似表达式很重要，应通过不断的试验找出较为合适的近似表达式，这样才能尽可能的提高拟合精度。

实验七 数值积分与数值微分

1．问题提出

  选用复合梯形公式，复合Simpson公式，Romberg算法，计算

 （1）

 （2）

 （3）

 （4）

2．要求

  1.编制数值积分算法的程序；

  2.分别用两种算法计算同一个积分，并比较其结果；

  3.分别取不同步长h=（b-a）/n，试比较计算结果（如n=10,20等）；

  4.给定精度要求，试用变步长算法，确定最佳步长。

3．实验目的和意义

  1.深刻认识数值积分法的意义；

  2.明确数值积分精度与步长的关系；

  3.根据定积分的计算方法，可以考虑二重积分的计算问题。

4．计算公式

  1.复化梯形公式

  2.复化辛普森公式

  3.龙贝格算法（参照教材P139）

5．结构程序设计

  首先对于积分（1）采用复合梯形公式，算法如下：

#define n 10

#include"stdio.h"

#include"math.h"

double f1(double x)

{

    return pow(4-pow(sin(x),2),1.0/2);

}

double T1()

{

    int k;

    double T=0;

    double a=0,b=1.0/4,h;

    h=(b-a)/n;

for(k=0;k<=n;k++)

    {

        if(k==0)

            T+=f1(a);

        else if(k==n)

            T+=f1(b);

        else

            T+=2*f1(a+k*h);

    }

    T=T*(h/2);

    return T;

}

void main()

{

    printf("该定积分的值为%lf（复化梯形公式）\\n",T1());

}

  然后，使用复化辛普森公式计算同一个定积分。只需在原来的基础上增加一个子函数S1，该函数建立如下：

double S1()

{

    int i,k;

    double S=0;

    double a=0,b=1.0/4,h;

    h=(b-a)/n;

for(i=1,k=0;k<=2*n;k++)

    {

        if(k==0)

            S+=f1(a);

        else if(k==2*n)

            S+=f1(b);

        else if(k==2*i-1)

        {S+=4*f1(a+k*(h/2));i++;}

        else

            S+=2*f1(a+k*(h/2));

    }

    S=S*(h/6);

    return S;

}

  接下来，我们再对积分（2）分别采用这两种方法编制程序，需要注意积分（2）对应的函数在x=0时函数值为1,。程序如下：

#define n 10

#include"stdio.h"

#include"math.h"

double f2(double x)

{

    return sin(x)/x;

}

double T2()

{

    int k;

    double T=0;

    double a=0,b=1,h;

    h=(b-a)/n;

for(k=0;k<=n;k++)

    {

        if(k==0)

            T+=1;

        else if(k==n)

            T+=f2(b);

        else

            T+=2*f2(a+k*h);

    }

    T=T*(h/2);

    return T;

}

double S2()

{

    int i,k;

    double S=0;

    double a=0,b=1,h;

    h=(b-a)/n;

for(i=1,k=0;k<=2*n;k++)

    {

        if(k==0)

            S+=1;

        else if(k==2*n)

            S+=f2(b);

        else if(k==2*i-1)

        {S+=4*f2(a+k*(h/2));i++;}

        else

            S+=2*f2(a+k*(h/2));

    }

    S=S*(h/6);

    return S;

}

void main()

{

    printf("该定积分的值为%lf（复化梯形公式）\\n",T2());

    printf("该定积分的值为%lf（复化辛普森公式）\\n",S2());

}

6．结果展示

    积分（3）和（4）的算法设计与上述算法差别仅在函数的建立上，此处算法从略，但仍给出执行结果。积分（1）、（2）、（3）和（4）使用复化梯形公式、复化辛普森公式得出的计算结果如图7-1、图7-2、图7-3和图7-4所示：

    图7-1                                    图7-2

     图7-3                                    图7-4

7．实验讨论与总结

1.数值积分精度与步长的关系。

通过简单的分析，我们认为通过对h值的改变，只要h值越小，即等分的区间越小，结果应该更加精确，精确度越高。经过实验的验证，也表明我们的推理正确，无论是复合梯形公式还是复合辛普森公式它们最终结果都会随着h值的减小而更加精确。复合梯形公式和复合辛普森公式计算出的结果进行比较，发现复合辛普森公式计算出的结果更加的精确。

2.变步长梯形公式

复化求积方法对于提高计算精度是行之有效的方法，但复化公式的一个主要缺点在于要先估计出步长。若步长太大，则难以保持计算精度，若步长太小，则计算量太大，并且积累误差也会增大。在实际计算中通常采用变步长的方法，即把步长逐次分半，直至达到某种精度为止。变步长梯形公式为：

    3.龙贝格算法

然而，变步长梯形求积法虽然算法简单，但收敛速度较慢。但可以利用梯形法算法简单的优点，形成一个新算法，这就是龙贝格算法（详见教材P139），龙贝格算法在区间逐次分半过程中，对用梯形法所获得的近似值进行多级“加工”，从而获得高精度的积分近似值的一种方法。它具有自动选取步长且精度高，计算量小的特点，便于在计算机上使用，是数值积分中较好的方法。下面以积分（2）为例，给出龙贝格算法的程序实现：

#include"iostream.h" 

#include"math.h"  

#define e  0.00000000000001  

double  f(double x)  

{   

    double y;   

    if (x==0)   

        return y=1.0;     

    else        

        y=sin(x)/x;   

    return y; 

}  

void romberg(double a,double b)  

{     

    int n=1,k=0;     

    double h,T2,S2=0,C2=0,R2=0,T1,C1,S1,R1;     

    h=(b-a)/2;     T2=h*(f(a)+f(b));   

while (fabs((R2-R1))>e)

    {     

        R1=R2; T1=T2; S1=S2; C1=C2;     

        double sum=0;     

        int i;     

for(i=1;i<=n;i++)

            sum=sum+f(a+(2*i-1)*h);          

        T2=T1/2+sum*h;    

        S2=(4*T2-T1)/3;    

        C2=(16*S2-S1)/15;     

        R2=(*C2-C1)/63;    

        n=n*2; k++; h=h/2;    

    }   

cout<<"I="<  

void main()   

{    

    double a,b;    

    a=0;b=1;  

    cout<<"龙贝格算法结果如下："<    romberg(a,b);      

}  

程序运行结果如图7-5所示：

              图7-5

参考文献

◆ 《数值计算方法》 合肥工业大学数学与信息科学系 编 合肥工业大学出版社 

◆ 《计算方法》 邓建中等编，西安交大出版社，1985。 

◆ 《数值计算和C程序集》蒋长锦编著，中国科学技术大学出版社，1998。 

◆ 《计算方法引论》徐萃薇编，高等教育出版社，1999。  

◆ 黄友谦，程诗杰，陈浙鹏，《数值试验》，北京：高等教育出版社，19 

◆ 蔡大用，《数值分析与实验学习指导》，北京：清华大学出版社与施普林格出版社，2001 

◆ 肖筱南，《值计算方法与上机实习指导》，北京：北京大学出版社，2004 

◆  A. Quarteroni, R. Sacco, F. Saleri,《Numerical Mathematics》, New York：Springer-Verlag, 2000