实验一 --插值方法

实验一—插值方法

实验学时：4

实验类型：设计 实验要求：必修

一 实验目的

通过本次上机实习，能够进一步加深对各种插值算法的理解；学会使用用三种类型的插值函数的数学模型、基本算法，结合相应软件（如VC/VB/Delphi/Matlab/JAVA/Turbo C ）编程实现数值方法的求解。并用该软件的绘图功能来显示插值函数，使其计算结果更加直观和形象化。

二 实验内容

通过程序求出插值函数的表达式是比较麻烦的，常用的方法是描出插值曲线上尽量密集的有限个采样点，并用这有限个采样点的连线，即折线，近似插值曲线。取点越密集，所得折线就越逼近理论上的插值曲线。本实验中将所取的点的横坐标存放于动态数组[]X n 中，通过插值方法计算得到的对应纵坐标存放

于动态数组[]Y n 中。

以Visual C++.Net 2005为例。

本实验将Lagrange 插值、Newton 插值和三次样条插值实现为一个C++类CInterpolation ，并在Button 单击事件中调用该类相应函数，得出插值结果并画出图像。CInterpolation 类为 class CInterpolation { public :

CInterpolation();//构造函数

CInterpolation(float *x1, float *y1, int n1);//结点横坐标、纵坐标、下标上限 ~ CInterpolation();//析构函数 ………… …………

int n, N;//结点下标上限，采样点下标上限

float *x, *y, *X;//分别存放结点横坐标、结点纵坐标、采样点横坐标

float *p_H,*p_Alpha,*p_Beta,*p_a,*p_b,*p_c,*p_d,*p_m;//样条插值用到的公有指针，分别存放

i h ，i α，i β，i a ，i b ，i c ，i d 和i m

};

其中，有参数的构造函数为

CInterpolation(float *x1, float *y1, int n1) {

//动态数组x1，y1中存放结点的横、纵坐标，n1是结点下标上限（即n1+1个结点） n=n1;

N=x1[n]-x1[0]; X=new float [N+1]; x=new float [n+1]; y=new float [n+1]; for (int i=0;i<=n;i++) {

x[i]=x1[i]; y[i]=y1[i]; }

for (int i=0;i<=N;i++) X[i]=x[0]+i; }

2.1 Lagrange 插值

()()n

n i i i P x y l x ==∑，其中0,()n

j i j j n

i j

x x l x x x =≠-=

-∏

对于一个自变量x ，要求插值函数值()n P x ，首先需要计算对应的Lagrange 插值基函数值()i l x float l(float xv,int i) //求插值基函数()i l x 的值 {

float t=1;

for (int j=0;j<=n;j++) if (j!=i)

t=t*(xv-x[j])/(x[i]-x[j]); return t; }

调用函数l(float x,int i)，可求出()n P x

float p_l(float x) //求()n P x 在一个点的插值结果 {

float t=0;

for (int i=0;i<=n;i++) t+=y[i]*l(x,i); return t; }

调用p_l(float x)可实现整个区间的插值

float *Lagrange() //求整个插值区间上所有采样点的插值结果 {

float *Y=new float [N+1]; for (int k=0;k<=N;k++) Y[k]=p_l(x[0]+k*h);

return Y; } 2.2

Newton 插值

010

()(,,)()n

n i i i P x f x x x x ω==∑ ，

其中1

01,0()(),0i i j j i x x x i ω-==??=?-≠??

∏，0100,()

(,,)()

i k i n

k k j j j k

f x f x x x x x ==≠=-∑∏

对于一个自变量x ，要求插值函数值()n P x ，首先需要计算出01(,,)i f x x x 和()i x ω

