数值计算方法实验报告

用两种方法求解以下积分： 

1.  ,  初值 

    2.  ,初值=1-.

1．1利用C语言编程求解，程序如下：

#include

void main()

{

    double a=100,b,t,I;

    int n;

for (n=50;n>=0;n--)

    { b=(1-a)/n;

      t=b;b=a;a=t;

    printf("n=%d,I=%lf     ",n,b);

    }

}

1.2.命令窗口中程序结果如下：

2.1利用C语言编程求解，程序如下：

#include

void  main()

{

    double a=0.63212056,b,t,I;

    int n;

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

{b=1-a*(n+1);

t=b;

b=a;

a=t;

printf("n=%d,I=%lf\\n",n,b);}

}

1.2.命令窗口中程序结果如下：

实验结果分析：

由运行就更可知，的结果不可靠。因为在计算时有舍入误差等，由于式子是要联想乘以n会将舍入误差放大。最后的结果由于误差放大了，所以结果不可靠。

   由运行可知，的结果可靠且接近直值，由于在式子中是乘以1/n，在连续的乘以1/n会将舍入误差放小，因而所得值可靠。

用牛顿迭代法和牛顿下山法求解方程f（x）=，取近似值。

1．牛顿迭代法：    ，直到。

2．牛顿下山法：   ，（λ=1至合适），且直到|

实验过程：

1．用牛顿迭代法求解

1.1利用C语言编程求解，程序如下：

#include

#include

double f(double x)

{

    double f;

    f=pow(x,3)/3-x;

    return f;

}

double g(double x)

{

    double g;

    g=pow(x,2)-1;

    return g;

}

void main()

{

    double a=-0.99,b,t;

    int n;

for (n=0;fabs(b-a)>1e-5;n++)

    {

        b=a-f(a)/g(a);

        t=b;

        b=a;

        a=t;

        printf("n=%d,x=%lf\\n",n,b);

    }

}

1.2.命令窗口结果截屏如下：

2．用牛顿下山法求解

2.1利用C语言编程求解，程序如下：

#include

#include

double f(double x)

{

    double f;

    f=pow(x,3)/3-x;

    return f;

}

double g(double x)

{

    double g;

    g=pow(x,2)-1;

    return g;

}

void main()

{

    double a=-0.99,b=a-f(a)/g(a),t,k;

    int n,m;

for (n=0;fabs(b-a)>1e-5;n++)

{for(m=0;fabs(f(b))>fabs(f(a));m++)

    {

        b=a-pow(0.5,m)*f(a)/g(a);

        k=pow(0.5,m);

        printf("k=%lf,f(b)=%lf,x=%lf,n=%d\\n",k,f(b),b,n);

    }

    t=b;b=a;a=t;

    }

}

2.2.命令窗口结果截屏如下：

实验结果分析：

由两个结果对比可知，牛顿迭代法的迭代次数较多，但是简单可靠。可见牛顿下山法要比牛顿迭代法的收敛速度快的多。但是对函数有一定的要求。牛顿迭代要求函数能求导。

∙实验三   解线性代数方程组——列主元消元法

用列主元法解线性方程组

=

1. 利用C语言编程求解，程序如下：

#include "stdio.h"

#include "math.h"

main()

