现将小生研一数值分析作业摘录如下，上机的代老师要求我们分为四个部分：解题理论依据、计算程序、计算结果截图和问题讨论，因公式输入不便所以第一部分的公式略去，需要的可去第一篇博客“研究生数值分析课后题（上机编程）-1 ”下载附件。因为本人编程水平有限，故使用的是最基本的函数，但保证原创和结果正确，所有程序在VC++6.0上验证通过，还望各位兄台不吝赐教。

 

            第四题：牛顿法求方程近似根

1.解题理论依据或方法应用条件:

依据牛顿迭代法进行计算

2.计算程序（使用软件：VC）：

#include<math.h>

void main()

{double x0,x1=1.9,f0,f1;

for(x0=1.0;fabs(x0-x1)>=0.00001;)

{x0=x1;

 f0=pow(x0,7)-28*pow(x0,4)+14;

 f1=7*pow(x0,6)-112*pow(x0,3);

 if(f1==0)break;

 else x1=x0-f0/f1;}

if(f1==0)printf("出错\n");

else printf("方程的解为：%.16lf\n",x1);

}

3.计算结果

 4.问题讨论（误差分析、上机出现情况等）

（1）在做除法运算的时候，应对除数是否等于零进行判断。        

（2）编程的过程出现问题如下：①在编写程序过程中，不小心将输入换成了全角输入造成错误②忘记头文件中增加#include<math.h>，导致pow函数无法使用。

 

     第五题：用Romberg算法求积分近似值

1.解题理论依据或方法应用条件:

依据Romberg法列出相应的方程，先求出第一列的相应元素，然后利用Romberg中三角法求出后面列的元素，将得到的第一行元素取差的绝对值，当小于误差要求的时候可以停止运算

2.计算程序（使用软件：VC）：

方法一：

#include<stdio.h>

#include<math.h>

#define f(x) pow(3,(x))*pow((x),1.4)*(5*(x)+7)*sin((x)*(x))

void main()

{double t[9][9];            

 double sum=0;

 int j=1,i=1,m=0,k=1;

 double a=1,b=3;

t[0][1]=(b-a)/2*(f(a)+f(b));

sum=f(a+(2*1-1)*(b-a)/pow(2,1));

 t[1][1]=0.5*(t[0][1]+sum*(b-a)/pow(2,0));

 t[0][2]=(pow(4,1)*t[1][1]-t[0][1])/(pow(4,1)-1);       /*先求出t[0][1]，t[0][2]的值*/

 for(j=1;fabs(t[0][j+1]-t[0][j])>0.00001;)              /*判断误差是否符合条件*/

 {j++;

  sum=0;                                      /*累加和清零*/

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

  sum+=f(a+(2*i-1)*(b-a)/pow(2,j));                 /*求出式子中的累加清零*/

  t[j][1]=0.5*(t[j-1][1]+sum*(b-a)/pow(2,j-1));         /*求出第一列元素的值*/

  for(k=j;k>=1;k--)

  {m++;t[k-1][m+1]=(pow(4,m)*t[k][m]-t[k-1][m])/(pow(4,m)-1);} 

/*依次用三角形法求出其他元素的值*/

   m=0;

 }

   printf("积分的值为：%.16lf\n",t[0][j+1]);

}

方法二：（逐列求值，较浪费电脑资源、时间）

#include<stdio.h>

#include<math.h>

#define f(x) pow(3,(x))*pow((x),1.4)*(5*(x)+7)*sin((x)*(x))

void main()

{double t[9][9];            

 double sum;

 int j=1,i=1,m=1,k=1,l=1;

 double a=1,b=3;

t[0][1]=(b-a)/2*(f(a)+f(b));

 for(l=1;l<=7;l++)

 {sum=0;

 for(i=1;i<=pow(2,l-1);i++)

 sum+=f(a+(2*i-1)*(b-a)/pow(2,l));

 t[l][1]=0.5*(t[l-1][1]+sum*(b-a)/pow(2,l-1));

 }

for(m=1;m<=7;m++)

 for(k=1;k<=(7-m+1);k++)

 t[k-1][m+1]=(pow(4,m)*t[k][m]-t[k-1][m])/(pow(4,m)-1);

for(j=1;fabs(t[0][j+1]-t[0][j])>0.00001;j++);

 printf("积分的值为：%.16lf\n",t[0][j+1]);

}

3.计算结果

4.问题讨论（误差分析、上机出现情况等）

（1）这道题目是前三道题目中耗时最长、但也是学到最多东西的题目，一开始使用的是第二个方法进行编程，VC报错，一直未找到哪里出错，最后发现原来是将

int j=1,i=1,m=1,k=1,l=1; 放在了 t[0][1]=(b-a)/2*(f(a)+f(b)); 之后，VC规定必须将所有的定义放在程序开头，这是自己以前所不知道的。

