北航数值分析实习第三题.docx
《北航数值分析实习第三题.docx》由会员分享，可在线阅读，更多相关《北航数值分析实习第三题.docx（40页珍藏版）》请在冰豆网上搜索。
北航数值分析实习第三题
数值分析计算实习报告
第三题
所在班级
A1班
学生姓名
学生学号
2015年11月
1算法设计方案
1.1求对应于
，
的
1.1.1通过求解非线性方程组获得x，y与t，u的关系
题目中给出的方程组建立了t，u和x，y的联系，而数表建立了t，u与z的联系。
为了求得x，y与z的联系，我们先求解方程组，获得t，u和x，y的联系。
将
和
代入非线性方程组中，用Newton法解出
和
。
1.1.2通过分片二次代数插值求得t，u和z的关系，进而得到x，y和z的关系
上一步，已经获得对应于
，
的
，
，现在只要再获得对应于
，
的z，就间接求得了对应于
，
的
。
根据题目给出的数表，进行分片二次代数插值，之后带入
，
，获得了
，它对应于
，
。
1.2曲面插值
上一步建立了数表
，利用该数表进行曲面插值，本质上就是求得系数矩阵
，就可获得
。
k从1开始尝试，每个k都进行曲面插值，满足精度要求时停止计算。
1.3观察
逼近
的效果
这里新给了点
，
和
将其重复上面步骤，得到与新的插值节点
对应的
。
再将
带入上一步求得的
即可得到。
比较时可以比较相对误差和绝对误差。
1.4补充—计算逆矩阵
计算系数矩阵
时，需要求解一些矩阵的逆矩阵。
这里采用列主元Gauss消去法来求解，即将原矩阵初等变换为单位阵，同样的变换可以把单位阵变为原矩阵的逆矩阵。
2C++程序
#include"stdio.h"
#include"math.h"
doublete[11][21]={0};
doubleue[11][21]={0};
doubleze[11][21]={0};
doublex[11]={0};
doubley[21]={0};
doubleinv[10][10]={0};
doubleC[10][10]={0};
doublers=0;
doublefanshu（double*p）//求向量的无穷范数
{
inti=0;
doublemax=fabs（p[1]）;
for（i=1;i<=4;i++）
{if（fabs（p[i]）>max）
max=fabs（p[i]）;
}
return（max）;
}
voidinverse（doubleX[10][10],intN）//采用列主元Gauss消去法求逆矩阵
{doublematrix[10][10];
doubleI[10][10]={0};
inti,j;
intk,ik,cnt,count;
doublem;
doubletemp;
for（i=1;i<=N;i++）
for（j=1;j<=N;j++）
matrix[i][j]=X[i][j];
for（i=1;i<=N;i++）
I[i][i]=1;//单位阵
for（k=1;k<=N-1;k++）
{
ik=k;
for（cnt=k;cnt<=N;cnt++）//选最大元素行号
{
if（fabs（matrix[cnt][k]）>fabs（matrix[ik][k]））
{
ik=cnt;
}
}
for（cnt=1;cnt<=N;cnt++）//交换
{
temp=matrix[k][cnt];
matrix[k][cnt]=matrix[ik][cnt];
matrix[ik][cnt]=temp;
temp=I[k][cnt];
I[k][cnt]=I[ik][cnt];
I[ik][cnt]=temp;
}
for（cnt=k+1;cnt<=N;cnt++）//消去了
{
m=matrix[cnt][k]/matrix[k][k];
for（count=1;count<=N;count++）
{
matrix[cnt][count]=matrix[cnt][count]-m*matrix[k][count];
I[cnt][count]=I[cnt][count]-m*I[k][count];
}
}
}
for（i=1;i<=N;i++）//把上对角矩阵化为单位阵的相同步骤就可把经过变形的I变成原矩阵的逆矩阵
{
inv[N][i]=I[N][i]/matrix[N][N];
for（k=N-1;k>=1;k--）
{
temp=0;
for（cnt=k+1;cnt<=N;cnt++）
temp+=matrix[k][cnt]*inv[cnt][i];
inv[k][i]=（I[k][i]-temp）/matrix[k][k];
}
}
}
voidnewton（doublex[11],doubley[11]）//牛顿法求解非线性方程组
{
doublet=0;
doubleu=0;
doublew=0;
doublev=0;
doubledF[5][5]={0};//F的导数
doubleF[5]={0};
doubled[5]={0};
inti,j,k;
intik;
intcnt,count;
doubletemp=0;
doublem=0;
doublejie[4]={0};
for（i=0;i<11;i++）
for（j=0;j<21;j++）
{t=1;
u=1;
w=1;
v=1;
do
{dF[1][1]=-0.5*sin（t）;
dF[1][2]=1;
dF[1][3]=1;
dF[1][4]=1;
dF[2][1]=1;
dF[2][2]=0.5*cos（u）;
dF[2][3]=1;
dF[2][4]=1;
dF[3][1]=0.5;
dF[3][2]=1;
dF[3][3]=-sin（v）;
dF[3][4]=1;
dF[4][1]=1;
dF[4][2]=0.5;
dF[4][3]=1;
dF[4][4]=cos（w）;
F[1]=-（0.5*cos（t）+u+v+w-x[i]-2.67）;//这里实际是-F
F[2]=-（t+0.5*sin（u）+v+w-y[j]-1.07）;
F[3]=-（0.5*t+u+cos（v）+w-x[i]-3.74）;
F[4]=-（t+0.5*u+v+sin（w）-y[j]-0.79）;
for（k=1;k<=3;k++）//列主元Gauss消去法求解dF*△x=-F
{
ik=k;
for（cnt=k;cnt<=4;cnt++）
{
if（fabs（dF[cnt][k]）>fabs（dF[ik][k]））
{
ik=cnt;
}
}
temp=F[k];
F[k]=F[ik];
F[ik]=temp;
for（cnt=k;cnt<=4;cnt++）
{
temp=dF[k][cnt];
dF[k][cnt]=dF[ik][cnt];
dF[ik][cnt]=temp;
}
for（cnt=k+1;cnt<=4;cnt++）
{
m=dF[cnt][k]/dF[k][k];
for（count=k+1;count<=4;count++）
dF[cnt][count]=dF[cnt][count]-m*dF[k][count];
F[cnt]=F[cnt]-m*F[k];
}
}
d[4]=F[4]/dF[4][4];
for（k=4-1;k>=1;k--）
{
temp=0;
for（cnt=k+1;cnt<=4;cnt++）
temp+=dF[k][cnt]*d[cnt];
d[k]=（F[k]-temp）/dF[k][k];
}
t=t+d[1];
u=u+d[2];
v=v+d[3];
w=w+d[4];
te[i][j]=t;
ue[i][j]=u;
jie[1]=t;
jie[2]=u;
jie[3]=v;
jie[4]=w;
}while（fanshu（d）/fanshu（jie）>1e-12）;//引用求范数的函数
}
}
voideryuanchazhi（doublete[11][21],doubleue[11][21]）//分片二次代数插值求z=f（x,y）
{inti,j,ii,jj,m;
intr[11][21]={0};
intk[11][21]={0};
doubletemp=0;
doubleu[6]={0,0.4,0.8,1.2,1.6,2};
doublet[6]={0,0.2,0.4,0.6,0.8,1.0};
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}};
for（i=0;i<=10;i++）
for（j=0;j<=20;j++）
{
if（te[i][j]<=0.3）//选择t的插值点，考察是边界还是中间的
r[i][j]=1;
elseif（te[i][j]>=0.7）
r[i][j]=4;
else
{for（m=2;m<=3;m++）
if（（te[i][j]>0.2*m-0.1）&&（te[i][j]<0.2*m+0.1））
r[i][j]=m;
}
if（ue[i][j]<=0.6）//选择u的插值点，考察是边界还是中间的
k[i][j]=1;
elseif（ue[i][j]>=1.4）
k[i][j]=4;
else
{for（m=2;m<=3;m++）
if（（ue[i][j]>0.4*m-0.2）&&（ue[i][j]<0.4*m+0.2））
k[i][j]=m;
}
ze[i][j]=0;
for（ii=k[i][j]-1;ii<=k[i][j]+1;ii++）
for（jj=r[i][j]-1;jj<=r[i][j]+1;jj++）//这两个循环是求和层次的
{
temp=1;
for（m=k[i][j]-1;m<=k[i][j]+1;m++）//这两个循环是求l的求积层次的
if（m!
=ii）
temp*=（te[i][j]-t[m]）/（t[ii]-t[m]）;
for（m=r[i][j]-1;m<=r[i][j]+1;m++）
if（m!
=jj）
temp*=（ue[i][j]-u[m]）/（u[jj]-u[m]）;
ze[i][j]+=temp*z[ii][jj];
}
}
FILE*fp=fopen（"shubiao-xyf.txt","a"）;
for（i=0;i<=10;i++）
for（j=0;j<=20;j++）
fprintf（fp,"x[%d]=%f\ty[%d]=%f\tf（x,y）=%13.11e\n",i,x[i],j,y[j],ze[i][j]）;
fclose（fp）;
}
voidqumianchazhi（）//曲面插值，求p（x,y）
{inti,j;
intii,jj;
intk,t;
doubletemp=0;
doubletemp1=0;
doublesigma=0;
doubleB[11+1][10]={0};
doubleG[21+1][10]={0};
doubleBB[10][10]={0};//B转置*B
doubleGG[10][10]={0};//G转置*G
doubleBBN[10][10]={0};//（B转置*B）的逆
doubleGGN[10][10]={0};//（G转置*G）的逆
doubleBBB[10][12]={0};//[（B转置*B）的逆]*B的转置
doubleBBBU[10][22]={0};//{[（B转置*B）的逆]*B的转置}*U
doubleBBBUG[10][10]={0};//{[（B转置*B）的逆]*B的转置}*U*G
doublep=0;
FILE*fp;
fp=fopen（"shubiao-ksic.txt","a"）;
k=0+1;//这里先加上1是为了矩阵运算方便
do{
//首先把系数矩阵C求出来
for（i=1;i<=11;i++）
for（j=1;j<=k;j++）
{
B[i][j]=pow（x[i-1],j-1）;
}
for（i=1;i<=21;i++）
for（j=1;j<=k;j++）
{
G[i][j]=pow（y[i-1],j-1）;
}
for（i=1;i<=k;i++）//求B转置*B
for（j=1;j<=k;j++）
{
temp=0;
for（t=1;t<=11;t++）
temp+=B[t][i]*B[t][j];
BB[i][j]=temp;
}
for（i=1;i<=k;i++）//求G转置*G
for（j=1;j<=k;j++）
{
temp=0;
for（t=1;t<=21;t++）
temp+=G[t][i]*G[t][j];
GG[i][j]=temp;
}
inverse（BB,k）;
for（i=1;i<=k;i++）//（B转置*B）的逆
for（j=1;j<=k;j++）
BBN[i][j]=inv[i][j];
inverse（GG,k）;
for（i=1;i<=k;i++）//（G转置*G）的逆
for（j=1;j<=k;j++）
GGN[i][j]=inv[i][j];
for（i=1;i<=k;i++）//[（B转置*B）的逆]*B的转置
for（j=1;j<=11;j++）
{
temp=0;
for（t=1;t<=k;t++）
temp+=BBN[i][t]*B[j][t];
BBB[i][j]=temp;
}
for（i=1;i<=k;i++）//{[（B转置*B）的逆]*B的转置}*U
for（j=1;j<=21;j++）
{
temp=0;
for（t=1;t<=11;t++）
temp+=BBB[i][t]*ze[t-1][j-1];
BBBU[i][j]=temp;
}
for（i=1;i<=k;i++）//{[（B转置*B）的逆]*B的转置}*U*G
for（j=1;j<=k;j++）
{
temp=0;
for（t=1;t<=21;t++）
temp+=BBBU[i][t]*G[t][j];
BBBUG[i][j]=temp;
}
for（i=1;i<=k;i++）//得到系数矩阵C
for（j=1;j<=k;j++）
{
temp=0;
for（t=1;t<=k;t++）
temp+=BBBUG[i][t]*GGN[t][j];
C[i][j]=temp;
}
//下面判断现在的k对应的C是否满足精度要求
for（i=0;i<11;i++）
for（j=0;j<21;j++）
{temp=0;
for（ii=0;iifor（jj=0;jjtemp+=C[ii+1][jj+1]*pow（x[i],ii）*pow（y[j],jj）;//这里的下标，矩阵下标与数组下标差1
p=temp;
temp1+=（ze[i][j]-p）*（ze[i][j]-p）;
}
sigma=temp1;
temp1=0;
fprintf（fp,"k=%d,σ=%13.11e\n",k-1,sigma）;
k++;
}while（sigma>=1e-7）;
k--;//减去最后一次判断循环条件前加的1
rs=k;//存储k值，用于比较函数中使用
fprintf（fp,"k=%d\n",k-1）;//减去一开始为了矩阵运算加的1
for（i=1;i<=k;i++）
for（j=1;j<=k;j++）
fprintf（fp,"C[%d][%d]=%13.11e\n",i-1,j-1,C[i][j]）;//这里为了匹配矩阵下标和书上所示的系数下标所以减1
}
voidcompare（）//比较f（x*,y*）与p（x*,y*）
{
doublepn[8+1][5+1]={0};//下标都加1，方便以后运算
inti,j;
doubletemp=0;
intii,jj;
for（i=1;i<=8;i++）
x[i]=0.1*i;
for（j=1;j<=5;j++）
y[j]=0.5+0.2*j;
for（i=1;i<=8;i++）
for（j=1;j<=5;j++）
{
temp=0;
for（ii=0;ii<=rs;ii++）
for（jj=0;jj<=rs;jj++）
temp+=C[ii+1][jj+1]*pow（x[i],ii）*pow（y[j],jj）;
pn[i][j]=temp;
}
newton（x,y）;
eryuanchazhi（te,ue）;
FILE*fp;
fp=fopen（"shubiao-compare.txt","a"）;
for（i=1;i<=8;i++）
for（j=1;j<=5;j++）
{
fprintf（fp,"x*[%d]=%f\ty*[%d]=%f\t\n",i,x[i],j,y[j]）;
fprintf（fp,"z*[%d][%d]=%13.11e\tp*[%d][%d]=%13.11e\t\n",i,j,ze[i][j],i,j,pn[i][j]）;
fprintf（fp,"绝对误差e=f*（x[%d],y[%d]）-p*（x[%d],y[%d]）=%13.11e\n",i,j,i,j,ze[i][j]-pn[i][j]）;
fprintf（fp,"相对误差er=e/f*（x[%d],y[%d]）=%13.11e\n",i,j,（ze[i][j]-pn[i][j]）/ze[i][j]）;
}
}
voidmain（）
{inti,j;
for（i=0;i<11;i++）
for（j=0;j<21;j++）
{
x[i]=0.08*i;
y[j]=0.5+0.05*j;
}
newton（x,y）;
eryuanchazhi（te,ue）;
qumianchazhi（）;
compare（）;
}
3
运行结果
3.1数表
x[0]=0.000000y[0]=0.500000f（x,y）=4.46504018480e-001
x[0]=0.000000y[1]=0.550000f（x,y）=3.24683262927e-001
x[0]=0.000000y[2]=0.600000f（x,y）=2.10159686683e-001
x[0]=0.000000y[3]=0.650000f（x,y）=1.03043608316e-001
x[0]=0.000000y[4]=0.700000f（x,y）=3.40189556266e-003
x[0]=0.000000y[5]=0.750000f（x,y）=-8.87358136380e-002
x[0]=0.000000y[6]=0.800000f（x,y）=-1.73371632750e-001
x[0]=0.000000y[7]=0.850000f（x,y）=-2.50534611467e-001
x[0]=0.000000y[8]=0.900000f（x,y）=-3.20276506388e-001
x[0]=0.000000y[9]=0.950000f（x,y）=-3.82668066110e-001
x[0]=0.000000y[10]=1.000000f（x,y）=-4.37795
升级会员 