计算方法实验报告高尚版.doc

1、实验一插值方法实验学时：4实验类型：设计实验要求：必修一 实验目的通过本次上机实习，能够进一步加深对各种插值算法的理解；学会使用用三种类型的插值函数的数学模型、基本算法，结合相应软件（如VC/VB/Delphi/Matlab/JAVA/Turbo C）编程实现数值方法的求解。并用该软件的绘图功能来显示插值函数，使其计算结果更加直观和形象化。二 实验内容通过程序求出插值函数的表达式是比较麻烦的，常用的方法是描出插值曲线上尽量密集的有限个采样点，并用这有限个采样点的连线，即折线，近似插值曲线。取点越密集，所得折线就越逼近理论上的插值曲线。本实验中将所取的点的横坐标存放于动态数组中，通过插值方法计算

2、得到的对应纵坐标存放于动态数组中。以Visual C+.Net 2005为例。本实验将Lagrange插值、Newton插值和三次样条插值实现为一个C+类CInterpolation，并在Button单击事件中调用该类相应函数，得出插值结果并画出图像。CInterpolation类为class CInterpolationpublic: CInterpolation();/构造函数CInterpolation(float *x1, float *y1, int n1);/结点横坐标、纵坐标、下标上限 CInterpolation();/析构函数int n, N;/结点下标上限，采样点下标上限f

3、loat *x, *y, *X;/分别存放结点横坐标、结点纵坐标、采样点横坐标float *p_H,*p_Alpha,*p_Beta,*p_a,*p_b,*p_c,*p_d,*p_m;/样条插值用到的公有指针，分别存放，和;其中，有参数的构造函数为CInterpolation(float *x1, float *y1, int n1) /动态数组x1，y1中存放结点的横、纵坐标，n1是结点下标上限（即n1+1个结点） n=n1; N=x1n-x10; X=new floatN+1; x=new floatn+1; y=new floatn+1; for(int i=0;i=n;i+) xi=x

4、1i; yi=y1i; for(int i=0;i=N;i+) Xi=x0+i;2.1 Lagrange插值，其中对于一个自变量，要求插值函数值，首先需要计算对应的Lagrange插值基函数值float l(float xv,int i) /求插值基函数的值 float t=1; for(int j=0;j=n;j+) if(j!=i) t=t*(xv-xj)/(xi-xj); return t;调用函数l(float x,int i)，可求出float p_l(float x) /求在一个点的插值结果 float t=0; for(int i=0;i=n;i+) t+=yi*l(x,i);

5、return t;调用p_l(float x)可实现整个区间的插值float *Lagrange() /求整个插值区间上所有采样点的插值结果 float *Y=new floatN+1; for(int k=0;k=N;k+) Yk=p_l(x0+k*h); return Y;2.2 Newton插值，其中，对于一个自变量，要求插值函数值，首先需要计算出和float *f()/该函数的返回值是一个长度为n1的动态数组，存放各阶差商float w(float x, int i)/该函数计算在求的函数中调用*f()得到各阶差商，然后在循环中调用w(float x)可得出插值结果float p_n(

6、float x)/该函数计算在一点的值调用p_n(float x)可实现整个区间的插值float *Newton()/该函数计算出插值区间内所有点的值2.3 三次样条插值三次样条插值程序可分为以下四步编写：（1） 计算结点间的步长、；（2） 利用、产生三对角方程组的系数矩阵和常数向量；（3） 通过求解三对角方程组，得出中间结点的导数；（4） 对自变量x，在对应区间上，使用Hermite插值；（5） 调用上述函数，实现样条插值。将每步写成函数：(1)void GetH(void)/该函数计算数组void GetAlpha(void)/该函数计算数组void GetBeta(void)/该函数计算

7、数组(2)void Geta(void)/该函数计算数组下对角线 void Getb(void)/该函数计算数组主对角线 void Getc(void)/该函数计算数组上对角线void Getd(void)/该函数计算方程组右端常数项(3)float *Chasing(float *pa,float *pb,float *pc,float *pd,int n)/追赶过程，计算各点斜率(4)float F0(float x)/该函数计算函数的值float F1(float x)/该函数计算函数的值float P0(float x)/该函数计算函数的值float P1(float x)/函数计算函

