北航研究生数值分析编程大作业3.docx

1、北航研究生数值分析编程大作业3数值分析计算实习三姓名： 学号： 学院：航空科学与工程学院 具体源程序见光盘一、开发平台的选择本题中选择VC+ 2008作为开发平台，利用VC+作为开发环境，编了一个基于MFC的程序。可以使用Microsoft Visual Studio 2008打开“源文件”文件夹中的“Third.sln”查看或运行源文件。本程序不仅仅是界面不同，编程的思想和运算的速度都有提高。其优点是通过设计类，将Gauss消元法、Newton迭代法、分片二次差值、最小二乘法等功能封装成类的成员函数，在定义了对象后可以任意调用。这都是利用了封装性。在定义类的对象之后可以方便地调用这些成员变量

2、和成员函数。应用面向对象的编程思想还有一个最大的优点是可扩展性。想添加或改变功能只需要改动成员函数就行了而不会影响到程序的其他部分。同时，利用MFC还可以开发出很友好的程序界面。本次程序的初始化界面如下：图1 程序初始化界面初始化的界面中会显示题中所给的数表。二、算法的设计方案1、求出， ()共231组，并将其带入题中方程组，使用Newton迭代法，分别求出与之对应的值。2、对于已求出的使用分片二次代数插值，对原题中关于的数表进行插值，从而得到相应的231个。将第一步中求出的与之对应，便得到了231个zi=f(xi,yi)的值。3、令k初始值为0并逐渐增大，并使用最小二乘法曲面拟合对zi=f(

3、xi,yi)进行拟合，得到每次的值。当时结束计算，输出拟合结果及相应值。4、针对，采用步骤2中的方法计算出，再用步骤3中的方法计算出的值。将这两组值列出并计算绝对误差，以观察逼近的程度。三、源代码在光盘中有全部的源代码，请用Microsoft Visual Studio 2008打开。同时还有一个可执行文件可以直接执行。由于有些程序段是MFC直接生成的，这里没有把那部分程序打出。主要源程序代码如下：1、Newton类的设计：（因为这里应用了Newton迭代法，故将类的名称暂时取名为Newton类）类的头文件为：#pragma once/使此文件只编译一次#define ACCURACY 1.0

4、e-12/定义精度#define sigma 1.0e-7/定义计算sigma时的精度#define n 4/方程组中方程的个数class CNewtonpublic:/定义成员变量 double x11,y21;/需要输出的xi,yi double m_Sigma100;/存储的的数组 double x18,y15;/存储x*,y*的数组public: CNewton(void); CNewton(void); /求解向量行范数 double Get_Norm(double an); /利用选主元高斯消去法求解Newton法中的线性方程组 void Gauss(double fn,doubl

5、e daofnn,double xn); /牛顿迭代法 void Newton(double x11,double y21,double u11*21,double t11*21); /分片二次差值 double Fenpian(double uu,double tt); /利用LU分解法求解逆矩阵列向量 void Doolittle(double A1010,double b10,double x10,int N); /矩阵求逆 void Matrix_Inversion(double a1010,int K); /最小二乘拟合 int Least_Squares_Collocation(d

6、ouble x11,double y21,double z1121,double C1021);类的实现文件为：#include StdAfx.h/程序总的头文件#include Newton.h/将Newton类的头文件包含进来#include CNewton:CNewton(void)/类的默认构造函数 int i,j; for(i=0;i11;i+) xi=0.08*i; for(i=0;i21;i+) yi=0.5+0.05*i; for (i=0;i100;i+) m_Sigmai=0; for(i=0;i8;i+) x1i=0.1*(i+1); for(j=0;j5;j+) y1j

7、=0.5+0.2*(j+1);CNewton:CNewton(void)/类的默认析构函数double CNewton:Get_Norm(double an)/求解向量行范数 int i; double max; max=fabs(a0); for(i=0;imax) max=fabs(ai); return(max);void CNewton:Gauss(double fn,double daofnn,double xn) /利用选主元高斯消去法求解Newton法中的线性方程组 int i,j,k; double ann+1; double m,q; double (*pn)n+1,(*P)

