北航数值分析第三次大作业.doc

《数值分析》第三次大作业

姓名:袁二凯
学号:SY0917145
联系方式:13488854452
642101136@

2011年12月8日

一、 算法的设计方案：
（一）、总体方案设计：
（1）解非线性方程组。将给定的当作已知量代入题目给定的非线性方程组，求得与相对应的数组t[i][j],u[i][j]。
（2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=。
（3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k值，并得到曲面拟合的系数矩阵C[r][s]。
（4）观察和的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点对应的，再与对应的比较即可，这里求解可以直接使用（3）中的C[r][s]和k。
（二）具体算法设计：
（1）解非线性方程组
牛顿法解方程组的解，可采用如下算法：
1）在附近选取，给定精度水平和最大迭代次数M。
2）对于执行
① 计算和。
② 求解关于的线性方程组

③ 若，则取，并停止计算；否则转④。
④ 计算。
⑤ 若，则继续，否则，输出M次迭代不成功的信息，并停止计算。
（2）分片双二次插值
给定已知数表以及需要插值的节点，进行分片二次插值的算法：
设已知数表中的点为： ，需要插值的节点为。
1) 根据选择插值节点：
若或，插值节点对应取或，
若或，插值节点对应取或。
若
则选择为插值节点。
2）计算

插值多项式的公式为：

注：本步进行插值运算的是，利用与的对应关系就可以得到与的对应关系。
（3）曲面拟合
根据插值得到的数表进行曲面拟合的过程：
1） 根据拟合节点和基底函数写出矩阵B和G：

2） 计算 。
在这里，为了简化计算和编程、避免矩阵求逆，记：
，
对上面两式进行变形，得到如下两个线性方程组：
，
通过解上述两个线性方程组，则有：
3） 对于每一个， 。
4） 拟合需要达到的精度条件为：
。
其中对应着插值得到的数表中的值。
5） 让k逐步增加，每一次重复执行以上几步，直到
成立。此时的k值就是要求解最小的k。

二、 源程序：
#include<stdio.h>
#include<iostream>
#include<stdlib.h>
#include<math.h>
#include<float.h>
#include<iomanip>

#define Epsilon1 1e-12 /*解线性方程组时近似解向量的精度*/
#define M 200 /*解线性方程组时的最大迭代次数*/
#define N 10 /*求解迭代次数时假设的k的最大值，用于定义包含k的存储空间*/
void Newton(); /*牛顿法求解非线性方程组子程序*/
void fpeccz(); /*分片二次代数插值子程序*/
void qmnh(); /*曲面拟合子程序*/
void duibi(); /*对比𝑓和p逼近效果的子程序*/
double x[11],y[21],t[11][21],u[11][21];/*定义全局变量*/
double z[11][21],C[10][10];
double kz;