8、数的值调用上述函数，实现三次样条插值float *Cubic_Spline()/该函数计算所有点的插值结果三 实验组织运行要求 实验前，由任课教师落实实验任务，每个学生事先编写好算法设计源程序代码。集中上机、调试并通过计算机图形可视化演示操作实例来测试、验证所学的数值分析理论。四 实验条件为每个学生提供一台具有WINDOWS 98/XP/NT/2000操作系统的计算机；同时提供VC+/VB/JAVA/TC等集成语言开发环境来编程设计计算方法的上机实验。五 实验步骤1 根据实验内容和算法流程图预先编好程序初稿，上机调试、运行， 输出正确的结果。2 三种类型的插值函数所生成的图形要显示在同一坐标轴

9、上，并用三种不同的颜色分别表示出来。3 分析它们的运行结果，并比较三种插值函数各自不同的特点。4 实验完毕后提交实验报告六 实验代码#pragma onceclass CInterpolationpublic: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,

10、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);flo

11、at 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.hCInterpolation:CInterpolation(void)CInterpolation:CInterpolation(void)CInterpolation:CInterpolation(float* x1, float* y1, int n1) n=n1; N=x1n-x10; X=n

12、ew floatN+1; x=new floatn+1; y=new floatn+1; for(int i=0;i=n;i+) xi=x1i; yi=y1i; for(int i=0;i=N;i+) Xi=x0+i;float CInterpolation:l(float xv, int i)float t=1; for(int j=0;j=n;j+) if(j!=i) t=t*(xv-xj)/(xi-xj); return t;float CInterpolation:p_l(float x) float t=0; for(int i=0;i=n;i+) t+=yi*l(x,i); ret

13、urn t;float* CInterpolation:lagrange(void)float *Y=new floatN+1; for(int k=0;k=N;k+) Yk=p_l(x0+k); return Y;float * CInterpolation:f()/该函数的返回值是一个长度为n1的动态数组，存放各阶差商dq=new floatn+1;for(int i=0;in;i+)dqi=0;for(int k=0;k=i;k+)float t=1;for(int j=0;j=i;j+)if(j!=k)t=t*(xk-xj);dqi=dqi+yk/t;return dq;float C

14、Interpolation:w(float x1,int i)/该函数计算w(x)if(i=0)return 1;elsefloat t=1;for(int j=0;ji;j+)t=t*(x1-xj);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+dqi*w(x1,i);return t;float* CInterpolation:Newton(vo

15、id)float *Y=new floatN+1;for(int i=0;i=N;i+)Yi = p_n(x0+i);return Y;void CInterpolation:GetH(void)/计算数组hip_H =new floatn;for(int i=0;in;i+)p_Hi=xi+1-xi;void CInterpolation:GetAlaph(void)/计算数组ap_Alpha =new floatn;for(int i=1;in;i+)p_Alphai=p_Hi-1/(p_Hi-1+p_Hi);void CInterpolation:GetBeta(void)/计算数组bp

16、Beta=new floatn;for(int i=1;in;i+)p_Betai=3*(1-p_Alphai)*(yi-yi-1)/p_Hi-1+p_Alphai*(yi+1-yi)/p_Hi);void CInterpolation:Geta(void)p_a=new floatn+1;for(int i=1;in+1;i+)p_ai=1-p_Alphai;void CInterpolation:Getb(void)p_b=new floatn+1;for(int i=1;in+1;i+)p_bi=2;void CInterpolation:Getc(void)p_c=new float

17、n+1;for(int i=1;in+1;i+)p_ci=p_Alphai;void CInterpolation:Getd(void)float m0=5.0,mn=20.0;p_d=new floatn+1;for(int i=2;in;i+)p_di=p_Betai;p_d1 =p_Beta1-(1-p_Alpha1)*m0;p_dn=p_Betan-p_Alphan*mn;void CInterpolation:Getm(void)/计算各点斜率p_m = new floatn + 1; p_m0 = 5.0; p_mn = 20.0;float *beta, *r, *s; beta

