研究生数值分析课后题（上机编程）-2

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

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

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

依据牛顿迭代法进行计算

2.[/b]计算程序（使用软件：VC[/b]）：[/b]

#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.[/b]计算结果[/b]


[/b]
4.问题讨论（误差分析、上机出现情况等）
（1）在做除法运算的时候，应对除数是否等于零进行判断。

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

第五题：用[/b]Romberg[/b]算法求积分近似值[/b][/b]

1.[/b]解题理论依据或方法应用条件:[/b]
依据Romberg法列出相应的方程，先求出第一列的相应元素，然后利用Romberg中三角法求出后面列的元素，将得到的第一行元素取差的绝对值，当小于误差要求的时候可以停止运算
2.[/b]计算程序（使用软件：VC[/b]）：[/b]
方法一：[/b]
#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]);
}
方法二：（逐列求值，较浪费电脑资源、时间）[/b]
#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.[/b]计算结果[/b]


[/b]
4.[/b]问题讨论（误差分析、上机出现情况等）[/b]
（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）这样程序的结果正确了，但是第二种方法较浪费电脑资源、时间，故对原程序更改有了第一种利用循环、用多少数据求多少数据的方法。

第六题：定步长四阶[/b]Runge-Kutta[/b]法求方程组[/b][/b]
1.[/b]解题理论依据或方法应用条件:[/b]
依据定步长四阶Runge-Kutta法进行计算

2.[/b]计算程序（使用软件：VC[/b]）：[/b]
#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.[/b]计算结果[/b]





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

（1）本题的编写过程相对前一题较为顺利，过程中只出现一个较大问题，原程序中我在for循环体内使用的是if(){} else if(){} else (){} 的结构，无错，但运行不出结果，改用if (){} if (){} if (){} 后出现结果，至今仍未找清楚原因。
（2）另一个小问题是，在计算第三行数据时，由于要使用前一列二行的数据，刚开始直接将复制进来，稍后发现这时的K值是二行中的数据，放在第三行中将出错，由于发现及时，未对后面程序编写造成影响。（故后续将K值都由具体的算式替代，避免出错）

本文出自 “马凯的技术博客” 博客，请务必保留此出处http://makai.blog.51cto.com/2205493/496814