void Newton(double x[11],double y[21])/*牛顿法求解非线性方程组子程序*/
{
double X[4],dx[4],F[4],dF[4][4],temp,m,fx,fX;
int i,j,k,l,p,ik,n;

for(i=0;i<=10;i++)
{
for(j=0;j<=20;j++)
{
X[0]=1; /*选取迭代初始向量，四个分别代表t，u，v，w*/
X[1]=1;
X[2]=1;
X[3]=1;
n=0;

loop1:{ F[0]=0.5*cos(X[0])+X[1]+X[2]+X[3]-x[i]-2.67;
F[1]=X[0]+0.5*sin(X[1])+X[2]+X[3]-y[j]-1.07;
F[2]=0.5*X[0]+X[1]+cos(X[2])+X[3]-x[i]-3.74;
F[3]=X[0]+0.5*X[1]+X[2]+sin(X[3])-y[j]-0.79;
/*求解F（x）*/

dF[0][0]=-0.5*sin(X[0]); /*求解F'（x）*/
dF[0][1]=1;
dF[0][2]=1;
dF[0][3]=1;
dF[1][0]=1;
dF[1][1]=0.5*cos(X[1]);
dF[1][2]=1;
dF[1][3]=1;
dF[2][0]=0.5;
dF[2][1]=1;
dF[2][2]=-sin(X[2]);
dF[2][3]=1;
dF[3][0]=1;
dF[3][1]=0.5;
dF[3][2]=1;
dF[3][3]=cos(X[3]);

/*高斯选主元消去法求解Δx*/
for(k=0;k<3;k++)
{
ik=k;
for(l=k;l<=3;l++)
{if(dF[ik][k]<dF[l][k])
ik=l;
} /*选主元*/

temp=0;
temp=F[ik];
F[ik]=F[k];
F[k]=temp;

for(l=k;l<=3;l++)
{ temp=0;
temp=dF[ik][l];
dF[ik][l]=dF[k][l];
dF[k][l]=temp;
}

for(l=k+1;l<=3;l++)
{
m=dF[l][k]/dF[k][k];
F[l]=F[l]-m*F[k];
for(p=k+1;p<=3;p++)
{dF[l][p]=dF[l][p]-m*dF[k][p];}
} /*消去过程*/
}

dx[3]=-F[3]/dF[3][3];
for(k=2;k>=0;k--)
{
temp=0;
for(l=k+1;l<=3;l++)
{temp=temp+dF[k][l]*dx[l]/dF[k][k];}
dx[k]=-F[k]/dF[k][k]-temp;
}

temp=0;
for(l=0;l<=3;l++) /*求解矩阵范数，用无穷范数*/
{
if(temp<fabs(dx[l]))
temp=fabs(dx[l]);
}
fx=temp;

temp=0;
for(l=0;l<=3;l++)
{
if(temp<fabs(X[l]))
temp=fabs(X[l]);
}
fX=temp;

if(fabs(fx/fX)<Epsilon1) /*判断是否成立*/
{ t[i][j]=X[0];
u[i][j]=X[1];
goto loop4;}
else
{ for(l=0;l<=3;l++)
{X[l]=X[l]+dx[l];}
n=n+1;
goto loop3;}
}
loop3:{if(n<M) /*判断是否超出规定迭代次数*/
goto loop1;
else
printf("迭代不成功\n");
goto loop4; }

loop4:{continue;}
}
}
}

void fpeccz(double t[11][21],double u[11][21])/*分片二次代数插值子程序*/
{
int s[11][21],r[11][21];
int i,j,i1,j1,m;
double z0[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}};
double t0[6]={0,0.2,0.4,0.6,0.8,1.0};
double u0[6]={0,0.4,0.8,1.2,1.6,2.0};
double temp1,temp2;

for(i=0;i<=10;i++) /*选取插值节点*/
{
for(j=0;j<=20;j++)
{
if(t[i][j]<=0.3) s[i][j]=1;
else if(t[i][j]>0.7) s[i][j]=4;
else
{
for(m=2;m<=3;m++)
{
if((t[i][j]>0.2*m-0.1)&&(t[i][j]<=0.2*m+0.1))
{
s[i][j]=m;
}
}
}
}
}
for(i=0;i<=10;i++)
{
for(j=0;j<=20;j++)
{
if(u[i][j]<=0.6) r[i][j]=1;
else if(u[i][j]>1.4) r[i][j]=4;
else
{
for(m=2;m<=3;m++)
{
if((u[i][j]>0.4*m-0.2)&&(u[i][j]<=0.4*m+0.2))
{
r[i][j]=m;
}
}
}
}
}

for(i=0;i<=10;i++) /*插值运算*/
{
for(j=0;j<=20;j++)
{
z[i][j]=0;
for(i1=s[i][j]-1;i1<=s[i][j]+1;i1++)
{
for(j1=r[i][j]-1;j1<=r[i][j]+1;j1++)
{
temp1=1.0;
for(m=s[i][j]-1;m<=s[i][j]+1;m++)
{ if(m!=i1) temp1*=(t[i][j]-t0[m])/(t0[i1]-t0[m]); }
temp2=1.0;
for(m=r[i][j]-1;m<=r[i][j]+1;m++)
{ if(m!=j1) temp2*=(u[i][j]-u0[m])/(u0[j1]-u0[m]); }
z[i][j]+=temp1*temp2*z0[i1][j1];
}
}
}
}
}

