北航数值分析实习第三题文档格式.docx
《北航数值分析实习第三题文档格式.docx》由会员分享，可在线阅读，更多相关《北航数值分析实习第三题文档格式.docx（30页珍藏版）》请在冰豆网上搜索。
比较时可以比较相对误差和绝对误差。
1.4补充—计算逆矩阵
计算系数矩阵时，需要求解一些矩阵的逆矩阵。
这里采用列主元Gauss消去法来求解，即将原矩阵初等变换为单位阵，同样的变换可以把单位阵变为原矩阵的逆矩阵。
2C++程序
#include"
stdio.h"
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;
=N;
for（j=1;
j<
j++）
matrix[i][j]=X[i][j];
I[i][i]=1;
//单位阵
for（k=1;
k<
=N-1;
k++）
{
ik=k;
for（cnt=k;
cnt<
cnt++）//选最大元素行号
{
if（fabs（matrix[cnt][k]）>
fabs（matrix[ik][k]））
{
ik=cnt;
}
}
for（cnt=1;
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++）//消去了
m=matrix[cnt][k]/matrix[k][k];
for（count=1;
count<
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++）//把上对角矩阵化为单位阵的相同步骤就可把经过变形的I变成原矩阵的逆矩阵
{
inv[N][i]=I[N][i]/matrix[N][N];
for（k=N-1;
k>
=1;
k--）
temp=0;
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;
11;
for（j=0;
21;
{ 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;
=3;
k++）//列主元Gauss消去法求解dF*△x=-F
ik=k;
for（cnt=k;
{
if（fabs（dF[cnt][k]）>
fabs（dF[ik][k]））
{
ik=cnt;
}
}
temp=F[k];
F[k]=F[ik];
F[ik]=temp;
for（cnt=k;
temp=dF[k][cnt];
dF[k][cnt]=dF[ik][cnt];
dF[ik][cnt]=temp;
for（cnt=k+1;
m=dF[cnt][k]/dF[k][k];
for（count=k+1;
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;
temp=0;
for（cnt=k+1;
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};
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}};
=10;
=20;
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<
m++）
升级会员 