{

    int i,j,p,n=0;

    double b[5],temp,x[4],max,

        a[4][5]={{1.1348,3.8326,1.1651,3.4017,9.5342},

                {0.5301,1.7875,2.5330,1.5435,6.3941},

                {3.4129,4.9317,8.73,1.3142,18.4231},

                {1.2371,4.9998,10.6721,0.0147,16.9237}};

    printf("请确定你的矩阵\\n");

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

    {

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

        {

            printf("%lf    ",a[i][j]);

        }

        printf("\\n");

    }

  //找第一列的最大值//

    max=a[0][0];

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

    {

        if(max        {

            max=a[i][0];

        }

    }    

//找第一列中最大值索在的行//

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

    {

        if(max==a[i][0])

        {

            break;

        }

    }n=i-1;

  //将最大值一行换为第一行//

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

    {

        b[i]=a[n][i];

        a[n][i]=a[0][i];

        a[0][i]=b[i];

    }

    //第一列为0//

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

    {   

        temp=a[0][0]/a[i+1][0];

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

        {

            a[i+1][j]=temp*a[i+1][j]-a[0][j];

        }

    }

    printf("第一次列主消元后的矩阵\\n");

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

    {

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

        {

            printf("%-9.5lf    ",a[i][j]);

        }

        printf("\\n");

    }

   //找第二列的最大值//

    max=a[1][1];

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

    {

        if(max        {

            max=a[i][1];

        }

    }

    //找第二列中最大值索在的行//

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

    {

        if(max==a[i][1])

        {

            break;

        }

    }n=i;

   //将最大值行换为第二行//

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

    {

        b[j]=a[n][j];

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

        a[1][j]=b[j];

    }

 //第二列为0//

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

    {   

        temp=a[1][1]/a[i+1][1];

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

        {

            a[i+1][j]=temp*a[i+1][j]-a[1][j];

        }

    }

    printf("第二次列主消元后的矩阵\\n");

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

    {

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

        {

            printf("%-9.5lf    ",a[i][j]);

        }

        printf("\\n");

    }

  //找第 三列的最大值//

    max= fabs(a[2][2]);

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

    {

        if(max        {

            max=a[i][2];

        }

    }

    //找第三列中最大值索在的行//

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

    {

        if(max==a[i][2])

        {

            break;

        }

    }n=i-1;

   //将最大值行换为第三行//

for(j=2;j<5;j++)

    {

        b[j]=a[n][j];

        a[n][j]=a[2][j];

        a[2][j]=b[j];

    }

  //第三列为0的矩阵//

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

    {   

        temp=a[2][2]/a[i+1][2];

        for(j=2;j<5;j++)

        {

            a[i+1][j]=temp*a[i+1][j]-a[2][j];

        }

    }

    printf("第三次列主消元后的矩阵\\n");

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

    {

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

        {

            printf("%-9.5lf    ",a[i][j]);

        }

        printf("\\n");

    }

  //计算x的值//

for(p=3;p>=0;p--)

    {

        x[p]=a[p][4];

        for(i=3;i>p;i--)

        {

                x[p]=x[p]-a[p][i];    

        }

        x[p]=x[p]/a[p][p];

    }

    printf("\\n计算所得x的值为：\\n");

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

    {

        printf("x[%d]=%lf\\n",i+1,x[i]);

    }

}

2.命令窗口结果截屏如下：

实验结果分析：

这种方法是将矩阵化为三角方程在求解。所以约化主元不为0，虽然可以完成方程的求解，但是无法保证结果的可靠性。如果我们在选择整个列子中块中最大的元作为约化列主，就能解决上面可能发生的矛盾，我们就将每次次消元时的最大值作为消元列主。这样就解决了计算过程中的舍入误差。这样消去的到的三角阵就较为正确，所以列主消元得到的值更可靠。

 用追赶法解三角方程组Ax=d，， d=

实验过程：

1.利用C语言编程求解，程序如下:

#include

void main()

{ a[5][5]={{2,2,0,0,0},{-1,1,2,0,0},{0,-1,1,2,0},{0,0,-1,1,2},{0,0,0,-1,1}},

f[5]={6,7,9,11,1},

c[4],d[5],e[4],y[5],x[5];

 int i,n,k;

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

    {d[0]=a[0][0];e[i-1]=a[i][i-1];c[i]=a[i-1][i]/d[i];d[i]=a[i][i]-a[i][i-1]*c[i-1];}

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

    {    c[0]=a[0][1]/a[0][0];

        c[i]=a[i-1][i]/(e[i]*c[i-1]);}

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

     {y[0]=f[0]/a[0][0];

     y[n]=(f[n]/e[n]*y[n-1])/(a[n][n]-e[n]*c[n-1]);}

for(k=4;k>0;k--)

    {    x[5]=y[5];

    x[k]=y[k]-c[k]*x[k+1];

    printf("%lf",x[k]);}}

2.命令窗口结果截屏如下:

实验结果分析：

追赶法的主要问题是把矩阵分解为两个对角阵，用循环把每一个元素表示出来是重点，所以要很清晰的理解追赶法求解的步骤和思路。

已知直升机旋转机翼外形曲线轮廓线上的某些点及端点处的一阶导数值，是S0=1.88,S19=-0.046115

试计算该曲线的横坐标为

2,4,6,12,16,30,60,110,180,280,400,515

该点处的纵坐标（要求该曲线具有二阶光滑度）。

k