计算方法实验报告(高尚版)

实验一——插值方法实验学时：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 xfloat 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.2 Newton 插值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 ) {//该函数计算数组下对角线i a }void Getb(void ) {//该函数计算数组主对角线i b }void Getc(void ) {//该函数计算数组上对角线i c }void Getd(void ) {//该函数计算方程组右端常数项i d }(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．实验完毕后提交实验报告六实验代码#pragma onceclass CInterpolation{public:CInterpolation(void);~CInterpolation(void);CInterpolation(float* x1, float* y1, int n1);float l(float xv, int i);float p_l(float x);float* lagrange(void);int n, N;//结点下标上限，采样点下标上限float *x, *y, *X, *Y,*dq;//分别存放结点横坐标、结点纵坐标、采样点横坐标float *p_H,*p_Alpha,*p_Beta,*p_a,*p_b,*p_c,*p_d,*p_m,*p_s;float *f();float w(float x1,int i);float p_n(float x);float* Newton(void);void GetH(void);void GetAlaph(void);void GetBeta(void);void Geta(void);void Getb(void);void Getc(void);void Getd(void);void Getm(void);float F0(float x);float F1(float x);float P0(float x);float P1(float x);float S3(float x1);float* Cubic_Spline(void);};#include"stdafx.h"#include"Lagrange.h"CInterpolation::CInterpolation(void){}CInterpolation::~CInterpolation(void){}CInterpolation::CInterpolation(float* x1, float* y1, int n1){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;}float CInterpolation::l(float xv, int i){float t=1;for(int j=0;j<=n;j++)if(j!=i)t=t*(xv-x[j])/(x[i]-x[j]);return t;}float CInterpolation::p_l(float x){float t=0;for(int i=0;i<=n;i++)t+=y[i]*l(x,i);return t;}float* CInterpolation::lagrange(void){float *Y=new float[N+1];for(int k=0;k<=N;k++)Y[k]=p_l(x[0]+k);return Y;}float * CInterpolation::f()//该函数的返回值是一个长度为n＋1的动态数组，存放各阶差商{for(int i=0;i<n;i++){dq[i]=0;for(int k=0;k<=i;k++){float t=1;for(int j=0;j<=i;j++)if(j!=k)t=t*(x[k]-x[j]);dq[i]=dq[i]+y[k]/t;}}return dq;}float CInterpolation::w(float x1,int i)//该函数计算w(x){if(i==0)return 1;else{float t=1;for(int j=0;j<i;j++)t=t*(x1-x[j]);return t;}}//在求的函数中调用*f()得到各阶差商，然后在循环中调用w(float x)可得出插值结果float CInterpolation::p_n(float x1)//该函数计算在一点的值{float *dq =f();float t=0;for(int i=0;i<=n;i++)t=t+dq[i]*w(x1,i);return t;}float* CInterpolation::Newton(void){float *Y=new float[N+1];for(int i=0;i<=N;i++)Y[i] = p_n(x[0]+i);return Y;}void CInterpolation::GetH(void)//计算数组hi{p_H =new float[n];for(int i=0;i<n;i++)p_H[i]=x[i+1]-x[i];}void CInterpolation::GetAlaph(void)//计算数组a{p_Alpha =new float[n];for(int i=1;i<n;i++)p_Alpha[i]=p_H[i-1]/(p_H[i-1]+p_H[i]);}void CInterpolation::GetBeta(void)//计算数组b{for(int i=1;i<n;i++){p_Beta[i]=3*((1-p_Alpha[i])*(y[i]-y[i-1])/p_H[i-1]+p_Alpha[i]*(y[i+1]-y[i])/p_H[i]) ;}}void CInterpolation::Geta(void){p_a=new float[n+1];for(int i=1;i<n+1;i++)p_a[i]=1-p_Alpha[i];}void CInterpolation::Getb(void){p_b=new float[n+1];for(int i=1;i<n+1;i++)p_b[i]=2;}void CInterpolation::Getc(void){p_c=new float[n+1];for(int i=1;i<n+1;i++)p_c[i]=p_Alpha[i];}void CInterpolation::Getd(void){float m0=5.0,mn=20.0;p_d=new float[n+1];for(int i=2;i<n;i++)p_d[i]=p_Beta[i];p_d[1] =p_Beta[1]-(1-p_Alpha[1])*m0;p_d[n]=p_Beta[n]-p_Alpha[n]*mn;}void CInterpolation::Getm(void)//计算各点斜率{p_m = new float[n + 1];p_m[0] = 5.0; p_m[n] = 20.0;float *beta, *r, *s;beta = new float[n + 1];r = new float[n];s = new float[n + 1];beta[1] = p_b[1]; r[1] = p_c[1] / beta[1];for(int i = 2; i < n; i++){beta[i] = p_b[i] - r[i - 1] * p_a[i];r[i] = p_c[i] / beta[i];}beta[n] = p_b[n] - r[n - 1] * p_a[n];s[1] = p_d[1] / p_b[1];for(int i = 2; i < n+1; i++)s[i] = (p_d[i] - s[i - 1] * p_a[i]) / beta[i];p_m[n - 1] = s[n - 1];for(int i = n - 2; i > 0; i--)p_m[i] = s[i] - r[i] * p_m[i + 1];}float CInterpolation::F0(float x){return (x - 1) * (x - 1) * (2 * x + 1);}float CInterpolation::F1(float x){return x * x * (-2 * x + 3);}float CInterpolation::P0(float x){return x * (x - 1) * (x - 1);}float CInterpolation::P1(float x){return x * x * (x - 1);}float CInterpolation::S3(float x1){p_s = new float[n];for (int i = 0; i < n; i ++){p_s[i] = y[i]*F0((x1 - x[i])/p_H[i]) + y[i + 1] * F1((x1 - x[i]) / p_H[i]) + p_H[i]*p_m[i]*P0((x1-x[i])/p_H[i])+ p_H[i] * p_m[i + 1] * P1((x1- x[i]) / p_H[i]);}int j;for(j = 0; j < n; j++){if(x[j] <= x1 && x1 <= x[j + 1]) break;}return p_s[j];}float* CInterpolation::Cubic_Spline(void){float *Y=new float[N+1];GetH();GetAlaph();GetBeta();Geta();Getb();Getc();Getd();Getm();for (int i = 0; i <= N; i++){Y[i] = S3(X[0]+i);}return Y;}void C实验一Dlg::OnPaint(){if (IsIconic()){CPaintDC dc(this); // 用于绘制的设备上下文SendMessage(WM_ICONERASEBKGND, reinterpret_cast<WPARAM>(dc.GetSafeHdc()), 0);int cxIcon = GetSystemMetrics(SM_CXICON);// 使图标在工作区矩形中居中int cyIcon = GetSystemMetrics(SM_CYICON);CRect rect;GetClientRect(&rect);int x = (rect.Width() - cxIcon + 1) / 2;int y = (rect.Height() - cyIcon + 1) / 2;dc.DrawIcon(x, y, m_hIcon); // 绘制图标}else{CDialogEx::OnPaint();}pWnd->Invalidate(); //使static控件区域无效pWnd->UpdateWindow();//更新窗口此时才真正向系统发送重绘消息没有这句绝对出问题CDC *pDC = pWnd->GetDC(); //获取控件的CDC指针pDC->Rectangle(0, 0, 500, 380); //画一个矩形框CPen *pPenRed = new CPen(); //创建画笔对象pPenRed ->CreatePen(PS_SOLID, 1, RGB(0, 0, 255)); //红色画笔CPen *pPen = NULL;CGdiObject*pOldPen = pDC ->SelectObject(pPenRed);//选中当前红色画笔并保存以前的画笔pDC ->MoveTo(20, 20); //绘制坐标轴pDC ->LineTo(20, 360); //竖起轴pDC ->LineTo(480, 360); //水平轴for(i = 0; i <= 10; i ++)//写X轴刻度{str.Format(_T("%ld"), i);pDC ->MoveTo(i * 40 + 20, 360); //绘制X轴刻度pDC ->LineTo(i * 40 + 20, 355);}for(i = 1; i <= 10; i ++) //写Y轴刻度值{str.Format(_T("%ld"), i);pDC ->MoveTo(25, 360 - 30 * i);////绘制Y轴刻度pDC ->LineTo(20, 360 - 30 * i);}pDC ->MoveTo(470, 357);//绘制X箭头pDC ->LineTo(480, 360);pDC ->LineTo(470, 363);pDC ->MoveTo(17, 30);//绘制Y箭头pDC ->LineTo(20, 20); //绘制左边箭头pDC ->LineTo(23, 30); //绘制右边箭头pDC ->SelectObject(pOldPen);//恢复以前的画笔delete pPenRed;if(pPen != NULL)delete pPen；ReleaseDC(pDC);}HCURSOR C实验一Dlg::OnQueryDragIcon(){return static_cast<HCURSOR>(m_hIcon);}void C实验一Dlg::OnLButtonDblClk(UINT nFlags, CPoint point)//双击取点函数{CDC *pDC=GetDC();pDC->Ellipse(point.x-3,point.y-3,point.x+3,point.y+3);g_arrPoint.push_back(point);ReleaseDC(pDC);CDialog::OnLButtonDblClk(nFlags, point);CDialogEx::OnLButtonDblClk(nFlags, point);}void C实验一Dlg::OnBnClickedButton1(){int count = g_arrPoint.size();float *x = new float[count];float *y = new float[count];int i = 0;std::vector<CPoint>::iterator iter;for (iter = g_arrPoint.begin(); iter != g_arrPoint.end(); iter++) {x[i] = (*iter).x;y[i] = (*iter).y;i++;}CInterpolation CIP = CInterpolation(x, y, count - 1);CIP.Y=grange();CDC *pDC = GetDC();pDC->SelectObject(new CPen(PS_SOLID, 0, RGB(255,100, 0)));pDC->MoveTo(CIP.X[0], CIP.Y[0]);for(int i = 1; i <= CIP.N; i++)pDC->LineTo(CIP.X[i], CIP.Y[i]);ReleaseDC(pDC);}void C实验一Dlg::OnBnClickedButton2(){Invalidate(TRUE);g_arrPoint.size();}void C实验一Dlg::OnBnClickedButton3(){int count = g_arrPoint.size();float *x = new float[count];float *y = new float[count];int i = 0;std::vector<CPoint>::iterator iter1;for (iter1 = g_arrPoint.begin(); iter1 != g_arrPoint.end(); iter1++) {x[i] = (*iter1).x;y[i] = (*iter1).y;i++;}CInterpolation ND = CInterpolation(x, y, count - 1);ND.Y=ND.Newton();CDC *pDC = GetDC();pDC->SelectObject(new CPen(PS_SOLID, 0, RGB(0,100, 255)));pDC->MoveTo(ND.X[0], ND.Y[0]);for(int i = 1; i <= ND.N; i++)pDC->LineTo(ND.X[i], ND.Y[i]);ReleaseDC(pDC);}void C实验一Dlg::OnBnClickedButton4(){int count = g_arrPoint.size();float *x = new float[count];float *y = new float[count];int i = 0;std::vector<CPoint>::iterator iter2;for (iter2 = g_arrPoint.begin(); iter2 != g_arrPoint.end(); iter2++) {x[i] = (*iter2).x;y[i] = (*iter2).y;i++;}CInterpolation CS = CInterpolation(x, y, count - 1);CS.Y=CS.Cubic_Spline();CDC *pDC = GetDC();pDC->SelectObject(new CPen(PS_SOLID, 0, RGB(0,255,0)));pDC->MoveTo(CS.X[0], CS.Y[0]);for(int i = 1; i <= CS.N; i++)pDC->LineTo(CS.X[i], CS.Y[i]);ReleaseDC(pDC);}七．实验结果八．实验二——数值积分实验学时：2实验类型：设计实验要求：必修一 实验目的通过该课程实习，学会使用数值积分的各种方法求解定积分计算的问题，体会各种方法的精度差异。