北航研究生數(shù)值分析上機(jī)作業(yè) 三 (報(bào)告+所有程序大全)

1、 數(shù)值分析上機(jī)作業(yè)3求解非線性方程組以及二元函數(shù)的插值擬合1.算法設(shè)計(jì)對(duì)于全部的插值節(jié)點(diǎn)(x,y),i二0,1,.,10,j二0,1,.,20,帶入非線性方程組中，用Newtonij迭代法解非線性方程組，得到(t,u),i二0,1,.,10,j二0,1,.,20。對(duì)(t,u)，在二維數(shù)表ijij中進(jìn)行插值，采用分片雙二次插值法。插值過(guò)程中，先選擇分片區(qū)域的中心節(jié)點(diǎn)，在數(shù)表中的列記為tt(0:5)，行記為uu(0:5)，中心節(jié)點(diǎn)記為(a,b)，生成向量t_temp(0:2)，t_temp(0)=(t一tt(a)(t一tt(a+1)/(tt(a-1)一tt(a)(tt(a-1)一tt(a+1)，i

2、it_temp(1)=(t一tt(a一1)(t一tt(a+1)/(tt(a)一tt(a一1)(tt(a)一tt(a+1)，iit_temp(2)=(t一tt(a一1)(t一tt(a)/(tt(a+1)一tt(a一1)(tt(a+1)一tt(a)，ii同理，生成向量u_temp(0:2)，u_temp(0)=(u一uu(a)(u一uu(a+1)/(uu(a一1)一uu(a)(uu(a一1)一uu(a+1)jju_temp(1)=(u一uu(a一1)(u一uu(a+1)/(uu(a)一uu(a一1)(uu(a)一uu(a+1)jju_temp(2)=(u一uu(a一1)(u一uu(a)/(uu(a

3、+1)一uu(a一1)(uu(a+1)一uu(a)jj記數(shù)表中以分片區(qū)域中心節(jié)點(diǎn)為中心的3X3的矩陣為T,對(duì)于(t,u)插值結(jié)果為ij(t_tempT)(T)(u_temp)。在擬合p(X,y)二丈cxrys時(shí)，需要計(jì)算C二(BtB)-1BtUG(GtG)j，令rsr,s=0W=BtB,M=GtG,計(jì)算W,M時(shí)，根據(jù)對(duì)稱性只需要計(jì)算對(duì)角線元素和對(duì)角線以上元素即可，節(jié)省運(yùn)算時(shí)間。于是WCM=BtUG=E，用選主元的LU分解法求解WF=E，再計(jì)算MCt=Ft，這里C,F只需要按行取元素進(jìn)行運(yùn)算即可，故不需要進(jìn)行轉(zhuǎn)置運(yùn)算。k從1到9依次增加，計(jì)算b的值，當(dāng)b10-7時(shí)，得到達(dá)到精度的最小的k。2.打

4、印輸出結(jié)果f(x,y):4.46504018480e-0013.24683262927e-0012.10159686683e-0013.40189556266e-003-8.87358136384e-0026.38015226510e-0015.06611755146e-0013.82176369277e-0011.54780200285e-0015.19926834902e-0028.40081395765e-0016.99764165673e-0015.66061442351e-0013.19242138041e-0012.06376192387e-0011.05151509180e+00

5、09.02927430830e-0017.60580266859e-0014.95519756001e-0013.73134042774e-0011.27124675148e+0001.11500201815e+0009.64607727215e-0016.82447678179e-0015.51085208597e-0011.49832105248e+0001.33499863207e+0001+0008.78960023174e-0017.39145108703e-0011.73189274038e+0001.56203457721e+0001.397216918

6、21e+0001.08408753268e+0009.36322772315e-0011.97122178640e+0001.79532959950e+0001.62406711323e+0001.29695464975e+0001+000k=1,sigma=3.22090897363e+000k=2,sigma=4.65996003320e-003k=3,sigma=1.72117537926e-004k=4,sigma=3.30953430190e-006k=5,sigma=2.54137773513e-008C_rs:2.02123036395e+000-3.6

7、6842591519e+0007.09246688452e-001-4-0016.74322022261e-0023.19192646165e+000-7.41209812962e-001-2.69690653026e+000-4.84601977266e-0016.05908514261e-0022.56706343343e-0011.58096413206e+000-4.65701259544e-0011.00853116187e-001-2.07687560816e-002-2.68608304872e-001-7.33963449843e-0011.08429

8、601112e+0003.07285137544e-001-4.68489486618e-0022.16521800059e-001-1.73026852965e-001-8.41324310602e-002-1.47683427939e-0012.77711894906e-0021.03043608316e-0012.64863491154e-0014.39171608118e-0016.24715198145e-0018.20347369475e-0011.02502405502e+0001.23780100674e+0001.45783058271e+0008.48607444215e-

9、0011.63095275395e+000-7.89186706497e-002-8-0012.55736987891e-0013.44966421960e-002-5.54328606191e-0021.40518220408e-001-1.30388672239e-0016.95959872154e-003-3.30022941240e-003f(x*,y*):0.1947200.4059790.6347770.8789601.1366111.4060421.6857841.974575p(x*,y*):0.1947300.4059900.6347870.8789

10、701.1366201.4060511.6857911.974581-0.183037-0.0225160.1588010.3586510.5749800.8059411.0498811.305335-0.183042-0.0225210.1587960.3586460.5749760.8059371.0498781.305332-0.445498-0.338221-0.207366-0.0552530.1159920.3044290.5082940.725992-0.445500-0.338224-0.207369-0.0552550.1159890.3044260.5082910.7259

11、89-0.597567-0.544438-0.465358-0.362680-0.238568-0.0950160.0661490.243254-0.597559-0.544430-0.465350-0.362671-0.238560-0.0950090.0661560.243261-0.646460-0.647361-0.620271-0.567565-0.491434-0.393902-0.276834-0.141950-0.646446-0.647348-0.620257-0.567551-0.491421-0.393890-0.276822-0.141939對(duì)f(x,y)插值的數(shù)值情況

12、:2*0.600.60.80.70.60.60.40.30.2f(x,y)p(x,y)的數(shù)值情況(z坐標(biāo)為10-7)iiii0.72*0.80.9*01.11.21.3OJ0.60.50.41.50.30.23.源程序1)主函數(shù)#include#include#definelength_x11/x向量長(zhǎng)度#definelength_y21/y向量長(zhǎng)度#definex_length8/x*向量長(zhǎng)度#definey_length5/y*向量長(zhǎng)度externdoubleinterpolation(doublett,doubleuu,double*table);/輸入一點(diǎn)，輸出一點(diǎn)插值externdo

13、ublerow_multi_sum(double*B,inti,intj,intlength);/矩陣列相乘intmain()doublett=0;doubleuu=0;doubletable66=0;doubleff=0;doublexxlength_x=0;doubleyylength_y=0;doublesum=0.0;doubleUlength_xlength_y=0;doubleBlength_x10=0;doubleGlength_y10=0;doubleW1010=0;doubleT1010=0;doubleM1010=0;doubleC1010=0;doubleF10lengt

14、h_y=0;doubleE1010=0;doubleXX1010=0;/doublefflength_xlength_y=0;/記錄f陣doublepplength_xlength_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;table00=-0.5;table01=-0.34;table02=0.14;table03=0.94;table04=2.06;table05=3.5;table10=-0.42;table11=-0.5;table12=-0.2

15、6;table13=0.3;table14=1.18;table15=2.38;table20=-0.18;table21=-0.5;table22=-0.5;table23=-0.18;table24=0.46;table25=1.42;table30=0.22;table31=-0.34;table32=-0.58;table33=-0.5;table34=-0.1;table35=0.62;table40=0.78;table41=-0.02;table42=-0.5;table43=-0.66;table44=-0.5;table45=-0.02;table50=1.5;table51

16、=0.46;table52=-0.26;table53=-0.66;table54=-0.74;table55=-0.5;for(ii=0;iilength_x;ii+)xxii=0.08*ii;for(ii=0;iilength_y;ii+)yyii=0.5+0.05*ii;for(ii=0;iilength_x;ii+)/生成U陣for(jj=0;jjlength_y;jj+)ite_equations(&tt,&uu,xxii,yyjj);/牛頓迭代解方程，輸入一點(diǎn)(x,y)輸出一點(diǎn)(t,u)Uiijj=(double)interpolation(tt,uu,(double*)table

17、);/輸入一點(diǎn)，輸出一點(diǎn)插值fprintf(fd,nf(x,y):n);for(ii=0;iix_length;ii+)for(jj=0;jj=y_length;jj+)fprintf(fd,%.11e,Uiijj);fprintf(fd,n);for(k_max=1;k_max=6;k_max+)for(kk=0;kk(k_max+1);kk+)/生成B陣for(ii=0;iilength_x;ii+)Biikk=pow(xxii,kk);for(kk=0;kk(k_max+1);kk+)for(ii=0;iilength_y;ii+)/生成G陣Giikk=pow(yyii,kk);for(

18、ii=0;ii=k_max;ii+)/BT*Bfor(jj=0;jj=k_max;jj+)Wiijj=row_multi_sum(double*)B,ii,jj,length_x);for(ii=0;ii=k_max;ii+)/GT*Gfor(jj=0;jj=k_max;jj+)Miijj=row_multi_sum(double*)G,ii,jj,length_y);for(ii=0;ii=k_max;ii+)/F=B*Ufor(jj=0;jjlength_y;jj+)Fiijj=0;for(kk=0;kklength_x;kk+)Fiijj=Fiijj+Bkkii*Ukkjj;for(ii

19、=0;ii=k_max;ii+)/E=F*Gfor(jj=0;jj=k_max;jj+)Eiijj=0;for(kk=0;kklength_y;kk+)Eiijj=Eiijj+Fiikk*Gkkjj;Matrix_LU(double*)W,(double*)E,(double*)M,(double*)C,k_max,0);/(double*)XX,for(ii=0;iilength_x;ii+)for(jj=0;jjlength_y;jj+)sum=0;for(r=0;r=k_max;r+)for(s=0;s=k_max;s+)sum=sum+Crs*pow(xxii,r)*pow(yyjj,

20、s);ppiijj=sum;sum=0;for(ii=0;iilength_x;ii+)for(jj=0;jjlength_y;jj+)sum=sum+(ppiijj-Uiijj)*(ppiijj-Uiijj);fprintf(fd,nk=%d,sigma=%.11en,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,Crs);fprintf(fd,n);break;for(ii=0;iix_length;ii+)xxii=0.1*(ii+1);f

21、or(ii=0;iiy_length;ii+)yyii=0.5+0.2*(ii+1);for(ii=0;iix_length;ii+)/生成f*陣for(jj=0;jjy_length;jj+)ite_equations(&tt,&uu,xxii,yyjj);/牛頓迭代解方程，輸入一點(diǎn)(x,y)輸出一點(diǎn)(t,u)Uiijj=(double)interpolation(tt,uu,(double*)table);/輸入一點(diǎn)，輸出一點(diǎn)插值fprintf(fd,nf(x*,y*):n);for(ii=0;iix_length;ii+)for(jj=0;jjy_length;jj+)fprintf(f

22、d,%-10f,Uiijj);fprintf(fd,n);for(ii=0;iix_length;ii+)/生成p*陣for(jj=0;jjy_length;jj+)sum=0;for(r=0;r=k_max;r+)for(s=0;s=k_max;s+)sum=sum+Crs*pow(xxii,r)*pow(yyjj,s);ppiijj=sum;fprintf(fd,np(x*,y*):n);for(ii=0;iix_length;ii+)for(jj=0;jjy_length;jj+)fprintf(fd,%-10f,ppiijj);fprintf(fd,n);fclose(fd);2)矩陣

23、2列相乘#include#include#definesize10#defineTTT(a,b)*(double*)TTT+(a)*(size)+b)/LU的宏定義，以方便操作二維數(shù)組doublerow_multi_sum(double*B,inti,intj,intlength)/矩陣列相乘doublesum=0;double*TTT;intkk=0;TTT=B;for(kk=0;kklength;kk+)sum=sum+(TTT(kk,i)*(TTT(kk,j);returnsum;3)9點(diǎn)分片雙二次插值#include#include#definesize6#definetable(a,

24、b)*(double*)table+(a)*(size)+b)/LU的宏定義，以方便操作二維數(shù)組doubleinterpolation(doublett,doubleuu,double*table)/輸入一點(diǎn)，輸出一點(diǎn)插值intii=0;intjj=0;intsection_t=0;intsection_u=0;doubleanswer=0;doubletemp=0;doubletemp_u3=0;doubletemp_t3=0;doublemy_t6=0;doublemy_u6=0;for(ii=1;ii=0)&(uu=0.6)&(uu=1.0)&(uu=1.4)&(uu=0)&(tt=0.

25、3)&(tt=0.5)&(tt=0.7)&(tt=1.0)section_t=4;else;temp_u0=(uu-my_usection_u)*(uu-my_usection_u+1)/(my_usection_u-1-my_usection_u)*(my_usection_u-1-my_usection_u+1);temp_u1=(uu-my_usection_u-1)*(uu-my_usection_u+1)/(my_usection_u-my_usection_u-1)*(my_usection_u-my_usection_u+1);temp_u2=(uu-my_usection_u-

26、1)*(uu-my_usection_u)/(my_usection_u+1-my_usection_u-1)*(my_usection_u+1-my_usection_u);temp_t0=(tt-my_tsection_t)*(tt-my_tsection_t+1)/(my_tsection_t-1-my_tsection_t)*(my_tsection_t-1-my_tsection_t+1);temp_t1=(tt-my_tsection_t-1)*(tt-my_tsection_t+1)/(my_tsection_t-my_tsection_t-1)*(my_tsection_t-m

27、y_tsection_t+1);temp_t2=(tt-my_tsection_t-1)*(tt-my_tsection_t)/(my_tsection_t+1-my_tsection_t-1)*(my_tsection_t+1-my_tsection_t);answer=0;for(ii=0;ii3;ii+)answer=answer+temp_tii*(temp_u0*table(section_t-1+ii),section_u-1)+temp_u1*table(section_t-1+ii),section_u)+temp_u2*table(section_t-1+ii),sectio

28、n_u+1);returnanswer;4)Newton迭代求解非線性方程組#include#includevoidite_equations(double*tt,double*uu,doublexx,doubleyy)/牛頓迭代解非線性方程組doublet=1;doubleu=1;doublev=1;doublew=1;doubled4=0;/增量doublex=xx;doubley=yy;doubleF_Jacob44=0;doubleF4=0;doubleL99=0;doubleU99=0;doubletemp1=0;doubletemp2=0;intii=0;intjj=0;intkk

29、=0;for(ii=0;ii=1e2;ii+)F_Jacob00=-0.5*sin(t);BY1102130劉晚春電子信息工程學(xué)院BY1102130劉晚春電子信息工程學(xué)院 F_Jacob01=1;F_Jacob02=1;F_Jacob03=1;F_Jacob10=1;F_Jacob11=0.5*cos(u);F_Jacob12=1;F_Jacob13=1;F_Jacob20=0.5;F_Jacob21=1;F_Jacob22=-sin(v);F_Jacob23=1;F_Jacob30=1;F_Jacob31=0.5;F_Jacob32=1;F_Jacob33=cos(w);F0=0.5*cos

30、(t)+u+v+w-x-2.67;F1=t+0.5*sin(u)+v+w-y-1.07;F2=0.5*t+u+cos(v)+w-x-3.74;F3=t+0.5*u+v+sin(w)-y-0.79;LU(double*)L,(double*)U,(double*)F_Jacob,(double*)F);huidai(double*)L,(double*)U,(double*)F,(double*)d);t=t+d0;u=u+d1;v=v+d2;w=w+d3;temp1=sqrt(t*t+u*u+v*v+w*w);temp2=sqrt(d0*d0+d1*d1+d2*d2+d3*d3);if(tem

31、p2/temp1)1e-12)*tt=t;*uu=u;break;5)解線性方程組用到LU分解#include#include#definesize4/#definehalfBW3#defineL(a,b)*(double*)L+(a)*(size)+b)/LU的宏定義，以方便操作二維數(shù)組#defineU(a,b)*(double*)U+(a)*(size)+b)#defineA(a,b)*(double*)A+(a)*(size)+b)/A的宏定義，以方便操作二維數(shù)組/#defineB(a,b)*(double*)B+(a)*(size)+b)/B的宏定義，以方便操作二維數(shù)組/#define

32、X(a,b)*(double*)X+(a)*(size)+b)/A的宏定義，以方便操作二維數(shù)組#defineb(ii)*(double*)b+ii)voidLU(double*LL,double*UU,double*AA,double*bb)/inta/double*L;double*U;double*A;doubleA_tempsizesize=0;double*b;doubleysize=0;doublexsize=0;intkk,ii,jj,pp,qq;doubles,s_temp=0;doublesum_answer=0;doubletemp=0;ints_count=0;intn=s

33、ize;/L=LL;U=UU;A=AA;b=bb;/kk=0;ii=0;jj=0;pp=0;qq=0;/for(ii=0;iin;ii+)/對(duì)A_temp賦值for(jj=0;jjn;jj+)A_tempiijj=A(ii,jj);for(ii=0;iin;ii+)L(ii,ii)=1.0;for(kk=0;kkn;kk+)/總循環(huán)求出LUsum_answer=0;s_count=kk;for(pp=kk;ppn;pp+)/doolitte選主元if(kk=0)sum_answer=0;elsefor(qq=0;qqs_temp)s_temp=s;s_count=pp;elseif(s(-1.

34、0)*s_temp)s_temp=(-1.0)*s;s_count=pp;else;for(pp=0;ppn;pp+)/換位A_temptemp=A_temps_countpp;A_temps_countpp=A_tempkkpp;A_tempkkpp=temp;for(pp=0;ppkk;pp+)/換位Ltemp=L(s_count,pp);L(s_count,pp)=L(kk,pp);L(kk,pp)=temp;temp=b(s_count);b(s_count)=b(kk);b(kk)=temp;for(ii=kk;iin;ii+)/先求一行Uif(kk=0)U(kk,ii)=A_tem

35、pkkii;elsesum_answer=0;for(qq=0;qqkk;qq+)sum_answer=sum_answer+L(kk,qq)*U(qq,ii);U(kk,ii)=A_tempkkii-sum_answer;for(ii=kk+1;iin;ii+)/求出一列Lsum_answer=0;for(qq=0;qqkk;qq+)sum_answer=sum_answer+L(ii,qq)*U(qq,kk);L(ii,kk)=(A_tempiikk-sum_answer)/U(kk,kk);6)回代求解方程組#include#include#definesize4#defineL(a,b

36、)*(double*)L+(a)*(size)+b)/LU的宏定義，以方便操作二維數(shù)組#defineU(a,b)*(double*)U+(a)*(size)+b)#definex(a)*(double*)x+a)/LU的宏定義，以方便操作二維數(shù)組#definey(a)*(double*)y+a)#defineb(a)*(double*)b+a)/#definexx(a)xxa/A*uu=yy=L*U*uu=yy=L*xx=yy&U*uu=xxvoidhuidai(double*LL,double*UU,double*bb,double*xx)/回代法求xx,LUx=bdouble*L;doub

37、le*U;double*b;double*x;doubleysize=0;doublesum_answer=0;intii=0;intjj=0;L=LL;U=UU;b=bb;x=xx;/迭代Ly=by0=-b(0);for(ii=1;iisize;ii+)sum_answer=0;for(jj=1;jj=0;ii-)sum_answer=0;for(jj=ii+1;jjsize;jj+)sum_answer=sum_answer+U(ii,jj)*x(jj);x(ii)=(yii-sum_answer)/U(ii,ii);7)求解C=(B*B)-1*(B*U*G)*(G*G)#include#

38、include#definesize10/#definehalfBW3#defineL(a,b)Lab/*(double*)L+(a)*(size)+b)/LU的宏定義，以方便操作二維數(shù)組#defineU(a,b)Uab/*(double*)U+(a)*(size)+b)#defineA(a,b)*(double*)A+(a)*(size)+b)/A的宏定義，以方便操作二維數(shù)組#defineB(a,b)*(double*)B+(a)*(size)+b)/B的宏定義，以方便操作二維數(shù)組#defineM(a,b)*(double*)M+(a)*(size)+b)/B的宏定義，以方便操作二維數(shù)組#d

39、efineC(a,b)*(double*)C+(a)*(size)+b)/B的宏定義，以方便操作二維數(shù)組/#defineX(a,b)*(double*)X+(a)*(size)+b)/A的宏定義，以方便操作二維數(shù)組#defineX(a,b)Xab/#defineb(ii)*(double*)b+ii)voidMatrix_LU(double*AA,double*BB,double*MM,double*CC,intk_max,intmode)double*A;double*B;double*M;double*C;doubleXsizesize=0;intn=k_max+1;doubleA_tem

40、psizesize=0;doubleLsizesize=0;doubleUsizesize=0;doubleYsizesize=0;double*b;intkk,ii,jj,pp,qq;doubles,s_temp=0;doublesum_answer=0;doubletemp=0;ints_count=0;/L=LL;/U=UU;A=AA;B=BB;M=MM;C=CC;/X=XX;/kk=0;ii=0;jj=0;pp=0;qq=0;/for(ii=0;iin;ii+)/對(duì)A_temp賦值for(jj=0;jjn;jj+)A_tempiijj=A(ii,jj);for(ii=0;iin;ii+

41、)L(ii,ii)=1.0;for(kk=0;kkn;kk+)/總循環(huán)求出LUsum_answer=0;s_count=kk;for(pp=kk;ppn;pp+)/doolitte選主元if(kk=0)sum_answer=0;elsefor(qq=0;qqs_temp)s_temp=s;s_count=pp;elseif(s(-1.0)*s_temp)s_temp=(-1.0)*s;s_count=pp;else;for(pp=0;ppn;pp+)/換位A_temptemp=A_temps_countpp;A_temps_countpp=A_tempkkpp;A_tempkkpp=temp;

42、for(pp=0;ppkk;pp+)/換位Ltemp=L(s_count,pp);L(s_count,pp)=L(kk,pp);L(kk,pp)=temp;for(pp=0;ppn;pp+)/換位B/temp=B(s_count,pp);B(s_count,pp)=B(kk,pp);B(kk,pp)=temp;for(ii=kk;iin;ii+)/先求一行Uif(kk=0)U(kk,ii)=A_tempkkii;elsesum_answer=0;for(qq=0;qqkk;qq+)sum_answer=sum_answer+L(kk,qq)*U(qq,ii);U(kk,ii)=A_tempkkii-sum_answer;for(ii=kk+1;iin;ii+)/求出一列Lsum_answer=0;for(qq=0;qqkk;qq+)sum_answer=sum_answer+L(ii,qq)*