数值分析实验报告

实验课成绩 武汉理工大学

学 生 实 验 报 告 书

实验课程名称： 数 值 分 析 开 课 名 称：计算机科学与技术学院 指导老师姓名： 熊盛武

学 生 姓 名： 学生专业班级：

20011—— 2012学年 第 二 学期

第一次试验

（1）二分法计算流程图：

开始 求根区间[a,b]循环次数i,i=0 方程x*x*x-x-1=0 否 b-a>e 是 i++ x=(a+b)/2 否 a=x f(a)f(x)<0 是 b=x 是 b-a>e 否 t=(a+b)/2 输出近似根t，循环次数i 结束 简单迭代法算法流程图：

开始 定义函数s(t)=t+1 方程x*x*x-x-1=0 输入迭代初始值x0,控制精 度e.迭代次数i，i= 0

x1=pow(s(x0),1.0/3.0); 否 fabs(x1-x0)>e 是 i++; x1=pow(s(x0),1.0/3.0); x0=x1; x1=pow(s(x0),1.0/3是 fabs(x1-x0)>ee 否 x=(x1+x0)/2 输出近似根 x,迭代次数i 结束 （2）

（3）牛顿迭代法流程图：

开始 定义函数s(t)=t-(t*t*t-t-1)/(3*t*t-1) 迭代初始值x0,控制精度e，迭代次数i,i=0 x1=x0 否 fabs(s(x1)-x1)>e 是 i++; x2=s(x1); x1=x2; 是 fabs(s(x1)-x1)>e 否 x=(x1+x2)/2; 输出近似根x,迭代次数i 结束

（4）弦截法算法程序流程图：

开始 定义函数s(t)= t*t*t-1 输入迭代初始值x0 x1,控制精度e,循环次数i，i=0 否 fabs(x1-x0)>e 是 i++; x2=x1-s(x1)*(x1-x0)/(s(x1)-s(x0)); x0=x1; x1=x2; 是 fabs(x1-x0)>e 否 x=(x1+x2)/2; 输出近似根x迭代次数i 结束

二分法程序源代码：

#include #include

#define f(x) (x*x*x-x-1) using namespace std; int main() { int i; float x,t,a,b,e; cout<<\"请输入求根区间a,b\"<<\"控制变量e\"<>a>>b>>e; i=0; while ((b-a)>e) { i++; x=(a+b)/2; if (f(a)*f(x)<0) b=x; if (f(a)*f(x)>0) a=x; } t=(a+b)/2; cout<<\"在求根区间a,b间近似根t=\"<调试过程，实验结果及分析：

计算 x*x*x-x-1=0 在[1,2]内的近似根。 精度达到0.0001时，程

序

运

行

结

果

如

下

图

：

当精度达到0.00001时，程序运行如下图：

调试过程中如果把while ((b-a)>e)改为while ((b-a)由运行结果看出：

精度达到0.0001时，二分法次数为14。 精度达到0.00001时，二分法次数为17。 精度要求越高，二分法次数越多。 二分法缺点：逼近速度较慢。

二分法对有根区间逐步逼近得到更精确的近似根，二分法不采用等步长扫描，而是在有根区间[a,b]中，取中点，然后计算中点值，判断中点对应函数值和端点函数值乘积的正负来确定新的端点，当有根区间长度小于精确度时得到近似根。

如果没有区间长度和精度大小的判断，程序是一个死程序，得不到近似根。

简单迭代法程序源代码：

#include #include double s(double t) { return (t+1); }

using namespace std; int main() { int i; double x,x0,x1,e; cout<<\"请输入迭代初始值x0\"<<\和控制精度e\"<>x0>>e; i=0;

while (fabs(x1-x0)>e) { i++;

x1=pow(s(x0),1.0/3.0); x0=x1;

x1=pow(s(x0),1.0/3.0); }

x=(x1+x0)/2;

cout<<\"近似根x=\"<调试过程，实验结果及分析：

计算 x*x*x-x-1=0 在[1,2]内的近似根。 精度达到0.0001时，程

序

运

行

结

果

如

下

图

：

当精度达到0.00001时，程序运行如下图：

当循环体程序代码如下： x1=pow(s(x0),1/3); x0=x1;

x1=pow(s(x0),1/3);

时会出现程序运行结果的界面，但是输出的近似根是1，迭代次数是1，这是因为x1=pow(s(x0),1/3)与三次求根函数式的数据类型不符。导致 x1=pow(s(x0),1/3)不再是三次求根的函数式

由运行结果看出：

精度达到0.0001时，简单迭代法迭代次数为5. 精度达到0.00001时，简单迭代法迭代次数为6。 精度要求越高，简单迭代法迭代次数越多。

简单迭代法归结于找直线和曲线的交点的横坐标，其中迭代法的效果并不是总能令人满意的。例如：该程序中简单迭代法另一种等价形式

x=x*x*x-1对此建立迭代公式： x(k+1)=x(k)*x(k)*x(k)-1

迭代初始值x0=1.5，则有x1=2.375, x2=12.39.继续迭代下去已经没有必要了，为 结果显然会越来越大，不可能区域某个极限。这种不收敛的迭代过程是发散的，一个发散的迭代过程，显然其结果是毫无价值的。

二分法和简单迭代法比较：

当精度都是0.0001，所求结果为：1.32474时，简单迭代法迭代次数比二分法少9次，当精度都是0.00001，所求结果都为：1.32472时，简单迭代法迭代次数比二分法少11次。

简单迭代法逼近速度比二分法好，计算的效率比二分法提高了很多精度要求越高，简单迭代法的优势也越明显。

二分法不能用于求偶数重根和复根，多用于为其他求根方法提供初始近似值。 简单迭代法的迭代过程可能是收敛的也可能是发散的，简单迭代法对迭代公式有要求，迭代公式必须是收敛的，即结果逐渐趋近于某个极限。

Aitken算法程序源代码：

#include #include double s(double t) {

return (t*t*t-1); }

using namespace std; int main() { int i; double x,x0,x1,x2,e; cout<<\"请输入迭代初始值x0\"<<\和控制精度e\"<>x0>>e; i=0;

while(fabs(x0*x0*x0-x0-1)>e) { i++; x1=s(x0); x2=s(x1); x0=x2-(x2-x1)*(x2-x1)/(x2-2*x1+pow((x1+1),1.0/3.0)); } x=x0;

cout<<\"近似根x=\"<调试过程，实验结果及分析：

计算 x*x*x-x-1=0 在[1,2]内的近似根

。

精度达到0.0001时，程序运行结果如下图：

当将x=x0;改成x=x1，输出结果近似根为1.32517，迭代次数为5，这与要求得到的试验结果1.32472有一些差距，这是有循环语句i++; x1=s(x0); x2=s(x1); x0=x2-(x2-x1)*(x2-x1)/(x2-2*x1+pow((x1+1),1.0/3.0));执行顺序决定的近似根式子应该是x=x0;

由运行结果看出：

精度达到0.0001时，Aitken迭代法迭代次数为5.

Aitken迭代法是将迭代值在迭代一次，此时对于发散的x(k+1)=x(k)*x(k)*x(k)-1迭代公式，经过Aitken迭代法处理后却获得了相当好的收敛性。

Aitken迭代法和二分法的比较：

当精度都是0.0001时，Aitken迭代法迭代次数为5.，比二分法少9次，而且Aitken迭代法所求出的近似根为：1，32472，更接近理想值，Aitken迭代法是将迭代值在迭代一次从迭代次数和所求根的差异上可以看出Aitken迭代法是一种加速迭代法。

牛顿迭代法程序源代码：

#include #include double s(double t) { return t-(t*t*t-t-1)/(3*t*t-1); }

using namespace std; int main() { int i; double x,x0,x1,x2,e; cout<<\"请输入迭代初始值x0\"<<\和控制精度e\"<>x0>>e; i=0; x1=x0;

while(fabs(s(x1)-x1)>e) { i++; x2=s(x1); x1=x2; }

x=(x1+x2)/2;

cout<<\"近似根x=\"<调试过程，实验结果及分析：

计算 x*x*x-x-1=0 在[1,2]内的近似根。 精度达到0.00001时，程序运行结果如下图：

由运行结果看出：

精度达到0。00001时，牛顿迭代法迭代次数为3

x2=s(x1);x1=x2;这两句顺序不能换过来，当换过来时，错误如下

c:\\program files\\microsoft visual studio\\myprojects\\erww\\q.cpp(26) : fatal error C1004: unexpected end of file found 执行 cl.exe 时出错.：

这是因为在执行循环体之前没有对x2赋值，系统无法识别赋值到x1的值是多少，因此也无法的到想要的结果。

从运行结果看出：牛顿迭代法有很好的收敛性。在精度要求相同时，牛顿迭代法迭代次数比二分法,简单迭代法，Aitken迭代法都要少，牛顿迭代法是一种更快的迭代法。

牛顿迭代法是将非线性方程逐步转化为某些线性方程来求解，这是牛顿迭代法的基本思想。 将牛顿迭代法应用于以下方程求根： 1 x=(2-e^x+x*x)/3 2 3*x*x-e^x=0

对于 x=(2-e^x+x*x)/3 应用牛顿迭代法程序运行结果如下图：

对于3*x*x-e^x=0 应用牛顿迭代法程序运行结果如下图：

弦截法程序源代码：

#include #include double s(double t) { return t*t*t-t-1; }

using namespace std; int main() { int i; double x,x0,x1,x2,e; cout<<\"请输入迭代初始值x0,x1\"<<\和控制精度e\"<>x0>>x1>>e; i=0;

while(fabs(x1-x0)>e)

{ i++;

x2=x1-s(x1)*(x1-x0)/(s(x1)-s(x0)); x0=x1; x1=x2; }

x=(x1+x2)/2;

cout<<\"近似根x=\"<调试过程，实验结果及分析：

对于x0=1,x1=2 利用弦截法程序运行结果如下图：

弦截法对精度0.00001，需要迭代7次才能得到近似根，收敛速度比牛顿法慢些，但它的计算量比牛顿迭代法少，特别是导数的计算比较复杂时，弦截法更现实了它的优越性。 对于 x=(2-e^x+x*x)/3 应用弦截运行结果如下图：

对于

3*x*x-e^x=0 应用牛顿迭代法程序运行结果如下图：

弦截法局部收敛的速度比牛顿法略慢一些，弦截法计算过程中没有导数的计算，弦截法的计算量比牛顿法少。

实验小结及体会：

试验中练习了方程求根的迭代法，在选择和使用各种迭代格式中，最重要的是判断其收敛性及了解其收敛速度。在实际计算中，不一定都要判断迭代格式的收敛性才予使用，也可以根据对问题解的了解及获得的迭代序列来判断是否收敛于所要的根，但要注意，对于局部收敛的迭代格式（如牛顿迭代法）迭代初值的选取非常重要。

尽管每种数值方法有它本身的局限性，但它们还是很有效的。因此应该了解每一种方法的来龙去脉，尤其是它的困难所在，并通过计算机实习来熟悉和掌权这些方法.。

第二次试验

（1）Lagrange插值公式算法流程图

开始 定义两个数组a[100], b[100] 插值节点控制数n n+1个插值点 f=0 for i 0 To n forj 0 To n 是 i=j 否 tmp=tmp*(x-a[j])/(a[i]-a[j]); f=f+tmp*b[i]; 输出f_x 开始 （2）Newton插值公式算法流程图

插值计算for i o to n-1 tmp=tmp*(x-a[i]); f=f+tmp*b[i+1]; 输入插值截点控制数n,插值序列，和要计算的函数点x 开始 设定初始值t=1,newton=f(x0) tmp=tmp*(x-a[i]);f=f+tmp*b[i+1]; 输出插值计算的函数值

结束 （三）分段插值公式算法流程图

开始 输入n来确定插值次数 根据区间[a,b]，对给定的节点划分为多个小区间，每个小区间n+1个节点 在每个小区间上构造插值多项式Pi（x） 判断x所属的区间 计算所求插值节点的函数值 输出所求函数值 结束

（四）Hermite公式流程图

开始 插值节点序列，对应的导数值 根据节点个数确定插值次数然后for i 0 to n a=(1-2*(x-xi))/(xi-xk)*Li(x)*Li(x)) b=(x-xi)* Li(x)*Li(x) 根据插值公式计算插值点的值for i 0 to n H(x)=yi((1-2（x-xi）*a）+mi*b*b) 输出插值点所对应的值 结束 编程计算f(x)的近似值 （一）Lagrange插值：

拉格朗日插值程序源代码

#include #include using namespace std; int main() {

int i,j,n;

double a[100],b[100];

cout<<\"请输入插值节点控制数:n\"<>n;

cout<<\"请输入n+1个插值点:\"<cin>>a[i]>>b[i]; }

double f,x,tmp;

cout<<\"请输入所求函数值对应的x点\"<>x; f=0;

for(i=0;i<=n;i++) {

tmp=1;

for(j=0;j<=n;j++) if(i==j)continue; else tmp=tmp*(x-a[j])/(a[i]-a[j]); f=f+tmp*b[i]; }

cout<<\"所求函数值为f_x=\"<（二）Newton插值： 牛顿插值法程序源代码

#include #include using namespace std; int main() {

int i,j,n;

double a[100],b[100];

cout<<\"请输入插值节点控制数：n\"<>n;

cout<<\"请输入n+1个插值点:\"<>a[i]>>b[i]; }

double f,x,tmp;

cout<<\"请输入所求函数值对应的x点\"<>x;

for(i=0;ii;j--) { b[j]=(b[j]-b[j-1])/(a[j]-a[j-1-i]); } tmp=1; f=b[0]; {

for(i=0;icout<<\"所求函数值为f_x=\"<return 0; }

（三）调试过程，实验结果及分析：

（1）利用拉格朗日插值 对于f(x)=e^x, x0=0,x1=0.5 利用线性插值计算f(0.25)

的近似值

运行结果如下图：

对于x0=0.5,x1=1.0 利用线性插值计算f(0.75)的近似值

运行结果如下图：

对于x0=0,x1=0.5，x2=2.0利用二次插值计算f(0.25)的近似值 运行结果如下图：

对于x0=0,x1=0.5，x2=2.0利用二次插值计算f(0.75)的近似值 运行结果如下图：

（2）利用牛顿插值对于f(x)=e^x, x0=0,x1=0.5 利用线性插值计算f(0.25)的近似值

运行结果如下图：

对于x0=0.5,x1=1 利用线性插值计算f(0.75)的近似值 运行结果如下图：

对于x0=0,x1=0.5，x2=2.0利用二次插值计算f(0.25)的近似值 运行结果如下图：

对于x0=0,x1=0.5，x2=2.0利用二次插值计算f(0.75)的近似值 运行结果如下图：

（3）运行结果的比较：

①拉格朗日插值和牛顿插值计算出的结果相同，但是牛顿插值法的计算量较少。拉格朗日法形式对称，结构简单，便于在计算机上实现，只是计算量大。

②二次插值因为在取点更准确和计算复杂，所得到的结果精确度更好。

编程计算f(x)=

（一）分段插值：

的近似值

分段线性插值法程序源代码

#include

#include using namespace std; int main() {

int i,j,n;

double a[10],b[10],x,f,tmp;

cout<<\"请输入插值节点的个数:n\"<>n;

cout<<\"请输入n个插值节点:\"<>a[i]>>b[i]; }

cout<<\"请输入所求函数值对应的x点\"<>x; f=0;

if(xfor(i=1;iif(xloop:for(i=0;i<=1;i++) {

tmp=1;

for(j=0;j<=1;j++) if(i==j)continue; else tmp=tmp*(x-a[j])/(a[i]-a[j]); f=f+tmp*b[i]; }

cout<<\"所求函数值为\"<分段抛物线插值法程序源代码

#include #include using namespace std; int main() {

int i,j,n;

double a[10],b[10],x,f,tmp;

cout<<\"请输入插值节点的个数:n\"<>n;

cout<<\"请输入n个插值节点:\"<>a[i]>>b[i]; }

cout<<\"请输入所求函数值对应的x点\"<>x; f=0;

if(xfor(i=1;iif(xloop:for(i=0;i<=2;i++) {

tmp=1;

for(j=0;j<=2;j++) if(i==j)continue; else tmp=tmp*(x-a[j])/(a[i]-a[j]); f=f+tmp*b[i]; }

cout<<\"所求函数值为\"<调试过程，实验结果及分析：

1、分段线性插值程序运行结果如下：

2、分段抛物线插值程序运行结果如下：

3、由运行结果可以看出：

分段抛物线插值比分段线性插值更准确，这是因为抛物线插值在所求点的附近选取的点，并且线性插值和抛物线插值在所考虑区间内插值得到的数值精确度比在考虑区间外所得到的数值的精确度要高，这说明了内插比外推插值效果要好。

第二次实验小结及体会

1、插值法是数值分析中的一类重要方法，有着广泛的应用。本次试验中几种数值方法及有关概念和理论是数值分析最基本的内容

2、拉格朗日插值多项式是研究数值微积分与微分方程数值解的重要工具，但该方法计算量很大，且误差估计很困难，理论上还是比较重要。

3、牛顿插值多项式是拉格朗日插值多项式的变形，具有继承性，比拉格朗日插值多项式节省计算量，试验中用到了均差的概念，当节点数增加或减少时，只需将原有的计算结果增加或减少一项，着也是牛顿插值比拉格朗日插值方便的地方，应当注意。

4、分段低次多项式插值由于具有良好的稳定性与收敛性，且算法简单，可根据节点处函数值的具体情况，采用不同的插值公式，解决了曲线平滑性的问题，更便于应用。

第三次试验

（一）复化梯形求算法积流程图

开始 输入区间a,b.和区间等分数n 求出步长h For i=0 to n-1 计算节点函数值f(a+i*h); 计算积分值T=T+2*temp; Tn=(f(a)+f(b)+T)*h/2; 输出积分值Tn 结束

（三）变步长梯形求积算法流程图

（二）复化Simpson求积算法流程图

开始 输入区间a,b.和区间等分数n求出步长h For i=0 to n-1temp1=f(a+i*h); temp2=f(a+(i+0.5)*h); S1=S1+2*temp1; S2=S2+4*temp2; 计算积分值Sn=(f(b)-f(a)+S1+S2)*h/6; 输出积分值 结束

开始 输入区间a,b.和区间等分数n求出步长h . 循环运算For i=1 to n-1; temp1=f(a+i*h);T0=T0+2*temp1; 计算n等分后的积分值T[n-1]=(f(a)+f(b)+T0)*h/2;Tn=T[n-1输出积分值判断是否等分每个子区间后再求出积分值，是输入1，否输入0 否 Dir=1 是 等分后积分N++; n=2*n; h=(b-a)/n; temp=0;For i=1to n/2 temp2=f(a+(2*i-1)*h); temp=temp+temp2; T[n-1]=T[n/2-1]/2+temp*h; Tn=T[n-1] 输出经过N-1次变步长所得积分值Tn 判断是否等分每个子区间后再求出积分值，是输入1，否输入0 是 Dir=1 否 输出经过N-1次变步长所得积分值Tn 结束

（四）Romberg公式求积流程图

开始 输入区间a,b.和区间等分数n求出步长h . 计算Tn 利用Tn分别求出T2n, T4n,T8n Sn=(4*T2n-Tn)/3, S2n=(4*T4n-T2n)/3, S4n=(4*T8n-T4n)/3, Cn=(16*S2n-Sn)/15 C2n=(16*S4n-S2n)/15 Rn=(*C2n-Cn)/63 输出Rn 结束

编程计算积分

（一）复化梯形：

复化梯形求积程序

#include #include double f(double t) {

return 1/(1+t*t); }

using namespace std; int main()

{

int i,n;

double h,a,b,T,temp,Tn;

cout<<\"请输入求积区间:\"<>a>>b;

cout<<\"请输入区间等分数:\"<>n; h=(b-a)/n; temp=0; T=0;

for(i=1;itemp=f(a+i*h); T=T+2*temp; }

Tn=(f(a)+f(b)+T)*h/2;

cout<<\"经过复化梯形求积所得积分值T\"<调试过程，实验结果及分析：

1/（x*x）在区间[0,1]上8等分后,程序运行结果如下:

对区间等分后求每个子区间上的积分值，然后将每个子区间的积分值相加，就得到整个积分区间上的积分值，梯形公式具有一次代数精度。

（二）复化Simpson：

复化Simpson算法程序源代码

#include #include double f(double t) {

return 1/(1+t*t); }

using namespace std; int main() {

int i,n;

double h,a,b,S1,S2,temp1,temp2,Sn; cout<<\"请输入求积区间:\"<>a>>b;

cout<<\"请输入区间等分数:\"<>n; h=(b-a)/n; temp1=0; temp2=0; S1=0; S2=0;

for(i=0;itemp1=f(a+i*h); temp2=f(a+(i+0.5)*h); S1=S1+2*temp1; S2=S2+4*temp2; }

Sn=(f(b)-f(a)+S1+S2)*h/6;

cout<<\"经过复化Simpson求积所得积分值S\"<调试过程，实验结果及分析：

1/（x*x）在区间[0,1]上8等分后,程序运行结果如下:

复化Simpson求积公式具有三次代数精度，因此在区间和区间等分数相同时，所取得的积分值结果比梯形公式更加精确。和梯形公式一样都是利用线性组合来计算积分的近似值。

（三）变步长梯形求积：

变步长梯形求积程序源代码

#include #include double f(double t) {

return 1/(1+t*t); }

using namespace std; int main() {

int i,n,N,dir;

double h,a,b,T0,temp,temp1,temp2,temp3,T[180]; cout<<\"请输入求积区间:\"<>a>>b;

cout<<\"请输入区间等分数:\"<>n; h=(b-a)/n; temp=0; temp1=0; temp2=0; temp3=0; T0=0; N=1;

for(i=1;itemp1=f(a+i*h); T0=T0+2*temp1; }

T[n-1]=(f(a)+f(b)+T0)*h/2;

cout<<\"经过\"<1)

{ temp=0;

for(i=1;i<(n/2+1);i++) {

temp2=f(a+(2*i-1)*h); temp=temp+temp2; }

T[n-1]=T[n/2-1]/2+temp*h;

cout<<\"经过\"<cout<<\"是否需要继续二等分每个子区间,如果是，请输入:1 如果否，请输

入:0\"<cin>>dir; while(dir>0) {

N++; n=2*n; h=(b-a)/n; goto loop; }

if (dir<1)

cout<<\"经过\"<调试过程，实验结果及分析：

1/（x*x）在区间[0,1]上2等分后,程序运行结果如下

变步长梯形求积中对于n等分后Tn不满足精度要求可以等分每个子区间来求T2n的值此时只需计算新增节点的函数值，因而计算量变小了，程序运行速度也提高了。

（四）Romberg：

Romberg公式求积程序源代码

#include #include double f(double t) {

return 1/(1+t*t); }

using namespace std; int main() {

int i,j,n;

double h,a,b,T0,temp,temp1,temp2,temp3,T[180],S[180],C[180],R[180]; cout<<\"请输入求积区间:\"<>a>>b;

cout<<\"请输入区间等分数:\"<>n; h=(b-a)/n; temp=0; temp1=0;

temp2=0; temp3=0; T0=0;

for(i=1;itemp1=f(a+i*h); T0=T0+2*temp1; }

T[n-1]=(f(a)+f(b)+T0)*h/2; for(j=0;j<3;j++) {

temp=0; n=2*n; h=(b-a)/n;

for(i=1;i<(n/2+1);i++) {

temp2=f(a+(2*i-1)*h); temp=temp+temp2; }

T[n-1]=T[n/2-1]/2+temp*h; }

S[n/8-1] =(4*T[n/4-1]-T[n/8-1])/3; S[n/4-1] =(4*T[n/2-1]-T[n/4-1])/3; S[n/2-1] =(4*T[n-1]-T[n/2-1])/3; C[n/4-1]=(16*S[n/2-1]-S[n/4-1])/15; C[n/8-1]=(16*S[n/4-1]-S[n/8-1])/15; R[n/8-1]=(*C[n/4-1]-C[n/8-1])/63;

cout<<\"经过Romberg公式所得积分值R\"<调试过程，实验结果及分析：

1/（x*x）在区间[0,1]上2等分后,程序运行结果如下:

Romberg积分公式在积分结果中加入了事后误差估计值进行补偿，因而积分计算的收敛加速，并且仅采用梯形求积法不断折半步长，就可以逐步加工出Sn,S2n,S4n,Cn,C2n,Rn等精度较高的积分结果。

第三次试验小结

1、通过试验进一步学习了积分和微分的数值计算方法，其基本原理主要是逼近论，即设法构造某个简单函数近似表示，然后对求积或求导得到的积分或导数的近似值。基于差值原理，推导出数值积分和数值微分的基本公式。

2、插值型求积公式有公式和公式两类。前者取等距节点，算法简单而且容易编制程序。但是，由于其收敛性和稳定性都没有保证，所以，常用的低阶复化公式。

3、公式不但具有最高代数精度，而且收敛性和稳定性都有保证，因此是高精度的求积公式。公式还可以通过选择恰当的权函数，用于计算奇异积分和广义积分，也可以是一些复杂的积分计算简化。公式的主要缺点是节点与系数无规律。所以高阶公式不便于上机使用。实际应用中可以把低阶公式进行复化。

4、算法是在积分区间逐次二分的过程中，通过对梯形值进行外推加速处理，从而获得高精度的积分近似值。它具有自动选取步长的特点，便于在计算机上使用。外推的思想在数值微分，微分方程数值解等方法中也可以得到应用。

第四次试验

（一）三角分解法程序源代码

#include using namespace std; int main() {

const int MAX_N=20;

static double a[MAX_N][MAX_N],b[MAX_N],x[MAX_N],y[MAX_N]; static double u[MAX_N][MAX_N],l[MAX_N][MAX_N]; int i,j,k,n,r;

cout<<\"请输入矩阵阶数：\"; cin>>n;

if(n>MAX_N) {

cout<<\"The input n is larger than MAX_N,please redefine the MAX_N\"<if(n<=0) {

cout<<\"please input n between 1 and\"<//输入 a[i][j],b[i]

cout<<\"now input the matrix a[i][j],i,j=0...\"<for(j=0;j>a[i][j]; }

cout<<\"now input the matrix b[i],i=0...\"<>b[i]; //分解过程 for(i=0;il[i][i]=1;//l矩阵对角元素为1 for(k=0;kfor(j=0;j{ u[k][j]=a[k][j]; for(r=0;r<=k-1;r++) u[k][j]-=l[k][r]*u[r][j]; }

for(i=k+1;kfor(r=0;r<=k-1;r++) l[i][k]-=l[i][r]*u[r][k]; l[i][k]/=u[k][k]; } } //ly=b

for(i=0;iy[i]=b[i];

for(j=0;j<=i-1;j++) y[i]-=l[i][r]*y[j]; } //ux=y

for(i=n-1;i>=0;i--) {

x[i]=y[i];

for(j=i+1;jcout<<\"方程组的解为：\";//输出x for(i=0;i（二）追赶法程序源代码

#include using namespace std; int main() {

const int MAX_N=300;

int a[MAX_N],b[MAX_N],c[MAX_N],d[MAX_N]; float l[MAX_N],r[MAX_N],y[MAX_N],x[MAX_N]; int i,n;

cout<<\"请输入矩阵阶数：\"; cin>>n;

if(n>MAX_N) {

cout<<\"The input n is larger than MAX_N,please redefine the MAX_N\"<if(n<=0) {

cout<<\"please input n between 1 and\"<//输入 a[n],b[n],c[n],d[n]

cout<<\"now input the matrix a[n],n=2...\"<>a[i];

cout<<\"now input the matrix b[n],n=1...\"<>b[i];

cout<<\"now input the matrix c[n],n=1...\"<>c[i];

cout<<\"now input the matrix d[n],n=1...\"<>d[i]; //追的过程 r[0]=0; y[0]=0; a[1]=0; c[n]=0;

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

l[i]=b[i]-a[i]*r[i-1]; r[i]=c[i]/l[i];

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

//赶的过程 cout<for(i=n-1;i>=1;i--) {

x[i]=y[i]-r[i]*x[i+1];\\ }

cout<<\"方程组的解为：\"; for(i=1;i<=n;i++) {

cout<cout<计算得解为：－2 －3 －3 －3 －2

（三）雅可比迭代法和高斯赛德尔迭代法程序源代码

#include #include using namespace std; int a,b,m; double *x0;

void Jacobi(double **c,double *d,int n,double eps); void Gauss(double **c,double *d,int n,double eps); void main() { int n;

double **A,*B; double e;

cout<<\"请选择求方程组的迭代方法!雅可比选0，高斯-赛德尔选1!\"<>n;

cout<<\"输入方程组的变量的个数以及方程的个数!\"<>a>>b;

A=new double*[b]; for(int i=0;iA[i]=new double[a]; }

B=new double[b]; x0=new double[a];

cout<<\"输入每个方程组的变量的系数以及方程右端的值!\"<for(int j=0;jcin>>A[k][j]; }

cin>>B[k]; }

cout<<\"输入方程组迭代的次数及所要求的精度!\"<>m>>e;

cout<<\"输入方程组迭代的初值!\"<{

cin>>x0[j]; }

switch (n) {

case 0:Jacobi(A,B,m,e); break;

case 1:Gauss(A,B,m,e); break;

default:cout<<\"你没有选择求解方程组的一种方法!!\"<void Jacobi(double **c,double *d,int n,double eps) {

int k,i;

double *y = new double[a],*x=new double[a],s,temp=0.0; k=1; while(1) {

temp = 0.0; for(i=0;ifor(int j=0;jif(j!=i) {

s+=c[i][j]*x0[j]; } }

s=(d[i]-s)/c[i][i]; y[i]=s;

if(fabs(x0[i]-s)>temp) {

temp=fabs(x0[i]-s); } }

if(tempcout<<\"迭代成功!迭代结果为：\"<cout<<\"y[\"<break; }

if(k==m) {

cout<<\"迭代失败!!\"<for(i=0;ix0[i]=y[i]; } } }

void Gauss(double **c,double *d,int n,double eps) {

int k,i;

double *y=new double[a],*x=new double[a],s,temp=0.0; for(i=0;ix[i]=x0[i]; y[i]=x[i]; } k=1; while(1) {

temp=0.0;

for(i=0;ifor(int j=0;jif(j!=i) {

s+=c[i][j]*y[j]; } }

s=(d[i]-s)/c[i][i]; y[i]=s;

if(fabs(x[i]-s)>temp) {

temp=fabs(x[i]-s);

} }

if(temp{ cout<<\"迭代成功!迭代结果为：\"<cout<<\"y[\"<cout<<\"迭代失败!!\"<for(i=0;ix[i]=y[i]; } } }

用雅可比迭代法，程序运行结果如下:

用高斯赛德尔迭代法程序运行结果如下:

通过运行结果可以看出高斯赛德尔迭代法比雅可比迭代法所求结果精度高，收敛性好

。

（四）超松弛迭代法程序源代码

#include using namespace std; #include #define MAX_N 20 #define MAXREPT 100 #define epsilon 0.00001 int main() { int n; int i,j,k;

double err,w;

static double a[MAX_N][MAX_N],b[MAX_N][MAX_N],c[MAX_N],g[MAX_N]; static double x[MAX_N],nx[MAX_N];

cout<<\"input n value (dim of Ax=c):\";//输入方程的维数 cin>>n;

if(n>MAX_N) {

cout<<\"The input n is larger than MAX_N,please redefine the MAX_N\"<if(n<=0) {

cout<<\"please input n between 1 and\"<//输入a[i][j],c[i]

cout<<\"now input the matrix a[i][j],i,j=0...\"<for(j=0;j>a[i][j]; }

cout<<\"now input the matrix c[i],i=0...\"<>c[i];

cout<<\"now input the w value:\"); cin>>w;

if(w<=1||w>=2) {

cout<<\"w must between 1 and 2.\"<}

for(i=0;ib[i][j]=-a[i][j]/a[i][i];

g[i]=c[i]/a[i][i];//为了简化程序，假设a[i][i]!=0 //否则要附加对a[i][i]的处理 }

for(i=0;jfor(j=0;jcout<<\"after \"<＝1.2时 解为：0232179 0.0662066 0.494194

＝1.4时 解为：023218 0.0662066 0.494194 ＝1.6时 解为：023218 0.0662074 0.494193 由以上3组结果可知，越小，精度越大

实验小结及体会

1、任何非奇异矩阵都可分解为三角形式。

2、如果需要重复地求解系数矩阵相同，而右端常数项不同的线性方程组时，三角分解法是很有效的

3、当系数矩阵是对角占优矩阵时，用追赶法求线性代数方程组是非常简单的。 4、高斯迭代法公式比较复杂，但在一般情况下比雅可比迭代法收敛快。

5、高斯迭代法是异步迭代法每次迭代时都会用到刚才已经迭代出来的数值，因此迭代对收敛情况有所改善。

6、迭代法是否收敛，与迭代矩阵密切相关，在用迭代法进行迭代计算时应先判断系数矩阵的谱半径的绝对值与1的大小关系，从而可以决定迭代法的收敛性。 7、超松弛迭代法的迭代矩阵与松弛因子w有关。