实验一——插值方法实验学时: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 插值()()nn i i i P x y l x ==∑,其中0,()nj i j j ni jx 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 xfloat 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.2Newton 插值010()(,,)()nn i i i P x f x x x x ω==∑,其中101,0()(),0i i j j i x x x i ω-==⎧⎪=⎨-≠⎪⎩∏,0100,()(,,)()ik i nk k j j j kf 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 ) {//该函数计算数组i hi}void GetAlpha(void){//该函数计算数组iα}void GetBeta(void){//该函数计算数组iβ}(2)void Geta(void){//该函数计算数组下对角线ia}void Getb(void){//该函数计算数组主对角线ib}void Getc(void){//该函数计算数组上对角线ic}void Getd(void){//该函数计算方程组右端常数项id}(3)float *Chasing(float *pa,float *pb,float *pc,float *pd,int n) {//追赶过程,计算各点斜率im}(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.实验完毕后提交实验报告六实验程序NInterpolation.cpp#include "stdafx.h"#include "NInterpolation.h"NInterpolation::NInterpolation(void){}NInterpolation::NInterpolation(double *x,double *y,int n){this->n=n;this->x=new double[n+1];this->y=new double[n+1];for(int i=0;i<=n;i++){this->x[i]=x[i];this->y[i]=y[i];}N=x[n]-x[0];X=new double[N+1];Y=new double[N+1];for(int i=0;i<=N;i++)X[i]=x[0]+i;a=new double[n+1];for(int k=0;k<=n;k++){double t=0;for(int i=0;i<=k;i++){double m=1;for(int j=0;j<=k;j++){if(j==i)continue;elsem*=(x[i]-x[j]);}t+=y[i]/m;}a[k]=t;}}NInterpolation::~NInterpolation(void) {}void NInterpolation::Netwon(void){for(int i=0;i<=N;i++){Y[i]=0;for(int j=0;j<=n;j++){Y[i]+=((l(X[i],j)*a[j]));}}}double NInterpolation::l(double mx,int p){double m=1;for(int i=0;i<p;i++){m*=(mx-x[i]);}return m;}ZInterpolation.cpp#include "StdAfx.h"#include "ZInterpolation.h"ZInterpolation::ZInterpolation(void){}ZInterpolation::ZInterpolation(double *x,double *y,int n) {this->n=n;this->x=new double[n+1];this->y=new double[n+1];for(int i=0;i<=n;i++){this->x[i]=x[i];this->y[i]=y[i];}N=x[n]-x[0];X=new double[N+1];Y=new double[N+1];for(int i=0;i<=N;i++)X[i]=x[0]+i;}ZInterpolation::~ZInterpolation(void){}double ZInterpolation::l(double mx, int i){double t=1;for(int j=0;j<=n;j++)if(j==i)continue;elset*=(mx-x[j])/(x[i]-x[j]);return t;}double ZInterpolation::L(double mx){double t=0;for(int i=0;i<=n;i++)t+=y[i]*l(mx,i);return t;}void ZInterpolation::lagrange(void){for(int i=0;i<=N;i++)Y[i]=L(X[i]);}x=[1 2 3 4];y=[12 -8 34 100];xi=1:0.1:4;yi=interp1(x,y,xi,'spline');plot(x,y,'o',xi,yi);七实验结果NInterpolation ZInterpolationCubic_Spline上述是第一个实验——插值方法程序的构思和编写过程。