void qmnh() /*曲面拟合子程序*/
{ int i,j,k,m,l,i1,j1,ik;
double A[N][21],D[N][21],B[11][N],Bt[N][11],BtB[N][N],BtU[N][21],BtB1[N][N];
double G[21][N],Gt[N][21],GtG[N][N],GtG1[N][N];
double sigma,p[11][21],temp,q;

printf("选择过程的k和sigma值:\n");
k=0;
sigma=1; /*选取初值，使循环开始*/

/*求解系数矩阵Crs*/
while(sigma>1e-7)
{
for(i=0;i<=10;i++)
{
for(j=0;j<=k;j++)
{
B[i][j]=pow(x[i],j);
Bt[j][i]=B[i][j];
}
}

for(i=0;i<=k;i++)
{
for(j=0;j<=k;j++)
{
temp=0;
for(l=0;l<=10;l++)
{ temp+=Bt[i][l]*B[l][j]; }
BtB[i][j]=temp;
}
}

for(i=0;i<=k;i++)
{
for(j=0;j<=20;j++)
{
temp=0;
for(l=0;l<=10;l++)
{ temp+=Bt[i][l]*z[l][j]; }
BtU[i][j]=temp;
}
}

for(l=0;l<=20;l++)
{
for(i=0;i<=k;i++)
{
for(j=0;j<=k;j++)
{ BtB1[i][j]=BtB[i][j]; }
}
for(m=0;m<=k-1;m++)
{
ik=m;
for(i=m;i<=k-1;i++)
{
if(fabs(BtB1[i][m])<fabs(BtB1[i+1][m]))
ik=i+1;
else ;
}
if(ik!=m)
{
for(i=m;i<=k;i++)
{
temp=BtB1[m][i];
BtB1[m][i]=BtB1[ik][i];
BtB1[ik][i]=temp;
}
temp=BtU[m][l];
BtU[m][l]=BtU[ik][l];
BtU[ik][l]=temp;
}

for(i=m+1;i<=k;i++)
{
q=BtB1[i][m]/BtB1[m][m];
for(j=m;j<=k;j++)
{ BtB1[i][j]=BtB1[i][j]-q*BtB1[m][j]; }
BtU[i][l]=BtU[i][l]-q*BtU[m][l];
}
}

A[k][l]=BtU[k][l]/BtB1[k][k];
for(m=k-1;m>=0;m--)
{
temp=0;
for(i=m+1;i<=k;i++)
{ temp+=A[i][l]*BtB1[m][i]; }
A[m][l]=(BtU[m][l]-temp)/BtB1[m][m];
}
}

for(i=0;i<=20;i++)
{
for(j=0;j<=k;j++)
{
G[i][j]=pow(y[i],j);
Gt[j][i]=G[i][j];
}
}

for(i=0;i<=k;i++)
{
for(j=0;j<=k;j++)
{
temp=0;
for(m=0;m<=20;m++)
{ temp+=Gt[i][m]*G[m][j]; }
GtG[i][j]=temp;
}
}

for(l=0;l<=20;l++)
{

for(i=0;i<=k;i++)
{
for(j=0;j<=k;j++)
{ GtG1[i][j]=GtG[i][j]; }
}

for(m=0;m<=k-1;m++)
{
ik=m;
for(i=m;i<=k-1;i++)
{
if(fabs(GtG1[i][m])<fabs(GtG1[i+1][m]))
ik=i+1;
else ;
}
if(ik!=m)
{
for(i=m;i<=k;i++)
{
temp=GtG1[m][i];
GtG1[m][i]=GtG1[ik][i];
GtG1[ik][i]=temp;
}
temp=Gt[m][l];
Gt[m][l]=Gt[ik][l];
Gt[ik][l]=temp;
}

for(i=m+1;i<=k;i++)
{
q=GtG1[i][m]/GtG1[m][m];
for(j=m;j<=k;j++)
{ GtG1[i][j]=GtG1[i][j]-q*GtG1[m][j]; }
Gt[i][l]=Gt[i][l]-q*Gt[m][l];
}
}

D[k][l]=Gt[k][l]/GtG1[k][k];
for(m=k-1;m>=0;m--)
{
temp=0;
for(i=m+1;i<=k;i++)
{ temp+=D[i][l]*GtG1[m][i]; }
D[m][l]=(Gt[m][l]-temp)/GtG1[m][m];
}
}

for(i=0;i<=k;i++)
{
for(j=0;j<=k;j++)
{
temp=0;
for(m=0;m<=20;m++)
{ temp+=A[i][m]*D[j][m]; }
C[i][j]=temp;
}
}

sigma=0;/*归零，开始计算sigma值*/

for(i=0;i<=10;i++)
{
for(j=0;j<=20;j++)
{
p[i][j]=0;
for(i1=0;i1<=k;i1++)
{
for(j1=0;j1<=k;j1++)
{ p[i][j]+=C[i1][j1]*pow(x[i],i1)*pow(y[j],j1); }
}
sigma+=pow(p[i][j]-z[i][j],2);
}
}
printf("k=%d sigma=%.12e\n",k,sigma);
k=k+1;
}
k--;
printf("达到精度要求时的k和sigma值：\n");
printf("k=%d sigma=%.12e\n",k,sigma);
kz=k; /*定义为全局变量，便于duibi()调用*/
}

