计算方法实验指导书.

计 算 方 法

实 验 指 导 书

彭彬

计算机技术实验中心

2012年3月

· 实验环境： VC++ 6.0

· 实验要求：在机房做实验只是对准备好的实验方案进行验证，因此上机前要检查实验准备情况，通过检查后

方可上机。没有认真准备的学生不能上机，本次实验没有分数。实验中要注意考察和体会数值计算中出现的一些问题和现象：误差的估计,算法的稳定性、收敛性、收敛速度以及迭代初值对收敛的影响等。

· 关于计算精度：如果没有特别说明，在计算的过程中，小数点后保留5位数字，最后四舍五入到小数点后四

位数字。迭代运算的结束条件统一为

51

102

-?。在VC++ 6.0中，可使用setprecision 在流的输出中控制浮点数的显示（缺省显示6位）。演示如下： # include # include # include

//输出6位精度，输出左对齐

cout<

· 关于图形绘制

本课程个别实验要求画出函数的曲线，所有画图题目均要求用MFC 完成。

利用VC++6.0的MFC 画图，先要建立一个工程，然后在***View 中加入自定义变量、自定义函数等，最后在OnDraw （）方法中调用自定义函数。也可以把代码直接写入OnDraw （）方法中。

画曲线有两种方法，（一）一句坐标逐个打点（用SetPixel()函数），（二）先把当前光标移动（MoveTo

（）函数）到曲线的始点，再用LineTo （）函数画线。线的样式由画笔决定。对封闭区域可以填充，填充的样式由画刷决定。

在VC++6.0中,先新建一个MFC AppWizard(exe)类型的工程(建立工程时，“应用程序类型”选择“单文档”；“是否包含数据库”选择“不包含数据库”；其它选择缺省),然后在“ClassView ”中选择XXView 类文件加以操作。图1是一个名为test 的工程，在CtestView 节点，点击右键，用户可以增加自定义变量，自定义函数。

图1 MFC 工程

void CTestView::drawOldLine(int N, CDC *pDC) //N —画线所用点的数目