8、n+1; for(i=0;in;i+) for(j=0;jn;j+) aij=daofij; for(i=0;in;i+) ain=-fi; for(i=0;in;i+) pi=a+i; for(i=0;in-1;i+) for(int s=i+1;s*(*pi+i) P=pi;pi=ps;ps=P; for(k=i+1;kn;k+) m=*(*pk+i)/(*(*pi+i); for(j=i;j=0;i-) q=0; for(j=i+1;jn;j+) q+=xj*(*(*pi+j); xi=(*(*pi+n)-q)/(*(*pi+i); void CNewton:Newton(double

9、x11,double y21,double u11*21,double t11*21)/牛顿迭代法 int i,j,k=0; double fn,v231,w231,deltan,mn; double daofnn; double E; for(i=0;i231;i+) ui=1; ti=2; vi=1; wi=2; for(i=0;i11;i+) for(j=0;j=ACCURACY); k+; double CNewton:Fenpian(double uu,double tt) /分片二次差值 double uc=0.4,tc=0.2,zz=0,lu3,lt3; double z66=

10、-0.5, -0.34, 0.14, 0.94, 2.06, 3.5, -0.42, -0.5, -0.26, 0.3, 1.18, 2.38, -0.18, -0.5, -0.5, -0.18, 0.46, 1.42, 0.22, -0.34, -0.58, -0.5, -0.1, 0.62, 0.78, -0.02, -0.5, -0.66, -0.5, -0.02, 1.5, 0.46, -0.26, -0.66, -0.74, -0.5, ; double u6=0, 0.4, 0.8, 1.2, 1.6, 2.0, t6=0, 0.2, 0.4, 0.6, 0.8, 1.0; int

11、 i=0,j=0; int up=0,tp=0; for(i=1;i5;i+) if(uc*i-0.2=uu)=2) up=1; break; else continue; for(j=1;j5;j+) if(tc*j-0.1=tt)=2) tp=1;break; else continue; if(up=0) if(uu1.8) i=4; else i=1; if(tp=0) if(tt0.9) j=4; else j=1; lu0=(uu-ui)*(uu-ui+1)/(ui-1-ui)/(ui-1-ui+1); lu1=(uu-ui-1)*(uu-ui+1)/(ui-ui-1)/(ui-u

12、i+1); lu2=(uu-ui-1)*(uu-ui)/(ui+1-ui-1)/(ui+1-ui); lt0=(tt-tj)*(tt-tj+1)/(tj-1-tj)/(tj-1-tj+1); lt1=(tt-tj-1)*(tt-tj+1)/(tj-tj-1)/(tj-tj+1); lt2=(tt-tj-1)*(tt-tj)/(tj+1-tj-1)/(tj+1-tj); for(int a=0;a3;a+) for(int b=0;b3;b+) zz+=lta*lub*za+j-1b+i-1; return(zz);void CNewton:Doolittle(double A1010,doub

13、le b10,double x10,int N)/利用LU分解法求解逆矩阵列向量 int i,j; double a1011; for(i=0;iN;i+) for(j=0;jN;j+) aij=Aij; for(i=0;iN;i+) aiN=bi; double p=0,q=0,z=0,w=0; for(j=0;jN+1;j+) a0j=a0j; for(i=1;iN;i+) ai0=ai0/a00; for(int m=1;mN;m+) for(j=m;jN+1;j+) for(int k=0;km;k+) p=p+amk*akj; amj=amj-p; p=0; for(i=m+1;iN

14、;i+) for(int k=0;k=0;i-) for(int t=N-1;ti;t-) w+=ait*xt; xi=(aiN-w)/aii; w=0; void CNewton:Matrix_Inversion(double a1010,int K) /矩阵求逆 int i,j,k; double b10,x10; double a_11010; for(k=0;kK;k+) for(i=0;iK;i+) bi=0; bk=1; Doolittle(a,b,x,K); for(i=0;iK;i+) a_1ik=xi; for(i=0;iK;i+) for(j=0;jK;j+) aij=a_

15、1ij;int CNewton:Least_Squares_Collocation(double x11,double y21,double z1121,double C1021)/最小二乘拟合 int K=0; int i,j,k,s,r; double BTB1010,BTB_BT1011,A1021,GTG1010; double he,delta,min; do K+; for(i=0;iK;i+) for(j=0;jK;j+) he=0; for(k=0;k11;k+) he+=pow(xk,i)*pow(xk,j); BTBij=he; Matrix_Inversion(BTB,K

16、); for(i=0;iK;i+) for(j=0;j11;j+) he=0; for(k=0;kK;k+) he+=BTBik*pow(xj,k); BTB_BTij=he; for(i=0;iK;i+) for(j=0;j21;j+) he=0; for(k=0;k11;k+) he+=BTB_BTik*zkj; Aij=he; for(i=0;iK;i+) for(j=0;jK;j+) he=0; for(k=0;k21;k+) he+=pow(yk,i)*pow(yk,j); GTGij=he; Matrix_Inversion(GTG,K); for(i=0;iK;i+) for(j

17、=0;jK;j+) he=0; for(k=0;k21;k+) he+=Aik*pow(yk,j); BTBij=he; for(i=0;iK;i+) for(j=0;jK;j+) he=0; for(k=0;kK;k+) he+=BTBik*GTGkj; Cij=he; min=0; for(i=0;i11;i+) for(j=0;j21;j+) delta=0; for(r=0;rK;r+) for(s=0;ssigma); return(K);2、对话框类及算法的驱动函数因为VC+程序是消息驱动的，真正的主函数是Windows调用的，它的功能是不断地从Windows系统的消息堆栈中读取消

18、息并用消息驱动相应的消息响应函数。所以在算法这一部分没有主函数，而消息响应函数起到了主函数的作用。每当我们激发一个消息（比如按下一个按钮）激发的消息将使消息响应函数开始运行，从而进行计算。窗体类的头文件：/ ThirdDlg.h : 头文件#pragma once/使此文件只编译一次#include afxcmn.h#include Newton.h/ CThirdDlg 对话框class CThirdDlg : public CDialog/ 构造public: CThirdDlg(CWnd* pParent = NULL); / 标准构造函数/ 对话框数据 enum IDD = IDD_T

19、HIRD_DIALOG ; protected: virtual void DoDataExchange(CDataExchange* pDX); / DDX/DDV 支持/ 实现protected: HICON m_hIcon; / 生成的消息映射函数 virtual BOOL OnInitDialog(); afx_msg void OnSysCommand(UINT nID, LPARAM lParam); afx_msg void OnPaint(); afx_msg HCURSOR OnQueryDragIcon(); DECLARE_MESSAGE_MAP()public: CLi

20、stCtrl m_Newton; CString m_Change; afx_msg void OnBnClickedButtonCul();/点击“计算”按钮后的驱动函数 CNewton m_a;/定义一个Newton类的对象 CListCtrl m_List_K_S; CListCtrl m_List_C; CListCtrl m_List_Compare;窗体类的实现文件：/ ThirdDlg.cpp : 实现文件#include stdafx.h#include Third.h#include ThirdDlg.h#include #ifdef _DEBUG#define new DE

21、BUG_NEW#endif/这里省略了默认的初始化程序 / TODO: 在此添加额外的初始化代码 /*初始化控件的程序段*/ m_Newton.SetView(LVS_REPORT);/初始化控件显示插值得到的数表 m_Newton.InsertColumn(0,m_Change,LVCFMT_LEFT,40,20);/显示标号栏最左上角的单元 for (int i=1;i=2;i+)/初始化列号 m_Change.Format(_T( %d),i);/将数字转化为字符（下同） m_Newton.InsertColumn(i,m_Change,LVCFMT_LEFT,70,20);/插入列并设置列号 m_Change.Format(_T( %d),3); m_Newton.InsertColumn(3,m_Change,LVCFMT_LEFT,206,20);/插入列并设置列号 for (int i=0;i231;i+)/初始化行号 m_Change.Format(_T(%d),231-i); m_Newton.InsertItem(0,m_Change);/插入行并设置行号 int t=0; for (int i=0;i11;i+) for (int j=0;j21;j+) m_Change.For