北航研究生数值分析上机作业 三 报告+所有程序大全.docx
《北航研究生数值分析上机作业 三 报告+所有程序大全.docx》由会员分享，可在线阅读，更多相关《北航研究生数值分析上机作业 三 报告+所有程序大全.docx（31页珍藏版）》请在冰豆网上搜索。
北航研究生数值分析上机作业三报告+所有程序大全
数值分析上机作业3——求解非线性方程组以及二元函数的插值拟合
1.算法设计
对于全部的插值节点
，带入非线性方程组中，用Newton迭代法解非线性方程组，得到
。
对
，在二维数表中进行插值，采用分片双二次插值法。
插值过程中，先选择分片区域的中心节点，在数表中的列记为
，行记为
，中心节点记为
，生成向量
，
，
，
，
同理，生成向量
，
记数表中以分片区域中心节点为中心的3×3的矩阵为
对于
插值结果为
。
在拟合
时，需要计算
，令
，计算
时，根据对称性只需要计算对角线元素和对角线以上元素即可，节省运算时间。
于是
，用选主元的LU分解法求解
，再计算
，这里
只需要按行取元素进行运算即可，故不需要进行转置运算。
从1到9依次增加，计算
的值，当
时，得到达到精度的最小的
。
2.打印输出结果
f（x,y）:
4.46504018480e-0013.24683262927e-0012.10159686683e-0011.03043608316e-0013.40189556266e-003-8.87358136384e-002
6.38015226510e-0015.06611755146e-0013.82176369277e-0012.64863491154e-0011.54780200285e-0015.19926834902e-002
8.40081395765e-0016.99764165673e-0015.66061442351e-0014.39171608118e-0013.19242138041e-0012.06376192387e-001
1.05151509180e+0009.02927430830e-0017.60580266859e-0016.24715198145e-0014.95519756001e-0013.73134042774e-001
1.27124675148e+0001.11500201815e+0009.64607727215e-0018.20347369475e-0016.82447678179e-0015.51085208597e-001
1.49832105248e+0001.33499863207e+0001.17712512374e+0001.02502405502e+0008.78960023174e-0017.39145108703e-001
1.73189274038e+0001.56203457721e+0001.39721691821e+0001.23780100674e+0001.08408753268e+0009.36322772315e-001
1.97122178640e+0001.79532959950e+0001.62406711323e+0001.45783058271e+0001.29695464975e+0001.14171810545e+000
k=1,sigma=3.22090897363e+000
k=2,sigma=4.65996003320e-003
k=3,sigma=1.72117537926e-004
k=4,sigma=3.30953430190e-006
k=5,sigma=2.54137773513e-008
C_rs:
2.02123036395e+000-3.66842591519e+0007.09246688452e-0018.48607444215e-001-4.15898479841e-0016.74322022261e-002
3.19192646165e+000-7.41209812962e-001-2.69690653026e+0001.63095275395e+000-4.84601977266e-0016.05908514261e-002
2.56706343343e-0011.58096413206e+000-4.65701259544e-001-7.89186706497e-0021.00853116187e-001-2.07687560816e-002
-2.68608304872e-001-7.33963449843e-0011.08429601112e+000-8.15635815143e-0013.07285137544e-001-4.68489486618e-002
2.16521800059e-001-1.73026852965e-001-8.41324310602e-0022.55736987891e-001-1.47683427939e-0012.77711894906e-002
-5.54328606191e-0021.40518220408e-001-1.30388672239e-0013.44966421960e-0026.95959872154e-003-3.30022941240e-003
f（x*,y*）:
0.194720-0.183037-0.445498-0.597567-0.646460
0.405979-0.022516-0.338221-0.544438-0.647361
0.6347770.158801-0.207366-0.465358-0.620271
0.8789600.358651-0.055253-0.362680-0.567565
1.1366110.5749800.115992-0.238568-0.491434
1.4060420.8059410.304429-0.095016-0.393902
1.6857841.0498810.5082940.066149-0.276834
1.9745751.3053350.7259920.243254-0.141950
p（x*,y*）:
0.194730-0.183042-0.445500-0.597559-0.646446
0.405990-0.022521-0.338224-0.544430-0.647348
0.6347870.158796-0.207369-0.465350-0.620257
0.8789700.358646-0.055255-0.362671-0.567551
1.1366200.5749760.115989-0.238560-0.491421
1.4060510.8059370.304426-0.095009-0.393890
1.6857911.0498780.5082910.066156-0.276822
1.9745811.3053320.7259890.243261-0.141939
对
插值的数值情况：
的数值情况（z坐标为10-7）：
3.源程序
1）主函数
#include
#include
#definelength_x11//x向量长度
#definelength_y21//y向量长度
#definex_length8//x*向量长度
#definey_length5//y*向量长度
externdoubleinterpolation（doublett,doubleuu,double**table）;//输入一点，输出一点插值
externdoublerow_multi_sum（double**B,inti,intj,intlength）;//矩阵列相乘
intmain（）
{
doublett=0;
doubleuu=0;
doubletable[6][6]={0};
doubleff=0;
doublexx[length_x]={0};
doubleyy[length_y]={0};
doublesum=0.0;
doubleU[length_x][length_y]={0};
doubleB[length_x][10]={0};
doubleG[length_y][10]={0};
doubleW[10][10]={0};
doubleT[10][10]={0};
doubleM[10][10]={0};
doubleC[10][10]={0};
doubleF[10][length_y]={0};
doubleE[10][10]={0};
doubleXX[10][10]={0};
//doubleff[length_x][length_y]={0};//记录f阵
doublepp[length_x][length_y]={0};//记录p阵
intr=0;
ints=0;
intii,jj,kk;
intk_max=0;
FILE*fd;
fd=fopen（"tmptmp.txt","w"）;
ii=0;
jj=0;
kk=0;
table[0][0]=-0.5;
table[0][1]=-0.34;
table[0][2]=0.14;
table[0][3]=0.94;
table[0][4]=2.06;
table[0][5]=3.5;
table[1][0]=-0.42;
table[1][1]=-0.5;
table[1][2]=-0.26;
table[1][3]=0.3;
table[1][4]=1.18;
table[1][5]=2.38;
table[2][0]=-0.18;
table[2][1]=-0.5;
table[2][2]=-0.5;
table[2][3]=-0.18;
table[2][4]=0.46;
table[2][5]=1.42;
table[3][0]=0.22;
table[3][1]=-0.34;
table[3][2]=-0.58;
table[3][3]=-0.5;
table[3][4]=-0.1;
table[3][5]=0.62;
table[4][0]=0.78;
table[4][1]=-0.02;
table[4][2]=-0.5;
table[4][3]=-0.66;
table[4][4]=-0.5;
table[4][5]=-0.02;
table[5][0]=1.5;
table[5][1]=0.46;
table[5][2]=-0.26;
table[5][3]=-0.66;
table[5][4]=-0.74;
table[5][5]=-0.5;
for（ii=0;iixx[ii]=0.08*ii;
for（ii=0;iiyy[ii]=0.5+0.05*ii;
for（ii=0;ii{
for（jj=0;jj{
ite_equations（&tt,&uu,xx[ii],yy[jj]）;//牛顿迭代解方程，输入一点（x,y）输出一点（t,u）
U[ii][jj]=（double）interpolation（tt,uu,（double**）table）;//输入一点，输出一点插值
}
}
fprintf（fd,"\nf（x,y）:
\n"）;
for（ii=0;ii{
for（jj=0;jj<=y_length;jj++）
{
fprintf（fd,"%.11e",U[ii][jj]）;
}
fprintf（fd,"\n"）;
}
for（k_max=1;k_max<=6;k_max++）
{
for（kk=0;kk<（k_max+1）;kk++）//生成B阵
{
for（ii=0;ii{
B[ii][kk]=pow（xx[ii],kk）;
}
}
for（kk=0;kk<（k_max+1）;kk++）
{
for（ii=0;ii{
G[ii][kk]=pow（yy[ii],kk）;
}
}
for（ii=0;ii<=k_max;ii++）//BT*B
{
for（jj=0;jj<=k_max;jj++）
{
W[ii][jj]=row_multi_sum（（double**）B,ii,jj,length_x）;
}
}
for（ii=0;ii<=k_max;ii++）//GT*G
{
for（jj=0;jj<=k_max;jj++）
{
M[ii][jj]=row_multi_sum（（double**）G,ii,jj,length_y）;
}
}
for（ii=0;ii<=k_max;ii++）//F=B'*U
for（jj=0;jj{
F[ii][jj]=0;
for（kk=0;kkF[ii][jj]=F[ii][jj]+B[kk][ii]*U[kk][jj];
}
for（ii=0;ii<=k_max;ii++）//E=F*G
for（jj=0;jj<=k_max;jj++）
{
E[ii][jj]=0;
for（kk=0;kkE[ii][jj]=E[ii][jj]+F[ii][kk]*G[kk][jj];
}
Matrix_LU（（double**）W,（double**）E,（double**）M,（double**）C,k_max,0）;//（double**）XX,
for（ii=0;iifor（jj=0;jj{
sum=0;
for（r=0;r<=k_max;r++）
for（s=0;s<=k_max;s++）
sum=sum+C[r][s]*pow（xx[ii],r）*pow（yy[jj],s）;
pp[ii][jj]=sum;
}
sum=0;
for（ii=0;iifor（jj=0;jj{
sum=sum+（pp[ii][jj]-U[ii][jj]）*（pp[ii][jj]-U[ii][jj]）;
}
fprintf（fd,"\nk=%d,sigma=%.11e\n",k_max,sum）;
if（sum<=1e-7）
{
fprintf（fd,"\nC_rs:
\n"）;
for（r=0;r<=k_max;r++）
{
for（s=0;s<=k_max;s++）
{
fprintf（fd,"%.11e",C[r][s]）;
}
fprintf（fd,"\n"）;
}
break;
}
}
for（ii=0;iixx[ii]=0.1*（ii+1）;
for（ii=0;iiyy[ii]=0.5+0.2*（ii+1）;
for（ii=0;ii{
for（jj=0;jj{
ite_equations（&tt,&uu,xx[ii],yy[jj]）;//牛顿迭代解方程，输入一点（x,y）输出一点（t,u）
U[ii][jj]=（double）interpolation（tt,uu,（double**）table）;//输入一点，输出一点插值
}
}
fprintf（fd,"\nf（x*,y*）:
\n"）;
for（ii=0;ii{
for（jj=0;jj{
fprintf（fd,"%-10f",U[ii][jj]）;
}
fprintf（fd,"\n"）;
}
for（ii=0;ii{
for（jj=0;jj{
sum=0;
for（r=0;r<=k_max;r++）
for（s=0;s<=k_max;s++）
sum=sum+C[r][s]*pow（xx[ii],r）*pow（yy[jj],s）;
pp[ii][jj]=sum;
}
}
fprintf（fd,"\np（x*,y*）:
\n"）;
for（ii=0;ii{
for（jj=0;jj{
fprintf（fd,"%-10f",pp[ii][jj]）;
}
fprintf（fd,"\n"）;
}
fclose（fd）;
}
2）矩阵2列相乘
#include
#include
#definesize10
#defineTTT（a,b）*（（double*）TTT+（a）*（size）+b）//LU的宏定义，以方便操作二维数组
doublerow_multi_sum（double**B,inti,intj,intlength）//矩阵列相乘
{
doublesum=0;
double**TTT;
intkk=0;
TTT=B;
for（kk=0;kk{
sum=sum+（TTT（kk,i））*（TTT（kk,j））;
}
returnsum;
}
3）9点分片双二次插值
#include
#include
#definesize6
#definetable（a,b）*（（double*）table+（a）*（size）+b）//LU的宏定义，以方便操作二维数组
doubleinterpolation（doublett,doubleuu,double**table）//输入一点，输出一点插值
{
intii=0;
intjj=0;
intsection_t=0;
intsection_u=0;
doubleanswer=0;
doubletemp=0;
doubletemp_u[3]={0};
doubletemp_t[3]={0};
doublemy_t[6]={0};
doublemy_u[6]={0};
for（ii=1;ii<6;ii++）
{
my_u[ii]=my_u[ii-1]+0.4;
my_t[ii]=my_t[ii-1]+0.2;
}
//判断输入点在哪个区域
//判断u用哪段插值
if（（uu>=0）&&（uu<0.6））
section_u=1;
elseif（（uu>=0.6）&&（uu<1.0））
section_u=2;
elseif（（uu>=1.0）&&（uu<1.4））
section_u=3;
elseif（（uu>=1.4）&&（uu<=2））
section_u=4;
else
;
//判断t用哪段插值
if（（tt>=0）&&（tt<0.3））
section_t=1;
elseif（（tt>=0.3）&&（tt<0.5））
section_t=2;
elseif（（tt>=0.5）&&（tt<0.7））
section_t=3;
elseif（（tt>=0.7）&&（tt<=1.0））
section_t=4;
else
;
temp_u[0]=（uu-my_u[section_u]）*（uu-my_u[section_u+1]）/（（my_u[section_u-1]-my_u[section_u]）*（my_u[section_u-1]-my_u[section_u+1]））;
temp_u[1]=（uu-my_u[section_u-1]）*（uu-my_u[section_u+1]）/（（my_u[section_u]-my_u[section_u-1]）*（my_u[section_u]-my_u[section_u+1]））;
temp_u[2]=（uu-my_u[section_u-1]）*（uu-my_u[section_u]）/（（my_u[section_u+1]-my_u[section_u-1]）*（my_u[section_u+1]-my_u[section_u]））;
temp_t[0]=（tt-my_t[section_t]）*（tt-my_t[section_t+1]）/（（my_t[section_t-1]-my_t[section_t]）*（my_t[section_t-1]-my_t[section_t+1]））;
temp_t[1]=（tt-my_t[section_t-1]）*（tt-my_t[section_t+1]）/（（my_t[section_t]-my_t[section_t-1]）*（my_t[section_t]-my_t[section_t+1]））;
temp_t[2]=（tt-my_t[section_t-1]）*（tt-my_t[section_t]）/（（my_t[section_t+1]-my_t[section_t-1]）*（my_t[section_t+1]-my_t[section_t]））;
answer=0;
for（ii=0;ii<3;ii++）
{
answer=answer+temp_t[ii]*（temp_u[0]*table（（section_t-1+ii）,section_u-1）+temp_u[1]*table（（s
升级会员 