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

北航数值分析第三次大作业————————————————————————————————作者：————————————————————————————————日期：数值分析第三次大作业一、算法的设计方案：（一）、总体方案设计：(1）解非线性方程组。将给定的(,)iixy当作已知量代入题目给定的非线性方程组,求得与(,)iixy相对应的数组t[i][j],ｕ［i][j］。(２)分片二次代数插值。通过分片二次代数插值运算，得到与数组ｔ[1１][21],u[１１][21］]对应的数组ｚ[11］［2１],得到二元函数z＝(,)iifxy。(３）曲面拟合。利用ｘ[i］,y［j],z[1１］[21］建立二维函数表,再根据精度的要求选择适当k值，并得到曲面拟合的系数矩阵C[r］[s]。（4)观察和(,)iipxy的逼近效果。观察逼近效果只需要重复上面(1)和(2）的过程,得到与新的插值节点(,)iixy对应的(,)iifxy,再与对应的(,)iipxy比较即可，这里求解(,)iipxy可以直接使用(３）中的C[r][s]和ｋ。（二)具体算法设计：(1)解非线性方程组牛顿法解方程组()0Fx的解*x，可采用如下算法：１)在*x附近选取(0)xD，给定精度水平0和最大迭代次数Ｍ。２)对于0,1,kM执行①计算()()kFx和()()kFx。②求解关于()kx的线性方程组()()()()()kkkFxxFx③若()()kkxx，则取*()kxx，并停止计算；否则转④。④计算(1)()()kkkxxx。⑤若kM,则继续，否则，输出M次迭代不成功的信息,并停止计算。(２)分片双二次插值给定已知数表以及需要插值的节点，进行分片二次插值的算法：设已知数表中的点为：00(0,1,,)(0,1,,)ijxxihinyyjjm，需要插值的节点为(,)xy。１)根据(,)xy选择插值节点(,)ijxy:若12hxx或12nhxx，插值节点对应取1i或1in，若12yy或12nyy，插值节点对应取1j或1im。若,2222,2222iijjhhxxxinyyyjm则选择(,)(1,,1;1,,1)krxykiiirjjj为插值节点。2)计算1111()(1,,1)()(1,,1)itktikttkjtrtjrttrxxlxkiiixxyylyrjjjyy插值多项式的公式为：1111(,)()()(,)jikrkrkirjpxylxlyfxy注：本步进行插值运算的是(,)tu，利用(,)ijxy与(,)tu的对应关系就可以得到z与(,)ijxy的对应关系。（3）曲面拟合根据插值得到的数表,,(,)ijijxyfxy进行曲面拟合的过程：1）根据拟合节点和基底函数写出矩阵Ｂ和G:010000111101()()()()()()()()()kkknnnxxxxxxBxxx010000111101()()()()()()()()()kkkmmmyyyyyyGyyy2）计算11()()TTTCBBBUGGG。在这里，为了简化计算和编程、避免矩阵求逆，记：1()TTABBBU,1()TTDGGG对上面两式进行变形,得到如下两个线性方程组:()TTBBABU,()TTGGDG通过解上述两个线性方程组,则有：TCAD3）对于每一个(,)ijxy,*00(,)()()kkrsijrsijrspxyCxy。4）拟合需要达到的精度条件为:*2700[(,)]10nmijijijpxyu。其中iju对应着插值得到的数表,,(,)ijijxyfxy中(,)ijfxy的值。5）让k逐步增加,每一次重复执行以上几步，直到*2700[(,)]10nmijijijpxyu成立。此时的k值就是要求解最小的ｋ。二、源程序：#ｉnｃｌudestｄio．ｈ＃inclｕｄeioｓｔreaｍ＞#incｌudestdlib.ｈ＃ｉｎcluｄematｈ.h＞#ｉncｌuｄefloat.ｈ＞#incｌude＜ｉomanip＞#dｅfｉｎeEpsilon11e－12/*解线性方程组时近似解向量的精度＊/#dｅfiｎeM20０/＊解线性方程组时的最大迭代次数＊/#definｅN1０／＊求解迭代次数时假设的k的最大值,用于定义包含ｋ的存储空间*/ｖoｉdNewｔｏｎ();/*牛顿法求解非线性方程组子程序*/ｖoidfpeccz(）；/*分片二次代数插值子程序＊/ｖｏidqmｎh();／＊曲面拟合子程序*／voiｄduiｂｉ(）；／*对比𝑓和p逼近效果的子程序*/doｕblex［11］，y[21],t[11］［２1］,u［11][２１];/＊定义全局变量*/dｏubleｚ[1１］［２1]，C［10］[10]；doubleｋｚ；voidＮewｔon(dｏubleｘ[11],doublｅy[２1］)/*牛顿法求解非线性方程组子程序*／{douｂlｅX［４］，dx[４]，Ｆ[4]，dＦ［4][4］,tｅmp,m,fx,ｆX;ｉnti,j,k，l，ｐ,ｉｋ,n;foｒ(i=0;i＝10;i++)｛fｏr（ｊ＝0;j=２0;j＋+）ﻩ｛ﻩX［0]=１；/*选取迭代初始向量,四个分别代表t，u，v，ｗ*/X[1]=1;X[2]=1;X[3]=１;ﻩn=０；ﻩlｏoｐ1:{F［0]=0．５*cｏs（Ｘ[0］)+X［1］+Ｘ［2]+Ｘ[３]－ｘ[i]-2．６7;F[1]=X［0]＋0.5*ｓｉn(Ｘ[１])+X[2]+X[3］－y［j]－1．07;F［2]=0．5*X[0]+X[1]＋cos（X[2])+X[3]－x[i］－３.74；Ｆ[3］=X[０]＋0.5＊X［1]+X[2]＋sin(X[3])-ｙ[j]-0．79;ﻩ/*求解F（x）*/ｄF［0］[０]=－0.５*sin(X[0])；/*求解Ｆ'（x）*/dF[0][1］=１;dF［0][2]=1；dＦ[0][3]=1；dＦ[1]［0］=１;dＦ[1]［1]＝０.5*coｓ(X［1]）；ｄF[１][2］=１;dF[１]［3]=1;dF［2]［0]=0．5;dF[2]［1]=１;dF[2][2]=-ｓin(X[2])；ｄF［2][3]=1;dF[3][0]＝1;ｄＦ[3][1]＝0．5；dF[3][2]=1;dF[３][3］=ｃoｓ(Ｘ[3]);ﻩﻩ/*高斯选主元消去法求解Δx*/ｆoｒ(ｋ=0;k＜3;k＋+)ﻩ{ﻩﻩiｋ=k；ﻩﻩfｏr（l＝k;l=3;l++)ﻩﻩﻩ{if(dF[ik][k］dＦ[l][ｋ]）ﻩﻩﻩﻩik＝l;}／*选主元*/ﻩﻩtemｐ=0;ﻩｔemp=F[ik]；Ｆ[ik］=F[k];ﻩＦ[k]=temｐ;ﻩﻩfｏr(l=k;l=3；l++)ﻩ{temp=0;ﻩﻩｔeｍp=dF[ｉk］[ｌ］；ﻩｄF［ik][l]=dF[k][l];dF[k][l］=temp;ﻩﻩﻩﻩ｝for（l=k+1；l=３;l+＋)ﻩ｛ﻩm=dF[l][k］/dF[k][ｋ］；F[l］=F［l]-ｍ*Ｆ[k];ﻩｆｏr（ｐ＝k+1;p=３;p+＋){dF［ｌ][ｐ］=dF［l][p］－ｍ*dＦ[ｋ]［p];}ﻩ}／*消去过程*/ﻩ}dx[3]=-Ｆ[3]／ｄF[3][3];ﻩfｏｒ(k=2;ｋ＝0;ｋ－－）ﻩ{ﻩteｍp=０；foｒ(ｌ=k+1；l=3;l++)ﻩﻩ{temp=temp＋ｄＦ[ｋ］[l]＊dｘ[l]/dF［k][ｋ];}ﻩdx[k]=-F［k］／dF［ｋ］[k]-temp;}temｐ=0;foｒ(l=0;l=3;l++)/*求解矩阵范数,用无穷范数*／{ｉf(temｐ＜fabs（dｘ［ｌ])）ｔｅｍp=ｆabs(dｘ[ｌ]);}fx＝ｔｅmp；ｔemp=0;for(l＝0;l＝3;ｌ+＋){if(tempｆａbs(X[l］)）teｍp＝fabｓ（X[l］);}fX=temp;iｆ(ｆaｂｓ(fx/fＸ）Epｓiｌｏｎ1）/*判断()()kkxx是否成立＊/{t[i][j］＝X[０］;u[ｉ］[ｊ］=X[1];ﻩｇotoloop4;}elsｅﻩ{fｏr(l＝0;l＝３;ｌ++){X［l］=X［l]+ｄx[l］;｝ﻩｎ=n＋１;goｔolooｐ3；}｝lｏop3:{if(ｎ＜M）/*判断是否超出规定迭代次数*/ﻩgoｔoloｏｐ１;ﻩelsｅﻩpriｎｔf(迭代不成功\ｎ);gｏtoloop4;}lｏｏp4:{contｉnue；}ﻩ}｝}vｏidｆpeccｚ（ｄｏｕｂlｅｔ[1１][21],ｄｏubleu[11］[2１])／*分片二次代数插值子程序*/{inｔs[11］[21]，r[１1]［21］;intｉ,j，i1,j1,m;ｄｏuｂｌｅz0[6]［６]={{-0.5,－0.３４,０.14，０.９4,2．０６,3.5},{-0.42，-0．5,-0.26,0．3,1.18,２.3８}，{－０.18，－0.5,－0.5，－0.18,0．46,１.42}，{０.22，-0.3４,-0．５8,-0.5,－0.1,0．６２｝,{0.７8，-0.02,-０．5,－0.６6,-0.5,-0.02｝,｛1.5，０.４6,-0.26,-0.66,－0.74,-0.５｝};dｏubｌet0[6］＝{０，0.２,０.4,０．6，0．8,1.0};dｏuｂｌeu0［６]={0,0．４,0.8,1.2,1.６,２．0};douｂleｔemp１,ｔemp2;for(i=0;ｉ＝10;i++)／*选取插值节点*/{fｏｒ(j=0;j＜=2０;j++)ﻩ{ｉf（t[i][j]＝0.3)ｓ[i][j]=1;elseｉf(ｔ[ｉ][j]＞0.７）s[i][j]=4；ｅlseﻩ{for(m=2;m=３;m++)ﻩ{if（（ｔ［i]［j]＞0.２*m－0．1）&&(t[i][j]=0.2*m+0.1））ﻩ｛s[i］［j]＝m；｝ﻩﻩ}}}}foｒ(i=0;i=１0;i++){for(ｊ=0；j＜=２0;j+＋)｛ｉf(u［ｉ][j］=0.６)r[ｉ][ｊ］＝1;ｅlseif（u[i][j]１.4)ｒ［i][j］＝4；ｅlseﻩ{for(m＝2;m＜=3;m++)ﻩ{if（(u[i][j］0.4＊m－0.2)&&(u[ｉ］[ｊ］=0.4*m＋0.2))ﻩ{ｒ[ｉ][ｊ]=m;}ﻩﻩ}｝｝｝fｏr(i=0；i=10;ｉ++）/*插值运算*/{foｒ(j=0；j=2０;ｊ+＋)ﻩ{z[ｉ］[j]=0;fｏr(i1=s［i][j]-1;i1＜=s［i］[ｊ]+１；i１++){fｏr(j１=r[i][ｊ]－１;j1＝r［i］［j]+1;j1++）ﻩ{ｔeｍp１＝１.0;ｆｏr（m=ｓ[i][j］－１;ｍ＜＝s[i][j]＋1;m++）{if（ｍ!=i１)ｔeｍp1*=(t［ｉ][j]-ｔ0［m])/(ｔ0[i1]－ｔ0［m]）；}teｍp２=1．0;ｆor(m=r[i]［j]-1;m＜=ｒ［ｉ][ｊ]＋1;m++）{if(m！＝j1)temp2*=(ｕ［i][j］-u0［m］)/（ｕ0［j1］-u0［m]）;}z[i][j]+=temｐ１*temp2*z0[i1］[j1］;｝ﻩ}ﻩ｝}ﻩ}voiｄqｍｎh()／*曲面拟合子程序*/{inti,j,k，m,ｌ,i1，ｊ1，ｉk;doubleA[N][２１],D[N][21］,B[11]［N］,Bｔ[N][1１],ＢtB[N][N],BtU[N][21］,BｔＢ1［N］[Ｎ]；ｄoublｅＧ［21］[N],Gt[N][2１],ＧtG[N][N］，ＧｔG1[N］[N]；douｂlesigma,p［１1][2１］,teｍp，q;printｆ(选