（2）改掉上面的错误后，系统不再报错，但是数值不对，逐行看程序发现原来是sum未清零，导致后面的数据出错。

（3）改掉上面的错误后，得到的数据依然不对，纠缠很久以及翻阅C语言书籍发现原来是#define f(x) pow(3,(x))*pow((x),1.4)*(5*(x)+7)*sin((x)*(x)) 语句中，X必须加上括号，这样才能保证运算顺序正确。

（4）这样程序的结果正确了，但是第二种方法较浪费电脑资源、时间，故对原程序更改有了第一种利用循环、用多少数据求多少数据的方法。

 

    

        第六题：定步长四阶Runge-Kutta法求方程组

1.解题理论依据或方法应用条件:

依据定步长四阶Runge-Kutta法进行计算

2.计算程序（使用软件：VC）：

#include<stdio.h>

main()

{float y[4][205]={{0},{0},{0},{0}};                  /*定义第一列值全为零*/

 int i,j;

 float k1,k2,k3,k4;

 float h=0.0005;                                                     

for(j=1;(j*0.0005)<=0.1;j++)

 for(i=1;i<=3;i++)                              /*本解法未使用0行，故从1开始取值*/

{if(i=1)

 {k1=1;k2=1;k3=1;k4=1;

  y[1][j]=y[1][j-1]+h/6*(k1+2*k2+2*k3+k4);}        /*计算第一行数值*/

 if(i=2)

 {k1=y[3][j-1];

  k2=y[3][j-1]+0.5*h*(1000-1000*y[2][j-1]-100*y[3][j-1]);

k3=y[3][j-1]+0.5*h*(1000-1000*(y[2][j-1]+0.5*h*y[3][j-1])-100*(y[3][j-1]+0.5*h*(1000-1000*y[2][j-1]-100*y[3][j-1])));

k4=y[3][j-1]+h*(1000-1000*(y[2][j-1]+0.5*h*(y[3][j-1]+0.5*h*(1000-1000*y[2][j-1]-100*y[3][j-1])))-100*(y[3][j-1]+0.5*h*(1000-1000*(y[2][j-1]+0.5*h*y[3][j-1])-100*(y[3][j-1]+0.5*h*(1000-1000*y[2][j-1]-100*y[3][j-1])))));

  y[2][j]=y[2][j-1]+h/6*(k1+2*k2+2*k3+k4);}         /*计算第二行数值*/

 if(i=3)

 {k1=1000-1000*y[2][j-1]-100*y[3][j-1];

  k2=1000-1000*(y[2][j-1]+0.5*h*y[3][j-1])-100*(y[3][j-1]+0.5*h*(1000-1000*y[2][j-1]-100*y[3][j-1]));

k3=1000-1000*(y[2][j-1]+0.5*h*(y[3][j-1]+0.5*h*(1000-1000*y[2][j-1]-100*y[3][j-1])))-100*(y[3][j-1]+0.5*h*(1000-1000*(y[2][j-1]+0.5*h*y[3][j-1])-100*(y[3][j-1]+0.5*h*(1000-1000*y[2][j-1]-100*y[3][j-1]))));

k4=1000-1000*(y[2][j-1]+h*(y[3][j-1]+0.5*h*(1000-1000*(y[2][j-1]+0.5*h*y[3][j-1])-100*(y[3][j-1]+0.5*h*(1000-1000*y[2][j-1]-100*y[3][j-1])))))-100*(y[3][j-1]+h*k3);

y[3][j]=y[3][j-1]+h/6*(k1+2*k2+2*k3+k4);}                        /*计算第三行数值*/

}

 printf("y1(0.025),y2(0.025),y3(0.025)值分别为：%f,%f,%f\n",y[1][50],y[2][50],y[3][50]);

 printf("y1(0.045),y2(0.045),y3(0.045)值分别为：%f,%f,%f\n",y[1][90],y[2][90],y[3][90]);

 printf("y1(0.085),y2(0.085),y3(0.085)值分别为：%f,%f,%f\n",y[1][170],y[2][170],y[3][170]);

 printf("y1(0.100),y2(0.100),y3(0.100)值分别为：%f,%f,%f\n",y[1][200],y[2][200],y[3][200]);                

}

3.计算结果

 4.问题讨论（误差分析、上机出现情况等）

（1）本题的编写过程相对前一题较为顺利，过程中只出现一个较大问题，原程序中我在for循环体内使用的是if(){}  else if(){}  else (){} 的结构，无错，但运行不出结果，改用if (){}  if (){}  if (){} 后出现结果，至今仍未找清楚原因。

（2）另一个小问题是，在计算第三行数据时，由于要使用前一列二行的数据，刚开始直接将复制进来，稍后发现这时的K值是二行中的数据，放在第三行中将出错，由于发现及时，未对后面程序编写造成影响。（故后续将K值都由具体的算式替代，避免出错）