北航研究生数值分析编程大作业3.docx
《北航研究生数值分析编程大作业3.docx》由会员分享，可在线阅读，更多相关《北航研究生数值分析编程大作业3.docx（36页珍藏版）》请在冰豆网上搜索。
北航研究生数值分析编程大作业3
数值分析
计算实习三
姓名：
学号：
学院：
航空科学与工程学院
具体源程序见光盘
一、开发平台的选择
本题中选择VC++2008作为开发平台，利用VC++作为开发环境，编了一个基于MFC的程序。
可以使用MicrosoftVisualStudio2008打开“源文件”文件夹中的“Third.sln”查看或运行源文件。
本程序不仅仅是界面不同，编程的思想和运算的速度都有提高。
其优点是通过设计类，将Gauss消元法、Newton迭代法、分片二次差值、最小二乘法等功能封装成类的成员函数，在定义了对象后可以任意调用。
这都是利用了封装性。
在定义类的对象之后可以方便地调用这些成员变量和成员函数。
应用面向对象的编程思想还有一个最大的优点是可扩展性。
想添加或改变功能只需要改动成员函数就行了而不会影响到程序的其他部分。
同时，利用MFC还可以开发出很友好的程序界面。
本次程序的初始化界面如下：
图1程序初始化界面
初始化的界面中会显示题中所给的数表。
二、算法的设计方案
1、求出
，
（
）共231组
，并将其带入题中方程组，使用Newton迭代法，分别求出与之对应的
值。
2、对于已求出的
使用分片二次代数插值，对原题中关于
的数表进行插值，从而得到相应的231个
。
将第一步中求出的
与之对应，便得到了231个zi=f（xi,yi）的值。
3、令k初始值为0并逐渐增大，并使用最小二乘法曲面拟合对zi=f（xi,yi）进行拟合，得到每次的
值。
当
时结束计算，输出拟合结果及相应
值。
4、针对
，采用步骤2中的方法计算出
，再用步骤3中的方法计算出
的值。
将这两组值列出并计算绝对误差，以观察
逼近
的程度。
三、源代码
在光盘中有全部的源代码，请用MicrosoftVisualStudio2008打开。
同时还有一个可执行文件可以直接执行。
由于有些程序段是MFC直接生成的，这里没有把那部分程序打出。
主要源程序代码如下：
1、Newton类的设计：
（因为这里应用了Newton迭代法，故将类的名称暂时取名为Newton类）
类的头文件为：
#pragmaonce//使此文件只编译一次
#defineACCURACY1.0e-12//定义精度
#definesigma1.0e-7//定义计算sigma时的精度
#definen4//方程组中方程的个数
classCNewton
{
public:
//定义成员变量
doublex[11],y[21];//需要输出的xi,yi
doublem_Sigma[100];//存储σ的的数组
doublex1[8],y1[5];//存储x*,y*的数组
public:
CNewton（void）;
~CNewton（void）;
//求解向量行范数
doubleGet_Norm（doublea[n]）;
//利用选主元高斯消去法求解Newton法中的线性方程组
voidGauss（doublef[n],doubledaof[n][n],doublex[n]）;
//牛顿迭代法
voidNewton（doublex[11],doubley[21],doubleu[11*21],doublet[11*21]）;
//分片二次差值
doubleFenpian（doubleuu,doublett）;
//利用LU分解法求解逆矩阵列向量
voidDoolittle（doubleA[10][10],doubleb[10],doublex[10],intN）;
//矩阵求逆
voidMatrix_Inversion（doublea[10][10],intK）;
//最小二乘拟合
intLeast_Squares_Collocation（doublex[11],doubley[21],doublez[11][21],doubleC[10][21]）;
};
类的实现文件为：
#include"StdAfx.h"//程序总的头文件
#include"Newton.h"//将Newton类的头文件包含进来
#include
CNewton:
:
CNewton（void）//类的默认构造函数
{
inti,j;
for（i=0;i<11;i++）
x[i]=0.08*i;
for（i=0;i<21;i++）
y[i]=0.5+0.05*i;
for（i=0;i<100;i++）
{m_Sigma[i]=0;}
for（i=0;i<8;i++）
x1[i]=0.1*（i+1）;
for（j=0;j<5;j++）
y1[j]=0.5+0.2*（j+1）;
}
CNewton:
:
~CNewton（void）//类的默认析构函数
{}
doubleCNewton:
:
Get_Norm（doublea[n]）//求解向量行范数
{
inti;
doublemax;
max=fabs（a[0]）;
for（i=0;iif（fabs（a[i]）>max）
max=fabs（a[i]）;
return（max）;
}
voidCNewton:
:
Gauss（doublef[n],doubledaof[n][n],doublex[n]）
//利用选主元高斯消去法求解Newton法中的线性方程组
{
inti,j,k;
doublea[n][n+1];
doublem,q;
double（*p[n]）[n+1],（*P）[n+1];
for（i=0;ifor（j=0;ja[i][j]=daof[i][j];
for（i=0;ia[i][n]=-f[i];
for（i=0;ip[i]=a+i;
for（i=0;i{
for（ints=i+1;sif（*（*p[s]+i）>*（*p[i]+i））
{P=p[i];p[i]=p[s];p[s]=P;}
for（k=i+1;k{
m=*（*p[k]+i）/（*（*p[i]+i））;
for（j=i;j<=n;j++）
*（*p[k]+j）=*（*p[k]+j）-m*（*（*p[i]+j））;
if（*（*p[k]+k）==0）
break;
}
}
if（*（*p[n-1]+n-1）==0）
;
else
{
x[n-1]=*（*p[n-1]+n）/（*（*p[n-1]+n-1））;
for（i=n-2;i>=0;i--）
{
q=0;
for（j=i+1;jq+=x[j]*（*（*p[i]+j））;
x[i]=（*（*p[i]+n）-q）/（*（*p[i]+i））;
}
}
}
voidCNewton:
:
Newton（doublex[11],doubley[21],doubleu[11*21],doublet[11*21]）//牛顿迭代法
{
inti,j,k=0;
doublef[n],v[231],w[231],delta[n],m[n];
doubledaof[n][n];
doubleE;
for（i=0;i<231;i++）
{
u[i]=1;
t[i]=2;
v[i]=1;
w[i]=2;
}
for（i=0;i<11;i++）
for（j=0;j<21;j++）
{
do
{
f[0]=0.5*cos（t[k]）+u[k]+v[k]+w[k]-x[i]-2.67;
f[1]=t[k]+0.5*sin
（u[k]）+v[k]+w[k]-y[j]-1.07;
f[2]=0.5*t[k]+u[k]+cos（v[k]）+w[k]-x[i]-3.74;
f[3]=t[k]+0.5*u[k]+v[k]+sin（w[k]）-y[j]-0.79;
daof[0][0]=1;
daof[0][1]=1;
daof[0][2]=1;
daof[0][3]=-0.5*sin（t[k]）;
daof[1][0]=0.5*cos（u[k]）;
daof[1][1]=1;
daof[1][2]=1;
daof[1][3]=1;
daof[2][0]=1;
daof[2][1]=-sin（v[k]）;
daof[2][2]=1;
daof[2][3]=0.5;
daof[3][0]=0.5;
daof[3][1]=1;
daof[3][2]=cos（w[k]）;
daof[3][3]=1;
m[0]=u[k];
m[1]=v[k];
m[2]=w[k];
m[3]=t[k];
Gauss（f,daof,delta）;
E=Get_Norm（delta）/Get_Norm（m）;
u[k]+=delta[0];
v[k]+=delta[1];
w[k]+=delta[2];
t[k]+=delta[3];
}while（E>=ACCURACY）;
k++;
}
}
doubleCNewton:
:
Fenpian（doubleuu,doublett）//分片二次差值
{
doubleuc=0.4,tc=0.2,zz=0,lu[3],lt[3];
doublez[6][6]={
{-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},
};
doubleu[6]={0,0.4,0.8,1.2,1.6,2.0},
t[6]={0,0.2,0.4,0.6,0.8,1.0};
inti=0,j=0;
intup=0,tp=0;
for（i=1;i<5;i++）
{
if（（uc*i-0.2=uu）==2）
{up=1;
break;}
else
continue;
}
for（j=1;j<5;j++）
{
if（（tc*j-0.1=tt）==2）
{tp=1;break;}
else
continue;
}
if（up==0）
{if（uu>1.8）
i=4;
else
i=1;}
if（tp==0）
{if（tt>0.9）
j=4;
else
j=1;}
lu[0]=（uu-u[i]）*（uu-u[i+1]）/（u[i-1]-u[i]）/（u[i-1]-u[i+1]）;
lu[1]=（uu-u[i-1]）*（uu-u[i+1]）/（u[i]-u[i-1]）/（u[i]-u[i+1]）;
lu[2]=（uu-u[i-1]）*（uu-u[i]）/（u[i+1]-u[i-1]）/（u[i+1]-u[i]）;
lt[0]=（tt-t[j]）*（tt-t[j+1]）/（t[j-1]-t[j]）/（t[j-1]-t[j+1]）;
lt[1]=（tt-t[j-1]）*（tt-t[j+1]）/（t[j]-t[j-1]）/（t[j]-t[j+1]）;
lt[2]=（tt-t[j-1]）*（tt-t[j]）/（t[j+1]-t[j-1]）/（t[j+1]-t[j]）;
for（inta=0;a<3;a++）
for（intb=0;b<3;b++）
zz+=lt[a]*lu[b]*z[a+j-1][b+i-1];
return（zz）;
}
voidCNewton:
:
Doolittle（doubleA[10][10],doubleb[10],doublex[10],intN）//利用LU分解法求解逆矩阵列向量
{
inti,j;
doublea[10][11];
for（i=0;ifor（j=0;ja[i][j]=A[i][j];
for（i=0;ia[i][N]=b[i];
doublep=0,q=0,z=0,w=0;
for（j=0;ja[0][j]=a[0][j];
for（i=1;ia[i][0]=a[i][0]/a[0][0];
for（intm=1;m{for（j=m;j{
for（intk=0;kp=p+a[m][k]*a[k][j];
a[m][j]=a[m][j]-p;
p=0;
}
for（i=m+1;i{for（intk=0;kq=q+a[i][k]*a[k][m];
a[i][m]=（a[i][m]-q）/a[m][m];
q=0;
}
}
x[N-1]=a[N-1][N]/a[N-1][N-1];
for（i=N-2;i>=0;i--）
{for（intt=N-1;t>i;t--）
w+=a[i][t]*x[t];
x[i]=（a[i][N]-w）/a[i][i];
w=0;
}
}
voidCNewton:
:
Matrix_Inversion（doublea[10][10],intK）//矩阵求逆
{
inti,j,k;
doubleb[10],x[10];
doublea_1[10][10];
for（k=0;k{
for（i=0;ib[i]=0;
b[k]=1;
Doolittle（a,b,x,K）;
for（i=0;ia_1[i][k]=x[i];
}
for（i=0;ifor（j=0;ja[i][j]=a_1[i][j];
}
intCNewton:
:
Least_Squares_Collocation（doublex[11],doubley[21],doublez[11][21],doubleC[10][21]）//最小二乘拟合
{
intK=0;
inti,j,k,s,r;
doubleBTB[10][10],BTB_BT[10][11],A[10][21],GTG[10][10];
doublehe,delta,min;
do
{
K++;
for（i=0;ifor（j=0;j{
he=0;
for（k=0;k<11;k++）
he+=pow（x[k],i）*pow（x[k],j）;
BTB[i][j]=he;
}
Matrix_Inversion（BTB,K）;
for（i=0;ifor（j=0;j<11;j++）
{
he=0;
for（k=0;khe+=BTB[i][k]*pow（x[j],k）;
BTB_BT[i][j]=he;
}
for（i=0;ifor（j=0;j<21;j++）
{
he=0;
for（k=0;k<11;k++）
he+=BTB_BT[i][k]*z[k][j];
A[i][j]=he;
}
for（i=0;ifor（j=0;j{
he=0;
for（k=0;k<21;k++）
he+=pow（y[k],i）*pow（y[k],j）;
GTG[i][j]=he;
}
Matrix_Inversion（GTG,K）;
for（i=0;ifor（j=0;j{
he=0;
for（k=0;k<21;k++）
he+=A[i][k]*pow（y[k],j）;
BTB[i][j]=he;
}
for（i=0;ifor（j=0;j{
he=0;
for（k=0;khe+=BTB[i][k]*GTG[k][j];
C[i][j]=he;
}
min=0;
for（i=0;i<11;i++）
for（j=0;j<21;j++）
{
delta=0;
for（r=0;rfor（s=0;sdelta+=C[r][s]*pow（x[i],r）*pow（y[j],s）;
delta=delta-z[i][j];
delta=delta*delta;
min+=delta;
}
m_Sigma[K-1]=min;
}while（min>sigma）;
return（K）;
}
2、对话框类及算法的驱动函数
因为VC++程序是消息驱动的，真正的主函数是Windows调用的，它的功能是不断地从Windows系统的消息堆栈中读取消息并用消息驱动相应的消息响应函数。
所以在算法这一部分没有主函数，而消息响应函数起到了主函数的作用。
每当我们激发一个消息（比如按下一个按钮）激发的消息将使消息响应函数开始运行，从而进行计算。
窗体类的头文件：
//ThirdDlg.h:
头文件
#pragmaonce//使此文件只编译一次
#include"afxcmn.h"
#include"Newton.h"
//CThirdDlg对话框
classCThirdDlg:
publicCDialog
{
//构造
public:
CThirdDlg（CWnd*pParent=NULL）;//标准构造函数
//对话框数据
enum{IDD=IDD_THIRD_DIALOG};
protected:
virtualvoidDoDataExchange（CDataExchange*pDX）;//DDX/DDV支持
//实现
protected:
HICONm_hIcon;
//生成的消息映射函数
virtualBOOLOnInitDialog（）;
afx_msgvoidOnSysCommand（UINTnID,LPARAMlParam）;
afx_msgvoidOnPaint（）;
afx_msgHCURSOROnQueryDragIcon（）;
DECLARE_MESSAGE_MAP（）
public:
CListCtrlm_Newton;
CStringm_Change;
afx_msgvoidOnBnClickedButtonCul（）;//点击“计算”按钮后的驱动函数
CNewtonm_a;//定义一个Newton类的对象
CListCtrlm_List_K_S;
CListCtrlm_List_C;
CListCtrlm_List_Compare;
};
窗体类的实现文件：
//ThirdDlg.cpp:
实现文件
#include"stdafx.h"
#include"Third.h"
#include"ThirdDlg.h"
#include
#ifdef_DEBUG
#definenewDEBUG_NEW
#endif
//这里省略了默认的初始化程序
//TODO:
在此添加额外的初始化代码
/*********************初始化控件的程序段************************/
m_Newton.SetView（LVS_REPORT）;//初始化控件显示插值得到的数表
m_Newton.InsertColumn（0,m_Change,LVCFMT_LEFT,40,20）;//显示标号栏最左上角的单元
for（inti=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（inti=0;i<231;i++）//初始化行号
{
m_Change.Format（_T（"%d"）,231-i）;
m_Newton.InsertItem（0,m_Change）;//插入行并设置行号
}
intt=0;
for（inti=0;i<11;i++）
{
for（intj=0;j<21;j++）
{
m_Change.For
升级会员 