void duibi() /*对比𝑓和p逼近效果的子程序*/
{ int i,j,i1,j1;
double p[8][5];

for(i=0;i<=7;i++)
{
for(j=0;j<=4;j++)
{
x[i]=0.1*(i+1);
y[j]=0.5+0.2*(j+1);}} /*重新输入节点*/

Newton(x,y);
fpeccz(t,u); /*求解𝑓（x*，y*）*/

for(i=0;i<=7;i++) /*求解p（x*，y*）*/
{
for(j=0;j<=4;j++)
{
p[i][j]=0;
for(i1=0;i1<=kz;i1++)
{
for(j1=0;j1<=kz;j1++)
{ p[i][j]+=C[i1][j1]*pow(x[i],i1)*pow(y[j],j1); }
}
printf("x[%d]=%.6f y[%d]=%.6f\n",i+1,x[i],j+1,y[j]);
printf("f(x[%d],y[%d])=%.12e\n",i+1,j+1,z[i][j]);
printf("p(x[%d],y[%d])=%.12e\n",i+1,j+1,p[i][j]);
printf("σ=%.12e\n",z[i][j]-p[i][j]);
}/*数表xi*,yi*,𝑓 (xi*,yi*),p（xi*,yi*）*/
}
}

void main()
{ int i,j;
for(i=0;i<=10;i++)
{
for(j=0;j<=20;j++)
{
x[i]=0.08*i;
y[j]=0.5+0.05*j;}} /*输入节点*/
Newton(x,y);
fpeccz(t,u);

for(i=0;i<=10;i++)
{
for(j=0;j<=20;j++)
{
printf("x[%d]=%.6f y[%d]=%.6f z[%d][%d]=%.12e \n",i,x[i],j,y[j],i,j,z[i][j]);
}} /*数表：xi,yi, 𝑓 (xi,yi)*/

qmnh();
for(i=0;i<=kz;i++)
{
for(j=0;j<=kz;j++)
{
printf("C[%d][%d]=%.12e \n",i,j,C[i][j]);
}}/*数表：C[r][s]*/

duibi();
}
三、运行结果输出
1、数表数表：xi,yi,𝑓(xi,yi)
x[0]=0.000000 y[0]=0.500000 z[0][0]=4.465040184807e-001
x[0]=0.000000 y[1]=0.550000 z[0][1]=3.246832629277e-001
x[0]=0.000000 y[2]=0.600000 z[0][2]=2.101596866827e-001
x[0]=0.000000 y[3]=0.650000 z[0][3]=1.030436083160e-001
x[0]=0.000000 y[4]=0.700000 z[0][4]=3.401895562675e-003
x[0]=0.000000 y[5]=0.750000 z[0][5]=-8.873581363800e-002
x[0]=0.000000 y[6]=0.800000 z[0][6]=-1.733716327497e-001
x[0]=0.000000 y[7]=0.850000 z[0][7]=-2.505346114666e-001
x[0]=0.000000 y[8]=0.900000 z[0][8]=-3.202765063876e-001
x[0]=0.000000 y[9]=0.950000 z[0][9]=-3.826680661097e-001
x[0]=0.000000 y[10]=1.000000 z[0][10]=-4.377957667384e-001
x[0]=0.000000 y[11]=1.050000 z[0][11]=-4.857589414438e-001
x[0]=0.000000 y[12]=1.100000 z[0][12]=-5.266672548835e-001
x[0]=0.000000 y[13]=1.150000 z[0][13]=-5.606384797965e-001
x[0]=0.000000 y[14]=1.200000 z[0][14]=-5.877965387677e-001
x[0]=0.000000 y[15]=1.250000 z[0][15]=-6.082697790899e-001
x[0]=0.000000 y[16]=1.300000 z[0][16]=-6.221894528764e-001
x[0]=0.000000 y[17]=1.350000 z[0][17]=-6.296883781856e-001
x[0]=0.000000 y[18]=1.400000 z[0][18]=-6.308997600028e-001
x[0]=0.000000 y[19]=1.450000 z[0][19]=-6.259561525454e-001
x[0]=0.000000 y[20]=1.500000 z[0][20]=-6.149885466094e-001
x[1]=0.080000 y[0]=0.500000 z[1][0]=6.380152265113e-001
x[1]=0.080000 y[1]=0.550000 z[1][1]=5.066117551467e-001
x[1]=0.080000 y[2]=0.600000 z[1][2]=3.821763692774e-001
x[1]=0.080000 y[3]=0.650000 z[1][3]=2.648634911537e-001
x[1]=0.080000 y[4]=0.700000 z[1][4]=1.547802002848e-001
x[1]=0.080000 y[5]=0.750000 z[1][5]=5.199268349094e-002
x[1]=0.080000 y[6]=0.800000 z[1][6]=-4.346804020490e-002
x[1]=0.080000 y[7]=0.850000 z[1][7]=-1.316010567885e-001
x[1]=0.080000 y[8]=0.900000 z[1][8]=-2.124310883088e-001
x[1]=0.080000 y[9]=0.950000 z[1][9]=-2.860045510580e-001
x[1]=0.080000 y[10]=1.000000 z[1][10]=-3.523860789794e-001
x[1]=0.080000 y[11]=1.050000 z[1][11]=-4.116554565222e-001
x[1]=0.080000 y[12]=1.100000 z[1][12]=-4.639049115188e-001
x[1]=0.080000 y[13]=1.150000 z[1][13]=-5.092367247005e-001
x[1]=0.080000 y[14]=1.200000 z[1][14]=-5.477611179623e-001
x[1]=0.080000 y[15]=1.250000 z[1][15]=-5.795943883391e-001
x[1]=0.080000 y[16]=1.300000 z[1][16]=-6.048572588895e-001
x[1]=0.080000 y[17]=1.350000 z[1][17]=-6.236734213318e-001
x[1]=0.080000 y[18]=1.400000 z[1][18]=-6.361682484133e-001
x[1]=0.080000 y[19]=1.450000 z[1][19]=-6.424676566901e-001
x[1]=0.080000 y[20]=1.500000 z[1][20]=-6.426971026996e-001
x[2]=0.160000 y[0]=0.500000 z[2][0]=8.400813957666e-001
x[2]=0.160000 y[1]=0.550000 z[2][1]=6.997641656732e-001
x[2]=0.160000 y[2]=0.600000 z[2][2]=5.660614423517e-001
x[2]=0.160000 y[3]=0.650000 z[2][3]=4.391716081176e-001
x[2]=0.160000 y[4]=0.700000 z[2][4]=3.192421380408e-001
x[2]=0.160000 y[5]=0.750000 z[2][5]=2.063761923874e-001
x[2]=0.160000 y[6]=0.800000 z[2][6]=1.006385238914e-001
x[2]=0.160000 y[7]=0.850000 z[2][7]=2.060740067837e-003
x[2]=0.160000 y[8]=0.900000 z[2][8]=-8.935402476698e-002
x[2]=0.160000 y[9]=0.950000 z[2][9]=-1.736269688648e-001
x[2]=0.160000 y[10]=1.000000 z[2][10]=-2.507999561599e-001
x[2]=0.160000 y[11]=1.050000 z[2][11]=-3.209322694446e-001
x[2]=0.160000 y[12]=1.100000 z[2][12]=-3.840977350046e-001
x