}
例2．已知函数表，用拉格朗日插值多项式求0.5，0.7，0.85三点处的函数值。
程序的结构：
在程序的说明部分将已知数据分别存入一维数组x和y。
在主函数main中提示人工键入插值点的数目m，继而在m次循环计算中分别提示键入插值点，计算插值，并输出结果。
程序中的主要常量，变量及函数：
n为插值基点的最大下标。
x[n+1],y[n+1]分别存放插值基点及其函数值。
M为插值点个数，X存放插值点的值。
程序清单：
#include"stdafx.h"
#include
#definen4/*插值基点的最大下标*/
voidmain（）
{
doublex[n+1]={0.4,0.55,0.8,0.9,1};/*插值基点的值*/
doubley[n+1]={0.41075,0.57815,0.88811,1.02652,1.17520};
intm,k,i,j;
floatX,L,P;
printf（"\nPleaseenterm="）;
scanf（"%d",&m）;/*键入插值点的个数*/
for（k=1;k<=m;k++）
{
printf（"\nPleaseenterX%d=",k）;
scanf（"%f",&X）;
P=0.0;
for（i=0;i<=n;i++）
{
L=1.0;
for（j=0;j<=n;j++）
{
if（j!
=i）
{
L=L*（X-x[j]）/（x[i]-x[j]）;
}
}
P=P+y[i]*L;
}
printf（"P（%f）=%f\n",X,P）;
}
}
计算结果：
（当提示“Pleaseenterm=”时，键入插值点个数3，再分别根据提示X1=或X2=，X3=依次键入插值点0.5,0.7,0.85），输出的结果分别为
P（0.500000）=0.521090P（0.700000）=0.758589P（0.850000）=0.956119
§2牛顿插值
1．方法概要：
牛顿n次插值多项式
Pn（x）=Pn-1（x）+an（x-x0）（x-x1）…（x-xn-1）
在计算机上计算均差表时，可用一维数祖x存放基点值，用二维数组F存放和计算由个阶均差构成的下三角矩阵。
计算公式如下：
Fi0=f（xi）,i=0,1,…n
Fij=（Fi,j-1-Fi-1,j-1）/（xi-xi-j）,j=1,2,…ni=j,j+1,…n
2.程序及例
例3.已知函数表同例2，试根据上面公式的算法计算均差表存入矩阵F，并用牛顿均差插值多项式求解上题。
注意到
，那么
计算
的流程为：
（1）
。
（2）对
，计算
。
程序的结构：
在程序开始时将初始数据分别存入一维数组x,y。
主函数main首先调用自定义函数junca计算均差表，再提示人工键入插值点个数m，继而在m次循环中分别提示键入插值点X，并计算相应的插值，输出值结果。
程序中的主要常量，变量及函数：
n,x[],y[],m,X均同上例。
Junca（doublex[],doubley[],doubleF[][n+1]）是自定义函数，计算均差表并打印输出，结果存于F中，P用于计算牛顿均差插值多项式的值。
#include"stdafx.h"
#include
constintn=4;
/*函数junca实现均差表的计算，将均差表存入数组F，并在屏幕上打印输出*/
voidjunca（floatx[],floaty[],floatF[][n+1]）
{
inti,j;
for（i=0;i<=n;i++）
{
F[i][0]=y[i];
}
for（j=1;j<=n;j++）
{
for（i=j;i<=n;i++）
{
F[i][j]=（F[i][j-1]-F[i-1][j-1]）/（x[i]-x[i-j]）;
}
}
printf（"\n%12s%12s\n","Xi","F（Xi）"）;
for（j=1;j<=38;j++）
printf（"--"）;
printf（"\n"）;
for（i=0;i<=n;i++）
{
printf（"%12f",x[i]）;
for（j=0;j<=i;j++）
printf（"%12f",F[i][j]）;
printf（"\n"）;
}
for（j=1;j<=38;j++）
printf（"--"）;
printf（"\n"）;
}
voidmain（）
{
floatx[n+1]={0.4,0.55,0.8,0.9,1};/*插值基点的值*/
floaty[n+1]={0.41075,0.57815,0.88811,1.02652,1.17520};
inti,m,k;
floatp,t,F[n+1][n+1];
floatX;
junca（x,y,F）;/*调用函数junca计算均差表*/
printf（"\nPleaseenterm="）;
scanf（"%d",&m）;/*键入插值点的个数*/
for（i=1;i<=m;i++）
{
scanf（"%f",&X）;/*键入值点*/
p=F[n][n];
for（k=n-1;k>=0;k--）
p=p*（X-x[k]）+F[k][k];
/*计算牛顿均差值多项式*/
printf（"P（%f）=%f\n",X,p）;
}
}
计算结果：
Xi
F（Xi）
1阶
2阶
3阶
4阶
0.400000
0.410750
0.550000
0.578150
1.116000
0.800000
0.888110
1.239840
0.309600
0.900000
1.026520
1.384100
0.412171
0.205143
1.000000
1.175200
1.486800
0.513500
0.225175
0.033386
（当提示Pleaseenterm=时，键入插值点个数3，再分别根据提示X1=或X2，X3=依次键入插值点0.5,0.7,0.85）输出的结果为
P[0.500000]=0.521090,P[0.700000]=0.758589,P[0.850000]=0.956119
练习题
1．已知函数表：
xi0.40.550.650.80.9
yi0.410750.578150.696750.888111.02652
分别用四次拉格朗日插值和牛顿插值，求x=0.596时的函数值。
第二章线性方程组的解法
本章的目的与要求：
1．熟悉求解线性方程组的有关理论与方法。
2．会编制列主元消去法、LU分解法与两种迭代法的计算程序，上机计算求得结果，
3．通过上机实习对计算数据的分析，进一步了解各种方法的功能，优缺点与适用范围，体会如何针对不同问题适当选择方法。
§1高斯消去法
1.方法概要
解线性方程组
（1）
为统一计算公式与简化计算程序，把方程组的增广矩阵记作A
（1）=（A,b）
（2）
解题分为消元与回代两个过程：
消元过程；
第1步，对第1列作消元，将
消为零。
第2步，对第2列作消元，将
消为零。
设
≠0（k=1,2,3,4….,n-1）。
一般地,对于第k步消元,采用如下公式:
（3）
（4）
（j=k,k+1,…,n+1）
直到作完ｎ-1步消元，相应的方程组化为与原来方程组等价的方程组：
（6）
回代过程：
回代过程是解方程组
，采用如下公式
（7）
2．程序流程图（见图2—1）
流程图说明
（1）首先输入方程组的阶数n及方程组的增广矩阵Ａ。
Ａ为一个ｎ×（ｎ+1）的二维数组：
所有运算的中间结果与最后结果全部放Ａ中相应的位置。
（2）本程序为了教学的需要，先屏幕输出原始增广矩阵，然后逐次输出每次消元后的矩阵以及每次所用的消元因子，最后输出解向量。
我们可以从中看出消元过程的具体情况。
在实用中无此必要，可以省略。
（3）这是顺序取主元的消去法。
为了检查
是否为零，设置了一个检查口，一旦遇到
，便使其输出一个记号，以示顺序取主元的方法失败.
3．程序及例
例1.求解方程组
2.5x１+2.3x２-5.1x3=3.7
5.3x１+9.6x２+1.5x３=3.8
8.lx１+1.7x2-4.3x３=5.5
（答：
x=（0.406705,0.237818,-0.419325）T）
#include"stdafx.h"
#include
#include
#definen3/*方程组的阶数*/
doubleaa[n][n+1]={{2.5,2.3,-5.1,3.7},{5.3,9.6,1.5,3.8},{8.1,1.7,-4.3,5.5}};
/*增广矩阵的初始数据*/
intgauss（doublea[][n+2],doublex[]）
{
inti,j,k;
floatc;
for（k=1;k<=n-1;k++）/*消元过程*/
{
for（i=k+1;i<=n;i++）
{
if（fabs（a[k][k]）<1e-12）
{
printf（"\ndet=0.fail!
\n"）;
return（0）;
}
c=a[i][k]/a[k][k];
for（j=k+1;j<=n+1;j++）
{
a[i][j]=a[i][j]-c*a[k][j];
}
}
}
for（k=n;k>=1;k--）/*回代过程*/
{
x[k]=a[k][n+1];
for（j=k+1;j<=n;j++）
{
x[k]=x[k]-a[k][j]*x[j];
}
x[k]=x[k]/a[k][k];
}
return
（1）;
}
voidmain（）
{
inti,j,det;
doublea[n+1][n+2],x[n+1];
for（i=1;i<=n;i++）
for（j=1;j<=n+1;j++）/*用a[1][1]~a[n][n+1]存放增广矩阵*/
a[i][j]=aa[i-1][j-1];
det=gauss（a,x）;/*调用函数gauss求解方程组，并获取返回标志值*/
if（det!
=0）
{
for（i=1;i<=n;i++）
printf（"\nx[%d]=%f",i,x[i]）;
printf（"\n--------------------------\n"）;
}
}
是
图2－1Gauss消去法程序流程图
§２ 列主元消去法
1．方法概要
作第K步消元时，从A（K）的子矩阵（见§1中式（3））的第一列中选取
｜,以
作为主元素，第
行称为主行。
将主行与第k行交换，然后再作这一步的消元。
2．计算框图
图2－2列主元消去法框图
3．程序及例
例2.用列主元消去法求解
程序的结构：
主函数main调用函数gaussl进行列主元高斯消去法计算，根据返回的函数值为1或0而输出计算结果或失败信息。
程序中的主要常量、变量及函数：
n为常量，指定方程组的阶数。
数组aa[n][n+1]存放方程组增广矩阵的原始数据并保留不变。
数组a[n+1][n+2]、x[n+1]分别用于计算增广和方程组的解。
intgauss1（a,x）是自定义函数，用列主元高斯消去法求解n阶线性方程，当系数矩阵奇异使求解失败时，输出信息“det=0.fail!
”并返回函数值0；否则求解成功，返回函数值1。
Intdet在主函数里接受gaussl的返回值，当det不为0时主函数打印出方程组的解。
#include"stdafx.h"
#include
#include
#definen3/*方程组的阶数*/
staticdoubleaa[n][n+1]={{1,2,-1,3},{1,-1,5,0},{4,1,-2,2}};/*增广矩阵的初始数据*/
intgauss1（doublea[][n+2],doublex[]）
{
inti,j,k,r;
doublec;
for（k=1;k<=n-1;k++）/*消元过程*/
{
r=k;
for（i=k;i<=n;i++）/*选列主元*/
{
if（fabs（a[i][k]）>fabs（a[r][k]））
{
r=i;
}
if（fabs（a[r][k]）<1e-12）
{
printf（"\ndet=0.fail!
\n"）;
return（0）;
}
if（r!
=k）
for（j=k;j<=n+1;j++）/*交换k、r两行*/
{
c=a[k][j];
a[k][j]=a[r][j];
a[r][j]=c;
}
}
for（i=k+1;i<=n;i++）/*进行消元计算*/
{
c=a[i][k]/a[k][k];
for（j=k+1;j<=n+1;j++）
{
a[i][j]=a[i][j]-c*a[k][j];
}
}
}
if（fabs（a[n][n]）<1e-12）
{
printf（"|ndet=0.fail!
\n"）;
return（0）;
}
for（k=n;k>=1;k--）/*回代过程*/
{
x[k]=a[k][n+1];
for（j=k+1;j<=n;j++）
{
x[k]=x[k]-a[k][j]*x[j];
}
x[k]=x[k]/a[k][k];
}
return
（1）;
}
voidmain（）
{
inti,j,det;
doublea[n+1][n+2],x[n+1];
for（i=1;i<=n;i++）
for（j=1;j<=n+1;j++）/*用a[1][1]~a[n][n+1]存放增广矩阵*/
a[i][j]=aa[i-1][j-1];
det=gauss1（a,x）;/*调用函数gauss1求解方程组，并获取返回标志值*/
if（det!
=0）
{
for（i=1;i<=n;i++）
{
printf（"\nx[%d]=%f",i,x[i]）;
}
printf（"\n--------------------------\n"）;
}
}
计算结果：
x1=0.250000x2=1.500000x3=0.250000
§3 线性方程组的迭代解法
1.雅可比迭代法
例3.根据雅可比迭代法编制程序求解方程组
。
其中
的计算流程为：
（1）
;
（2）对
做：
①
②若
则
。
程序中的主要常量，变量：
n，eps是常量，表示方程组的阶数及指定的精度。
aa[n][n],bb[n]分别存放方程组的系数矩阵及常数项。
intNO是最大迭代次数，当迭代达到NO次以上时程序输出失败信息，并结束。
x[n+1],y[n+1]分别用于存放相继两次的迭代向量。
用变量max计算。
程序清单:
//Definestheentrypointfortheconsoleapplication.
//
#include"stdafx.h"
#include
#include
#definen3/*方程组的阶数*/
#defineeps0.5e-4/*给定精度要求*/
staticdoubleaa[n][n]={{10,-1,-2},{-1,10,-2},{-1,-1,5}};
staticdoublebb[n]={7.2,8.3,4.2};
voidmain（）
{
intk=0,NO,i,j;
doubled,s,max;
doublea[n+1][n+1],b[n+1],x[n+1],y[n+1];
for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
a[i][j]=aa[i-1][j-1];
b[i]=bb[i-1];
}
printf（"\nPleaseenterNo:
"）;
scanf（"%d",&NO）;/*键入最大迭代次数NO*/
for（i=1;i<=n;i++）
y[i]=0;
do{
k++;
for（i=1;i<=n;i++）
x[i]=y[i];
if（k>=NO）
{
printf（"\nFail!
"）;
break;
}/*当k>=NO时输出失败信息，结束迭代*/
max=0.0;
for（i=1;i<=n;i++）
{
s=0.0;
for（j=1;j<=n;j++）
if（j!
=i）
{
s=s+a[i][j]*x[j];
}
y[i]=（b[i]-s）/a[i][i];
d=fabs（y[i]-x[i]）;
if（max}
}while（max>=eps）;
if（max{
printf（"k=%d\n",k）;
for（i=1;i<=n;i++）
printf（"x[%d]=%f\n",i,x[i]）;
}
}
运行结果：
产生输入提示PleaseenterNO：
后，键入适当大的正整数，可行到结果数据为
k=11
x[1]=1.099979
x[2]=1.199979
x[3]=1.299975
2.Gauss-Seidel迭代法
例4.根据高斯-赛德尔迭代法编制程序求解上题的方程组。
程序中的主要常量及变量：
n,eps,aa[][],bb[],a[][],b[],x[]，NO均同上例。
每次迭代计算时，用变量s临时存放x[i]的旧值。
用d计算x[i]与其旧值之差最大绝对值。
程序清单:
//Definestheentrypointfortheconsoleapplication.
#include"stdafx.h"
#include
#include
#definen3
#defineeps0.5e-4
staticdoubleaa[n][n]={{10,-1,-2},{-1,10,-2},{-1,-1,5}};
staticdoublebb[n]={7.2,8.3,4.2};
voidmain（）
{
intk=0,NO,i,j;
doubled,sum,s;
doublea[n+1][n+1],b[n+1],x[n+1];
for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
a[i][j]=aa[i-1][j-1];
b[i]=bb[i-1];
}
printf（"\nPleaseenterNO:
"）;
scanf（"%d",&NO）;
for（i=1;i<=n;i++）
{
x[i]=0;
}
do{
k++;
d=0.0;
if（k>=NO）
{
printf（"\nFail!
"）;
break;
}
for（i=1;i<=n;i++）
{
sum=0.0;
s=x[i];/*s临时存放x[i]的旧值*/
for（j=1;j<=n;j++）
if（j!
=i）
{
sum=sum+a[i][j]*x[j];
}
x[i]=（b[i]-sum）/a[i][i];/*计算x[i]的新值*
升级会员 