北航研究生数值分析上机作业 三 报告+所有程序大全.docx
- 文档编号:8155805
- 上传时间:2023-01-29
- 格式:DOCX
- 页数:31
- 大小:153.58KB
北航研究生数值分析上机作业 三 报告+所有程序大全.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;ii xx[ii]=0.08*ii; for(ii=0;ii yy[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;kk F[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;kk E[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;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; } sum=0; for(ii=0;ii for(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;ii xx[ii]=0.1*(ii+1); for(ii=0;ii yy[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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航研究生数值分析上机作业 报告+所有程序大全 北航 研究生 数值 分析 上机 作业 报告 所有 程序 大全