float *f() {

//该函数的返回值是一个长度为n ＋1的动态数组，存放各阶差商 }

float w(float x, int i) {

//该函数计算()i x ω }

在求()n P x 的函数中调用*f()得到各阶差商，然后在循环中调用w(float x)可得出插值结果 float p_n(float x) {

//该函数计算()n P x 在一点的值 }

调用p_n(float x)可实现整个区间的插值 float *Newton() {

//该函数计算出插值区间内所有点的值 }

2.3 三次样条插值

三次样条插值程序可分为以下四步编写： （1） 计算结点间的步长i hi 、i α、i β；

（2） 利用i hi 、i α、i β产生三对角方程组的系数矩阵和常数向量； （3） 通过求解三对角方程组，得出中间结点的导数i m ； （4）

对自变量x ，在对应区间1[,]i i x x +上，使用Hermite 插值；

（5）调用上述函数，实现样条插值。将每步写成函数：

(1)

void GetH(void)

{

hi

//该函数计算数组

i

}

void GetAlpha(void)

{

α

//该函数计算数组

i

}

void GetBeta(void)

{

β

//该函数计算数组

i

}

(2)

void Geta(void)

{

a

//该函数计算数组下对角线

i

}

void Getb(void)

{

b

//该函数计算数组主对角线

i

}

void Getc(void)

{

c

//该函数计算数组上对角线

i

}

void Getd(void)

{

d

//该函数计算方程组右端常数项

i

}

(3)

float *Chasing(float *pa,float *pb,float *pc,float *pd,int n) {

//追赶过程，计算各点斜率

i

m

}

(4)

float F0(float x)

{

//该函数计算函数

0()x

?的值}

float F1(float x)

{

//该函数计算函数

1()x

?的值}

float P0(float x)

{

//该函数计算函数

0()x

ψ的值}

float P1(float x)

{

//该函数计算函数

1()x

ψ的值

}

调用上述函数，实现三次样条插值

float *Cubic_Spline()

{

//该函数计算所有点的插值结果

}

三实验组织运行要求

实验前，由任课教师落实实验任务，每个学生事先编写好算法设计源程序代码。集中上机、调试并通过计算机图形可视化演示操作实例来测试、验证所学的数值分析理论。

四实验条件

为每个学生提供一台具有WINDOWS 98/XP/NT/2000操作系统的计算机；同时提供VC++/VB/JAVA/TC等集成语言开发环境来编程设计计算方法的上机实验。

五实验步骤

1．根据实验内容和算法流程图预先编好程序初稿，上机调试、运行，输出正确的结果。

2．三种类型的插值函数所生成的图形要显示在同一坐标轴上，并用三种不同的颜色分别表示出来。

3．分析它们的运行结果，并比较三种插值函数各自不同的特点。

4．实验完毕后提交实验报告

六实验报告

实验报告的内容包括：实验目的、实验内容、源程序清单和结果分析及心得体会。并分析它们的运行结果，并比较三种插值函数各自不同的特点。

上述是第一个实验——插值方法程序的构思和编写过程。将其简单概括，可分为以下几步：

（1）按照本实验指导书的开篇——“集成开发环境Visual C++.Net的使用”中的第一部分——“开发环境介绍”建立项目。

（2）为项目添加一个类，在该类中编写插值算法。该类需要一个有参数的构造函数将类外的插值条件数据传递到类中。

（3）按照“集成开发环境Visual C++.Net的使用”所述为项目添加对话框双击事件响应程序，将鼠标双击的坐标存放于数组或线性表中，用于插值，并将双击的点用小圆圈标出。

（4）在对话框上添加一个Button（按钮），双击改按钮，添加响应程序，调用拉格朗日插值函数，并将得到的相邻插值点用直线连接起来。连线的方法在“集成开发环境Visual C++.Net的使用”有介绍。

（5）用与（4）类似的方法添加牛顿插值和三次样条插值的Button和相应的响应程序。

对于插值结点的存储，可以采用线性表。为了通过插值结点确定插值区间，可将这些点按照横坐标递增排序。

务必熟练实验一程序编写的每一步骤，做到能够从头至尾独立编写。在项目的建立、画图、类的添加和编程思维各方面，后面的四个实验与实验一都是类似的。本实验指导书对后面的四个实验不再如此赘述。

实验二——数值积分

实验学时：2

实验类型：设计

实验要求：必修

一实验目的

通过该课程实习，学会使用数值积分的各种方法求解定积分计算的问题，体会各种方法的精度差异。二实验内容

本实验将梯形法的递推化和龙贝格算法编写为类CIntegration，在工程中调用该类得出积分结果并画出区间二分图像。首先画出被积函数图像，并储存被积函数在图像中每个像素对应的采样点的横坐标（存放于float型指针p_X中）和纵坐标（存放于int型指针p_Y中）；然后调用CIntegration类中相应函数，计算积分结果；最后将x轴上的区间二分点与其对应的被积函数值的点连线。

CIntegration有公有变量和指针float a,b,h,*x;int n和带参数的构造函数

CIntegration::CIntegration(float ma,float mb)

{

a=ma; //a、b分别为积分区间左右端点

b=mb;

n=1; //区间二分个数

h=b-a; //步长

x=new float [n+1]; //存放区间二分点 x[0]=a; x[1]=b;

}

另需被积函数

float CIntegration::f(float x) { float mTemp=sin(x); //可替换为其他函数 return mTemp;

}

绘制函数图像时，需要将图像按照比例放大。 2.1 梯形法的递推化

梯形法递推化公式为

1

2102

1()22n n n k k h T T f x -+==+∑

该过程是下列过程的往复循环 ①

求已有各区间中点函数值的和

1

1

2

()n k k f x

-+

=∑；

② 求积分值2n T 。

③

区间二分，即更新结点，并将步长h 减半。

要实现该算法过程，需要先编写下面两个函数 float CIntegration::Sum(void ) { //该函数求

1

10

2

()n k k f x

-+

=∑的值

}

void CIntegration::Divide(void ) {

//该函数实现区间二分和步长更新 float *px=new float[2*n+1]; for(int i=0;i

px[2*i+1]=(x[i]+x[i+1])/2;

}

px[2*n]=x[n]; delete x; x=px;

n=n*2;

h=h/2;

}

为了后面编程的方便，可将单次二分过程写为函数 float CIntegration::T1ToT2(float T1) { }

最后，调用T1ToT2，可实现梯形法递推化的函数 float CIntegration::Trapezium_Recurrence(float e) { }

在后面的Romberg 函数中，还需要用到迭代一定次数的梯形法递推函数 float CIntegration::Trapezium_Recurrence(int N) { }

2.2 龙贝格算法

其中

1

2102

1()22n n n k k h T T f x -+==+∑

243n n n T T S -=

21615n n

n S S C -=

26463n n

n C C R -=

调用梯形法递推函数Trapezium_Recurrence(int N)可以很方便的实现单次龙贝格算法的函数 float CIntegration::R1(float * T) { }

调用函数R1可实现迭代的龙贝格函数 float CIntegration::Romberg(float e) {

}

三实验组织运行要求

实验前，由任课教师落实实验任务，每个学生事先编写好算法设计源程序代码。集中上机、调试并通过计算机图形可视化演示操作实例来测试、验证所学的数值分析理论。

四实验条件

为每个学生提供一台具有WINDOWS 98/XP/NT/2000操作系统的计算机；同时提供VC++/VB/JAVA/TC等集成语言开发环境来编程设计计算方法的上机实验。

五实验步骤

1、根据实验内容和算法流程图预先编好程序初稿，上机调试、运行。

2、用绘图函数绘出曲边梯形的面积并计算其值并输出正确的结果。

3、测试、分析它们的运行结果，比较两种方法的精度差异。

4、实验完毕后提交实验报告。

六实验报告

实验报告的内容包括：实验目的、实验内容、源程序清单和结果分析及心得体会。并分析它们的运行结果，并比较两种数值积分不同的特点和精度差异。

实验三——微分方程

实验学时：2

实验类型：综合/设计

实验要求：必修

一实验目的

通过本次实验，熟悉求解常微分方程初值问题的有关方法和理论，主要是欧拉法、改进欧拉法、四阶龙格库塔法，学会编制这三种方法的计算程序。体会这三种解法的功能、优缺点及适用场合。

二实验内容

本实验将欧拉法、改进欧拉法和经典四阶龙哥库塔法编写为一个类CDifferential_Equation，调用该类实现一阶微分方程初值问题的求解，并画出解函数的图像。CDifferential_Equation的框架如下class CDifferential_Equation

{

public:

CDifferential_Equation(float ma,float mb,float mh,float my0);

~CDifferential_Equation(void);

float f(float x, float y);

void Euler(void); //欧拉法

void Euler_Improved(void); //改进欧拉法

void Classical_Forth_Order_Runge_Kutta(void); //经典四阶龙哥库塔法

float a,b,h,y0; //分别为区间左端点、区间右端点、步长

float *x,*y; //存放各个结点和横、纵坐标

int n; //结点个数，构造函数中得到

};

其中有参数的构造函数为

CDifferential_Equation::CDifferential_Equation(float ma,float mb,float mh,float my0)

{

a=ma;b=mb;h=mh;

n=(b-a)/h;

x=new float[n+1];

y=new float[n+1];

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

x[i]=a+h*i;

y[0]=my0;

}

三实验组织运行要求

实验前，由任课教师落实实验任务，每个学生事先编写好算法设计源程序代码。集中上机、调试并通过计算机图形可视化演示操作实例来测试、验证所学的数值分析理论。

四实验条件

为每个学生提供一台具有WINDOWS 98/XP/NT/2000操作系统的计算机；同时提供VC++/VB/JAVA/TC等集成语言开发环境来编程设计计算方法的上机实验。

五实验步骤

1、根据实验内容和算法流程图预先编好程序初稿，上机调试、运行。

2、在屏幕上按适当比例和位置画出坐标轴及解的函数曲线输出正确的结果。

3、测试、分析它们的运行结果，比较三种方法的精度和误差。

4、实验完毕后提交实验报告。

六实验报告

实验报告的内容包括：实验目的、实验内容、源程序清单和结果分析及心得体会。并分析它们的运行结果，比较三种常微分方程求初值的精度和误差。

实验四——非线性方程

实验学时：2

实验类型：设计

实验要求：必修

一实验目的

通过本次实验，熟练的掌握方程求根的最基本、常用的运算方法和理论。主要有二分法、牛顿法、弦截法，并体会它们各自不同的特点及收敛速率。

二 实验内容

2.1 算法原理与底层计算程序

非线性方程f(x)=0。

（1）二分法

对于单根区间[a,b]，f(a)f(b)<0：

① 区间中点m 为(a+b)/2，如果|b-a|

② 如果f(a)=0（或f(b)=0），则a （或b ）即为方程的解（注意：程序中，两浮点数不可判断相等，即使两浮点数a 和b 理论上相等，但“a==b;”一般不成立）；

③ 如果f(a)f(m)<0，则对区间[a,m]调用函数自身（递归），或将m 赋值给b （循环），跳转①，

如果f(m)f(b)<0，则对区间[m,b]调用函数自身（递归），或将m 赋值给a （循环），跳转①。 float Bisection(float a, float b,float e)

{//a 和b 分别为区间左右端点，e 为容许误差，返回值为方程的解 //建议该函数用递归实现

//函数在[a,b]区间须连续、端点函数值异号且函数在该区间有唯一解 //每步的区间中点存入顺序表中，以备后面画图之用 }

（2）牛顿法 1()

'()

k k k k f x x x f x +=-

，0x 已知

float Newton(float x, float e)

{//x 为初值x0，e 为容许误差，返回值为方程的解

//每次得到的近似解存入顺序表中，以备后面画图之用 }

（3）快速弦截法

1

11()

()()

k k k k k k k x x x x f x f x f x -+--=--，0x 和1x 已知

float Secant(float x0, float x1, float e)

{//x0，x1分别为前两个近似解，用于启动迭代过程，e 为容许误差，返回值为方程的解 //将每次迭代结果存入顺序表中，以备后面画图之用 }

也可以定义多个待求解方程为函数，在二分法、牛顿法和弦截法的函数中增加函数指针参数，用以调用不同的待求解方程。

2.2 可视化交互式对话框设计

（1）绘制直角坐标系，并画出待求解方程的函数图象。

（2）使用文本编辑框向程序传递a 、b 、x 、x0、x1、e 等条件初始值。

（3）使用三个按钮分别调用三种算法的函数，用静态文本框显示计算结果。对于区间二分法，用小圆圈在x 轴标出每次得到的区间中点；对于牛顿法，对每个k ，用直线连接(,

())k k x f x 和+1(,0)k x ；对于弦

截法，对每个k ，用直线连接(,())k k x f x 和2(,0)k x +。

图4 非线性方程

三实验组织运行要求

实验前，由任课教师落实实验任务，每个学生事先编写好算法设计源程序代码。集中上机、调试并通过计算机图形可视化演示操作实例来测试、验证所学的数值分析理论。

四实验条件

为每个学生提供一台具有WINDOWS 98/XP/NT/2000操作系统的计算机；同时提供VC++/VB/JAVA/TC等集成语言开发环境来编程设计计算方法的上机实验。

五实验步骤

1、根据实验内容和算法流程图预先编好程序初稿，上机调试、运行。

2、按适当的比例在屏幕上画出f(x)的函数曲线以及根在x轴上的逼近过程。

3、测试、分析它们的运行结果并比较三种方法的运算速率及精度误差。

4、实验完毕后提交实验报告。

六实验报告

实验报告的内容包括：实验目的、实验内容、源程序清单和结果分析及心得体会。并分析它们的运行结果，比较三种方程求根的运算速率及精度误差。

实验五——线性方程组

实验学时：2

实验类型：设计

实验要求：必修

一实验目的

通过本课程的实习，学会编写全主元消去法的计算程序。掌握解线性方程组的最基本算法及其运用，

进一步了解该解法的功能、优缺点，领会系数矩阵对解的影响。

二 实验内容

求解n 元线性方程组

1112131112122232223132333331

2

3n n n n n n nn n n a a a a x b a a a a x b a a a a x b a a a a x b ??????

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

??????=??????

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

即AX B =

其中第k 个方程表示为

),2,1(1

n i b x a

i

n

j j ij

==∑=

2.1 算法原理与底层计算程序

表9 线性方程组初始数据

（1）高斯——赛德尔迭代法

(1)()()()()

1

1221331441111

(1)(1)()()()

2

2112332442222

(1)(1)(1)()()3

3113223443333

(1)(1)4411441()1()1()1(k k k k k n n k k k k k n n k k k k k n n k k x a x a x a x a x b a x a x a x a x a x b a x a x a x a x a x b a x a x a a ++++++++=------+=------+=------+=-- (1)(1)()42243344(1)(1)(1)(1)(1)

1122331,11)

1()k k k n n k k k k k n n n n n n n n nn x a x a x b x a x a x a x a x b a +++++++---???

????????----+?

??

?=------+??

对于第i 个方程，表达式为

)(11

)(11)

1()

1(∑∑+=-=++--=n

i j k j ij i j k j ij i ii k i

x a x a b a x

该算法可以通过编写三个函数实现。

float *Adjust(float *pX0) 该函数计算一次迭代，pX0数组为上次迭代结果，返回值为本次迭代结果（数组） { //已有条件数据：系数矩阵数组A ，常数数组B ，未知数个数n

//计算结果数据：根据参数pX0数组（表示某步迭代结果）计算下一步迭代数组 float *pX1=(float *)malloc((n+1)*sizeof(float)); //不使用第0个元素

for(int i=1;i<=n;i++) //求每个(1)

k i

x +

{

//)(11

)(11)

1()1(∑∑+=-=++--=n

i j k j ij i j k j ij i ii k i

x a x a b a x

，两个一重循环 }

return pX1;

}

float Distance(float *pX0,float *pX1) //该函数计算两个向量的距离（即计算两个向量之差的范数） { //已有条件数据：向量维数n

}

float *Gauss_Seidel(float *pX0,float e) //该函数计算以pX0为初值，精度为e 的迭代解 { //已有条件数据：系数矩阵数组A ，常数数组B ，未知数个数n

//计算结果数据：根据参数pX0数组和e （表示某步迭代结果）计算迭代结果 //该函数调用Adjust(float *pX0)和Distance(float *pX0,float *pX1) }

（2）高斯列主元消去法

将线性方程组写成增广矩阵形式

11

12131121222322313233331

2

3n n n n n n nn n a a a a b a a a a b a a a a b a a a a b ??

????????

??????

将该增广矩阵作为全局数据或类的数据成员。

表10 高斯列主元消去法数据

所谓消去法，就是通过行初等变换，将增广矩阵化为对角矩阵的方法，最后一列就是方程组的解。因

此，首先编写对数组G 进行行初等变换的三个函数。 void Exchange(int i,int j) //交换G 的第i 行和第j 行 {

}

void Mult(int i,float r) //将G 的第i 行乘以r { }

void AddTo(int i,int j,float r) //将第i 行乘以r 加到第j 行 { } 另外还需选主元函数。

int GetMaxElem(int i) //查找G[i][i]……G[n][i]中最大元素的行标

{ }

有了这些函数，可以很方便的实现高斯列主元消去法。

void GaussEliminate()

{

for(int y=1;y<=n;y++) //将G化为上三角矩阵

{

//消掉G[y][y]以下的元素（简单起见，可以先将第y行除以G[y][y]）

}

for(int y=n;y>=1;y--) //将上三角矩阵化为对角矩阵

{

//消掉G[y][y]以上的元素

}

//G的最后一列就是方程组的解

}

对于C++语言，为了方便，使用二维动态数组。可调用下面的函数为二维指针分配动态存储空间和销毁空间。为了方便程序的移植和重复利用，可将这两个函数写在一个文件MatrixAllocate.h中。C#和Java可以直接使用二维数组。

template

bool New(T **&ppData,int H,int W)

{//该函数为二维指针ppData分配动态存储空间

T *ppData1 = new T[W*H];

ppData=new T*[H];

for(int i=0;i

{

ppData[i]=ppData1+i*W;

}

return true;

}

template

bool Delete(T **&ppData)

{//该函数销毁二维指针ppData

T *ppData1=ppData[0];

delete ppData1;

delete ppData;

ppData=NULL;

return true;

}

2.2 可视化交互式对话框设计

该实验不需要对话框，建立控制台应用程序，将结果输出即可。

三实验组织运行要求

实验前，由任课教师落实实验任务，每个学生事先编写好算法设计源程序代码。集中上机、调试并通过计算机图形可视化演示操作实例来测试、验证所学的数值分析理论。

四实验条件

为每个学生提供一台具有WINDOWS 98/XP/NT/2000操作系统的计算机；同时提供VC++/VB/JAVA/TC等集成语言开发环境来编程设计计算方法的上机实验。