{

pDC->TextOut(250,10,"龙格现象");//文本输出 //设置画笔，将影响画线的样式 CPen *pnewPen,*poldPen; pnewPen=new CPen(); pnewPen->CreatePen(PS_SOLID,10,RGB(0,0,0)); poldPen=pDC->SelectObject(pnewPen);

//画坐标系

pDC->MoveTo(0,380); pDC->LineTo(640,380); pDC->MoveTo(320,0); pDC->LineTo(320,480); //画曲线,先移动到本曲线的第一个点 pDC->MoveTo(20,380); for(int i=0;i<=N;i++){ //原始坐标 float x0=-1+2.0/N*i;float y0=1.0/(1+4*x0*x0); //坐标放大 x0=x0*300;y0=-y0*300; //坐标转化为整型 int x1=(int)(320+x0);int y1=(int)(440+y0); pDC->LineTo(x1,y1); } //重置画笔 pDC->SelectObject(poldPen); }

图 2 自定义函数画函数2

1

()14f x x =+的图形

图3是一个根据Lagrange 插值多项式求函数值的自定义函数Lagrange(float x, int n, float x1[], float y1[])，其中n 是插值节点的个数，数组X1,y1存放插值节点的坐标，x 是待求点的x 坐标，函数根据插值多项式返回对应的y 坐标。

float CTestView::Lagrange(float x, int n, float x1[], float y1[]) {

float y=0;//存放函数值

int k=0;//控制变量,求lagrange 基函数的值 for( k=0;k

return y; }

图3 自定义函数，根据插值多项式求值

图4是onDraw 函数，这是MFC 画图的主要部分。相关代码如下： void CTestView::OnDraw(CDC* pDC) { CTestDoc* pDoc = GetDocument(); ASSERT_VALID(pDoc); // TODO: add draw code for native data here //画原始图形

float px[11],py[11];

图4 OnDraw()函数的内容

//N--N等分

for(int N=2;N<=10;N=N+2)

{

//求插值节点

for(int i=0;i<=N;i++)

{

px[i]=-1+2.0/N*i;

py[i]=1.0/(1+4*px[i]*px[i]);

}

//移动到第一个点

pDC->MoveTo(20,380);

//设置画笔

CPen *pnewPen,*poldPen;

pnewPen=new CPen();

pnewPen->CreatePen(PS_DASHDOT,1+N/2,RGB(255-20*N,0+10*N,20*N));

poldPen=pDC->SelectObject(pnewPen);

//按照插值节点求插值多项式的值,并描点

for(int k=0;k<=2000;k++)

{

float x=-1+2.0/2000*k;

float y=Lagrange(x,N+1,px,py);//注意N+1个节点

//坐标变换

x=x*300;

//Y坐标放大,并且反号,因为Y坐标轴是向下的

y=-y*300;

//转化为整型

int x1=(int)(320+x);

int y1=(int)(440+y);

if (N==2 && x1==80)

pDC->TextOut(x1,y1,"二等分");

if (N==4 && x1==50)

pDC->TextOut(x1,y1,"四等分");

if (N==6 && x1==60)

pDC->TextOut(x1,y1,"六等分");

if (N==8 && x1==580)

pDC->TextOut(x1,y1,"八等分");

if (N==10 && x1==600)

pDC->TextOut(x1,y1,"十等分");

pDC->LineTo(x1,y1);

}

//重置画笔

pDC->SelectObject(poldPen);

}

}

实验项目一览表

实验一 一元非线性方程求根的算法

【实验性质】

综合性实验。 【实验目的】

了解非线性方程求根的基本方法；掌握二分法和牛顿法。

【实验内容】

选择合适的初值，用二分法和牛顿迭代法求出一元非线性方程3

()sin 41f x x x x =--+的全部实根，

610ε-=。

【理论基础】

数值法求非线性方程的跟一般分三步： ① 判定根的存在性； ② 确定根的分布范围 ③ 根的精确化

常用的数值方法有二分法和迭代法。

二分法的基本思想是：先确定有根区间，然后通过逐步缩小（二分区间）有根区间的长度，以求根的近似值。二分法只能求实根，不能求复根及重根。

迭代法通过建立一个迭代序列来逼近根。迭代效果与迭代初值、迭代公式的收敛性有关。

迭代初始值一般根据高等数学的知识或作图的方法来确定；迭代公式通过对已知方程的某种变换得到。迭代的收敛性以及收敛速度不仅与迭代公式有关，可能还与迭代初值有关。

牛顿迭代是常用的迭代方法。把()f x 在k x 处按泰勒公式展开为：

()()k f x f x =+2()

()()()2!

...k k k f f x x x x x ξ'''-+

-+，

忽略高次项，得到()()()()k k k f x f x f x x x '≈+-。如果x *

是方程的根，则有：

由此得：

从而建立起牛顿迭代公式：

当相邻两次迭代结果满足1k k x x ε+-≤时停止计算。

如果非线性方程为一元多项式方程，则其牛顿解法的关键在于求()k f x 与()k f x '，这时可以用如下的快速算法（秦九韶算法）。

设n 次多项式()f x 为：

1n -次多项式()Q x 为： 如果()f x 、()Q x 满足

()()()k n f x x x Q x b =-+ （1）

则有

于是求()k f x 就转换为求系数n b 。

将式（1）展开，有如下递推关系：

因此可以用上述递推公式求出n b ，从而计算出()k f x 。 同理，设2n -次多项式 满足

1()()()k n Q x x x H x c -=-+ （2）

则有

由（2）式有递推关系：

因此求()k f x '可以用上式递推计算。 【实验过程】

1. 手工或用MFC 画出3

()sin 41f x x x x =--+的图形。 【MFC 画图代码】

2.估计三个有根区间，用二分法求出三根，并用下述样式的表描述求根过程。对每个求根过程，如果二分区间次数超过10次，则只记录最前面和最后面各5次。

【实验代码】

3.取三个初值，用牛顿迭代法求根，并用下述样式的表描述求根过程。对每个求根过程，如果迭代次数超过30次，则只记录最前面和最后面各15次

【实验代码】

4.比较两种方法，分析实验出现的问题，并总结解决办法。【实验心得】

【实验性质】

验证性实验。 【实验目的】

掌握Gauss 消元法和追赶法求解线性方程组。 【实验内容】

应用列主元消去法和追赶法迭代法求解下方程组： 【理论基础】

线性方程组的数值解法分直接算法和迭代算法。

直接算法是经过有限次的运算即求得（如果没有舍入误差）方程组精确解的方法，比如高斯消元法、矩阵分解法、追赶法等。在有限步内，迭代法得不到精确解。

直接法程序复杂，但运算量较小。由于受到计算机字长的限制，直接法也得不到精确解，仅适用于系数矩阵阶数不太高的线性方程组。 1.Gauss 算法：

Gauss 消元法分消元和回代两个过程，消元和回代过程都要求主元非零。如果经过k-1步后的主元10k kk a -≠，则可进行下一步：

(1)

(1)k ik

ik k kk

a l a

--=

，(1,2,...,i K k n =++)

()(1)(1)k k k ij ij ik ik a a a l --=-?（,1,2,...,i j k k n =++）

()(1)(1)k k k i i i ik b b b l --=-?（1,2,...,i k k n =++）

经过n-1步后，消元结束，然后回代：

本算法最的特点是顺序消元，又称顺序高斯消元法，与线性代数的处理方法有显著的不同，请同学们牢记。计算过程中，如果出现主元为0，将使Gauss 算法的计算终止。

在使用Gauss 消元法时，可以先判断系数矩阵是否为严格对角占优（即

1n

ii ij

j j i

a a =≠>∑(1,2,,)i n =），如果是就开始计算，否则就输出“不能保证主元非零，计算结束”。

在线性代数中，只要系数行列式det()0A A =≠，线性方程组就可求解，但Gauss 方法当主元为0

时就不能继续进行，因此顺序高斯消元法有缺陷。另一方面，即使顺序高斯消元法可行，但主元(1)

k kk a -很

小，运算中用它作为除法的分母（/ij ij ik kj kk a a a a a =-?），也会导致其它元素数量级的严重增长和舍入误差的扩散。这是Gauss 法的另一缺陷。

主元素消去法是对上述算法的改进，其中列主元消去法较为常用：在第i 步，从第i 列的第i 行至最后一行中选取绝对值最大的元素（注意，不是该列的所有元素，而是部分元素），通过行互换，将其调到主元位置，然后再做消元。用列选主元方法可以克服高斯消元法的额外限制，只要方程组有解，列选主元消元法就能畅通无阻地顺利求解，同时提高了解的精确度。 2、追赶法

在实践中，常遇到(3,5)m m =如对角方程组的求解问题。这类方程组的系数矩阵一般非奇异，可分解为特殊的矩阵，基于相关矩阵可得到求解的快速算法。 【实验过程】

1.画出选列主元的算法框图。

2.分别用高斯消元法和列主元消去法求解，用表格比较两种算法的结果与精度，分析实验出现的问题，并总结解决办法。

3.用追赶法求解，分别给出两个方程组Ly f =和Ux y =，分析实验出现的问题，并总结解决办法。

4.从使用条件、运算量比较列主元消去法和追赶法。 【实验心得】

【实验性质】

验证性实验。 【实验目的】

掌握迭代法求解线性方程组。 【实验内容】

应用雅可比迭代法和Gauss-Sediel 迭代法求解下方程组： 【理论基础】

线性方程组的数值解法分直接算法和迭代算法。迭代法将方程组的求解转化为构造一个向量序列，其极限就是方程组的解。迭代法程序简单，但有时工作量较大，在有限步内，得不到精确解，适宜某些系数矩阵阶数较高的问题。

迭代的基本思想是构造一个关于解向量的迭代序列{}

()K X ，使得这个序列随k 的增大而逐渐逼近准确解*X ，常见的有Jacobi 迭代算法和Gauss-Seidel 加速迭代算法。

设1112121

22

21

,1

000

000n n n n n nn a a a a a a A a a a -=++??????

??????

?

??????????????????

?

?

???＝L+D+U Jacobi 迭代公式：()

(1)

1()/(1,2,...,)n k k i ij j

i ii

j j i

x a x b a i n -=≠=-+=∑

Gauss-Seidel 迭代公式为：1

()()(1)1

1

1[()]/(1,2,...,)i n

k k k i

i ij j

ij j ii

j j i ii

x b a x

a x a i n a --==+=

-+=∑∑

迭代收敛的充分条件是：系数矩阵A 严格对角占优，或者迭代矩阵G 满足11G <或者1G ∞

<

本题的参考答案：0.2857、 -0.3265、0.0408、2.9796 【实验过程】

1.取三个不同初值，分别用雅可比法和高斯－塞得尔法求解，用表格记载求解过程。

【实验心得】

实验四 插值多项式的求法

【实验性质】

验证性实验。 【实验目的】

掌握Lagrange 插值算法、Newton 插值算法；理解Newton 插值算法相对于Lagrange 插值算法的优点。 【实验内容】

先用C 语言自带的系统函数sin x 求出sin x 在0,

,,,6432

x ππππ

=的值，然后分别用Lagrange 、Newton 方法求出sin()5π的值，并与用C 语言sin x 函数计算出的sin()5

π

作比较。

【理论基础】

https://www.360docs.net/doc/485428174.html,grange 插值公式

经过1n +个插值节点0011(,),(,),(,)n n x y x y x y 的n 次Lagrange 插值多项式为： 2.NewTon 插值公式

经过1n +个插值节点0011(,),(,),

(,)n n x y x y x y 的n 次Newton 插值多项式为：

上式中插商[,,

,]f x x x 的计算常用到如下的Newton 差商表：

由上表可知：每增加一个插值节点，都需从一阶插商算起。但在具体应用时，有两种方案： （一）当节点数n +1固定时，使用两个数组[1],[1]x n y n ++来存放插值计算的中间结果。其中[1]x n +存放插值节点的X 坐标；[1]y n +开始存放插值节点的Y 坐标，计算过程中存放差商值。

计算差商时，一次计算差商表的一列。计算01[,]f x x 时用到y[1]和y[0]，计算12[,]f x x 时用到y[2]和y[1]，不必再用到y[0]，由此可把01[,]f x x 放到y[0]，把12[,]f x x 放到y[1]，…。计算二阶差商时，把二阶差商放到y[0]、y[1]、…。这样牛顿插值公式中用到的f[x0]、f[x0,x1]、f[x0,x1,x2]、……,依次均存放在Y[0]上。

算法实现参考： Y ＝y[0]，t=1

for (int j=1;j<=n;j++) /计算J 阶插商 {

t=t*(x-x[j-1]);

for(int i=0;i

y[i]=[1][]

[][]

y i y i x i j x i +-+-;//计算j 阶插商的各种结果

}

Y=Y+y[0]*t;

}

f[x0,x1,x2]。当增加第n+1个节点时，用x[n+1]存放该点X 坐标，用y[n+1][0]存放点的Y 坐标，用y[n+1][i]存放本行的i 阶差商（1≤i ≤n ）。Y[n][i]＝(Y 的本行前1列―Y 的前1行前1列)/(X 的本行

1.分别用Lagrange 插值和Newton 插值求解。

2.用表格Lagrange 法、Newton 法以及C 语言自带函数对sin()5

的求值结果。

3.分析实验出现的问题，总结解决办法。 【实验心得】

实验五 数值积分

【实验性质】

验证性实验。 【实验目的】

理解插值型积分法；掌握复化积分法算法。 【实验内容】

对1

20

1

1I dx x =+?

，用复化梯形积分和变步长梯形积分求值（截断误差不超过5

1102

-?）。 【理论基础】

积分在工程中有重要的应用，数值积分的基本思想是用被积函数()f x 在区间[,]a b 上的一些点

(0,1,

,)i x i n =处的值()i f x 的线性组合作为积分的近似值：

实际应用中()f x 是未知的，一般用()f x 的次数不超过n 的插值多项式()n p x 来代替（插值型求积方

法）：

()()(())()n

b

b

b

n i i a

a

a

i f x dx p x dx l x dx f x =≈≈∑?

??。

如果考虑等距节点，即()/h b a n =-，(0,1,

,)i x a ih i n =+=，作变换x a th =+，有：

()()()n

b

k k a

k f x dx b a C f x =≈-∑?

，

上式称为Newton-Cotes 求积公式，其中0

0(1)

()!()!

n k

n

n

k j j k

C t j dt nk n k -=≠-=

--∏?称为Cotes 系数： 1n =:Cotes 系数为11

,

22 ； 2n =:Cotes 系数为141

,,666

3n =:Cotes 系数为1331

,,,8888

4n =: Cotes 系数为

73212327,,,,9090909090

1.Newton-Cotes 求积公式

Newton-Cotes 求积公式中，令1n =有梯形求积公式：[()()]2

b a

T f a f b -=+； Newton-Cotes 求积公式中，令2n =有抛物线公式（Simpson 公式）：

[()4(

)()]6

2

b a a b S f a f f b -+=

++；

Newton-Cotes 求积公式中，令4n =有Cotest 公式：

01234[7()32()12()32()7()]90

b a C f x f x f x f x f x -=

++++。

上述公式的余项分别为：

3

()(),[,]12

b a f a b ηη-''-∈、5

(4)

4

()(),[,]1802

b a f

a b ηη--

∈?和7

(6)

5

()(),[,]9452

b a f

a b ηη--

∈?

2.复化求积公式：

为提高精度，常把区间[,]a b 等分成n 个小区间，在每个小区间上使用梯形公式、Simpson 公式、Cotes 公式，就得到复化积分公式:

1

1

[()()2()]2n n k h T f a f b f a kh -==+++∑ ―复化梯形公式()b a h n -=

11

1/20

1

[()()4()2()]6

n n n k k k k h S f a f b f x f x --+===

+++∑∑ ―复化Simpson 公式()b a h n -=

―复化Cotes 公式()b a

h n

-=

为了能够控制计算精度，实践中常用变步长（每次区间减半）的梯形积分（辛普森积分、柯特斯积分）；它是充分利用复化积分中2n T 与n T 的关系得到的：

如果计算过程中，2n n T T ε-≤，就停止计算，用2n T 作为计算的结果。 对于复化Simpson 、复化Cotes 积分，用同样的方法处理。 【实验过程】

1.用复化梯形积分求解，用表格记载求解过程。

2.用变步长梯形积分求解，用表格记载求解过程。

3.比较复化积分与变步长积分，分析实验出现的问题，总结解决办法。【实验心得】

实验六 常微分方程组的数值解

【实验性质】

验证性实验。 【实验目的】

理解欧拉法的特点、几何意义；掌握Runge-Kutta 法的特点、算法。 【实验内容】

分别取h=0.10、0.02，用Euler 法和四阶Runge-Kutta 法求解初值问题 要求：

（1）要求计算结果小数点后保留6位有效数字。

（2）用表格比较数值解与精确解()cos sin y x x x =+。 （3）在同一坐标系画出数值解曲线与精确解曲线。 【理论基础】

1.Euler 公式

2.四阶Runge-Kutta 公式

00

11234()(22)6n n y x y h

y y k k k k +==++++?????，11132

43

(,)11(,)22

11(,)22(,)n n n n n n n n k f x y k f x h y hk k f x h y hk k f x h y hk ==++=++=++????????? 【实验过程】

1.分别用欧拉法和四阶Runge-Kutta 法求解。

2.用表格记载程序运行的过程。

3.比较两种方法的优劣；分析实验出现的问题，总结相关解决办法。【实验心得】