18、 = new floatn + 1; r = new floatn; s = new floatn + 1; beta1 = p_b1; r1 = p_c1 / beta1;for(int i = 2; i n; i+) betai = p_bi - ri - 1 * p_ai;ri = p_ci / betai; betan = p_bn - rn - 1 * p_an;s1 = p_d1 / p_b1; for(int i = 2; i 0; i-) p_mi = si - ri * p_mi + 1;float CInterpolation:F0(float x)return (x -

19、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 floatn; for (int i = 0; i n; i +) p_si = yi*F0(x1

20、 xi)/p_Hi) + yi + 1 * F1(x1 - xi) / p_Hi) + p_Hi*p_mi*P0(x1-xi)/p_Hi)+ p_Hi * p_mi + 1 * P1(x1 - xi) / p_Hi); int j; for(j = 0; j n; j+) if(xj = x1 & x1 = xj + 1) break; return p_sj; float* CInterpolation:Cubic_Spline(void)float *Y=new floatN+1; GetH();GetAlaph(); GetBeta(); Geta(); Getb(); Getc();

21、 Getd(); Getm(); for (int i = 0; i = N; i+) Yi = S3(X0+i); return Y; void C实验一Dlg:OnPaint()if (IsIconic()CPaintDC dc(this); / 用于绘制的设备上下文SendMessage(WM_ICONERASEBKGND, reinterpret_cast(dc.GetSafeHdc(), 0);int cxIcon = GetSystemMetrics(SM_CXICON);/ 使图标在工作区矩形中居中int cyIcon = GetSystemMetrics(SM_CYICON);

22、CRect rect;GetClientRect(&rect);int x = (rect.Width() - cxIcon + 1) / 2;int y = (rect.Height() - cyIcon + 1) / 2;dc.DrawIcon(x, y, m_hIcon);/ 绘制图标elseCDialogEx:OnPaint();pWnd-Invalidate(); /使static控件区域无效 pWnd-UpdateWindow();/更新窗口此时才真正向系统发送重绘消息没有这句绝对出问题CDC *pDC = pWnd-GetDC(); /获取控件的CDC指针pDC-Rectangl

23、e(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 MoveTo(i

24、 40 + 20, 360);/绘制X轴刻度pDC -LineTo(i * 40 + 20, 355);for(i = 1; i 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

25、pOldPen);/恢复以前的画笔delete pPenRed;if(pPen != NULL)delete pPen；ReleaseDC(pDC);HCURSOR C实验一Dlg:OnQueryDragIcon()return static_cast(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);Rel

26、easeDC(pDC);CDialog:OnLButtonDblClk(nFlags, point);CDialogEx:OnLButtonDblClk(nFlags, point);void C实验一Dlg:OnBnClickedButton1()int count = g_arrPoint.size(); float *x = new floatcount; float *y = new floatcount; int i = 0; std:vector:iterator iter;for (iter = g_arrPoint.begin(); iter != g_arrPoint.end

27、); iter+) xi = (*iter).x; yi = (*iter).y; i+; CInterpolation CIP = CInterpolation(x, y, count - 1); CIP.Y=CIP.lagrange();CDC *pDC = GetDC(); pDC-SelectObject(new CPen(PS_SOLID, 0, RGB(255,100, 0); pDC-MoveTo(CIP.X0, CIP.Y0); for(int i = 1; i LineTo(CIP.Xi, CIP.Yi); ReleaseDC(pDC); void C实验一Dlg:OnBn

28、ClickedButton2()Invalidate(TRUE);g_arrPoint.size();void C实验一Dlg:OnBnClickedButton3()int count = g_arrPoint.size(); float *x = new floatcount; float *y = new floatcount; int i = 0; std:vector:iterator iter1;for (iter1 = g_arrPoint.begin(); iter1 != g_arrPoint.end(); iter1+) xi = (*iter1).x; yi = (*it

29、er1).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.X0, ND.Y0); for(int i = 1; i LineTo(ND.Xi, ND.Yi); ReleaseDC(pDC); void C实验一Dlg:OnBnClickedButton4()int count = g_arrPoint.size(

30、); float *x = new floatcount; float *y = new floatcount; int i = 0; std:vector:iterator iter2;for (iter2 = g_arrPoint.begin(); iter2 != g_arrPoint.end(); iter2+) xi = (*iter2).x; yi = (*iter2).y; i+; CInterpolation CS = CInterpolation(x, y, count - 1); CS.Y=CS.Cubic_Spline();CDC *pDC = GetDC(); pDC-

31、SelectObject(new CPen(PS_SOLID, 0, RGB(0,255,0); pDC-MoveTo(CS.X0, CS.Y0); for(int i = 1; i LineTo(CS.Xi, CS.Yi); ReleaseDC(pDC);七 实验结果八实验二数值积分实验学时：2实验类型：设计实验要求：必修一 实验目的通过该课程实习，学会使用数值积分的各种方法求解定积分计算的问题，体会各种方法的精度差异。二 实验内容本实验将梯形法的递推化和龙贝格算法编写为类CIntegration，在工程中调用该类得出积分结果并画出区间二分图像。首先画出被积函数图像，并储存被积函数在图像中每

32、个像素对应的采样点的横坐标（存放于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 floatn+1; /存放区间二分点x0=a;x1=b;另需被积函数float CInteg

33、ration:f(float x)float mTemp=sin(x); /可替换为其他函数return mTemp;绘制函数图像时，需要将图像按照比例放大。2.1 梯形法的递推化 梯形法递推化公式为该过程是下列过程的往复循环 求已有各区间中点函数值的和； 求积分值。 区间二分，即更新结点，并将步长减半。要实现该算法过程，需要先编写下面两个函数float CIntegration:Sum(void)/该函数求的值void CIntegration:Divide(void)/该函数实现区间二分和步长更新float *px=new float2*n+1;for(int i=0;in;i+)px2*

34、i=xi;px2*i+1=(xi+xi+1)/2;px2*n=xn;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 龙贝格算法其中调用梯形法递推函数Trapezium_R

35、ecurrence(int N)可以很方便的实现单次龙贝格算法的函数float CIntegration:R1(float * T)调用函数R1可实现迭代的龙贝格函数float CIntegration:Romberg(float e)三 实验组织运行要求实验前，由任课教师落实实验任务，每个学生事先编写好算法设计源程序代码。集中上机、调试并通过计算机图形可视化演示操作实例来测试、验证所学的数值分析理论。四 实验条件为每个学生提供一台具有WINDOWS 98/XP/NT/2000操作系统的计算机；同时提供VC+/VB/JAVA/TC等集成语言开发环境来编程设计计算方法的上机实验。五 实验步骤1、

36、 根据实验内容和算法流程图预先编好程序初稿，上机调试、运行。2、 用绘图函数绘出曲边梯形的面积并计算其值并输出正确的结果。3、 测试、分析它们的运行结果，比较两种方法的精度差异。4、 实验完毕后提交实验报告。六 实验代码#pragma onceclass CIntegrationpublic:CIntegration(void);CIntegration(int ax,int ay,int len);/原点坐标, CIntegration(void);CIntegration(double ma, double mb);double *p_X;/像素对应的采样点的横坐标int* p_Y;/像素

37、对应的采样点的纵坐标double a,b,h,*x,*sT,*tempx;int n,tsize,count;double f(double xx);double Sum(void);void Divide(void);double T1ToT2(double T1);double Trapezium_Recurrence(double e);double Trapezium_Recurrence(int N);double R1(double* T,int N);double Romberg(double e);void init(void);#include stdafx.h#includ

38、e Integration.h#includemath.h#includeiostreamusing namespace std;CIntegration:CIntegration(void)CIntegration:CIntegration(void)CIntegration:CIntegration(int ax,int ay,int len)/原点坐标, p_X = new doublelen; p_Y = new intlen; for(int i=0;ilen;i+) p_Xi = i/80.0; CIntegration:CIntegration(double ma, double mb)a=ma; /a、b分别为积分区间左右端点b=mb;n=1; /区间二分个数h=b-a; /步长x=new doublen+1; /存放区间二分点x0=a;x1=b;