北航數(shù)值分析A大作業(yè)3

1、北航數(shù)值分析a大作業(yè)3一、算法設(shè)計(jì)方案1、解非線性方程組將各擬合節(jié)點(diǎn)（xi，yj）分別帶入非線性方程組，求出與相對(duì)應(yīng)的數(shù)組teij，ueij，求解非線性方程組選擇newton迭代法，迭代過(guò)程中需要求解線性方程組，選擇選主元的doolittle分解法。2、二元二次分偏插值對(duì)數(shù)表z(t,u)進(jìn)行分片二次代數(shù)插值，求得對(duì)應(yīng)（tij，uij）處的值，即為 的值。根據(jù)給定的數(shù)表，可將整個(gè)插值區(qū)域分成 16 個(gè)小 的區(qū)域，故先判斷(t ij, u ij ) 所在，的區(qū)域，再作此區(qū)域的插值，計(jì)算 z ij，相應(yīng)的lagrange形式的插值多項(xiàng)式為：其中 (k=m-1, m, m+1) (r=n-1, n,

2、n+1)3、曲面擬合從k=1開(kāi)始逐漸增大k的值，使用最小二乘法曲面擬合法對(duì)z=f(x,y)進(jìn)行擬合，當(dāng)時(shí)結(jié)束計(jì)算。擬合基函數(shù)r(x)s(y)選擇為r(x)=xr，s(y)=ys。擬合系數(shù)矩陣c通過(guò)連續(xù)兩次解線性方程組求得。， 其中，4、觀察比較計(jì)算的值并輸出結(jié)果，以觀察逼近的效果。其中。二、全部源程序2e |",fji);printf("n");printf("-n");printf("n");k=0;printf(" 不同k對(duì)應(yīng)的精度 ");printf("n-");doc=(dou

3、ble*) calloc(k+1)*(k+1),sizeof(double);nihe(*f,x,y,11,21,k+1,k+1,c);2e",k,d);if(d>=det)free(c);elseprintf("n-");printf("n 故可知k=k=%d<%.0e，滿(mǎn)足題設(shè)要求",det);break;while(+k<11);2e ",ci*(k+1)+j);printf("n");printf("-n");2e | %+.12e | %f |n",*i,+

4、*j,d,p,abs(d-p);printf("-nnn");6e方陣奇異n",maxs);ak*n+k=sk;for(j=k+1;(j<n)&&(k<n-1);j+)for(t=0;t<k;t+)ak*n+j-=ak*n+t*at*n+j;aj*n+k=sj/ak*n+k;/解方程luxbvoid solve_lu(double* a,int n,double* b,double* x)int i,t;for(i=0;i<n;i+)xi=bi;for(t=0;t<i;t+)xi-=ai*n+t*xt;for(i=n-

5、1;i>-1;i-)for(t=i+1;t<n;t+)xi-=ai*n+t*xt;xi/=ai*n+i;/解線性方程組axbvoid solve_lin(double* a,int n,double* b,double* x,int m)/b為n×m矩陣int* m,i,j;m=(int*) calloc(n,sizeof(int);double *bt,*xt,temp;bt=(double*) calloc(n*m,sizeof(double);xt=(double*) calloc(n*m,sizeof(double);transpose(b,n,m,bt); do

6、olittle(a,n,m);/將a三角分解for(i=0;i<m;i+)for(j=0;j<n-1;j+)temp=bti*n+j;bti*n+j=bti*n+mj;bti*n+mj=temp;/將b轉(zhuǎn)置，使得同一方程組對(duì)應(yīng)的系數(shù)連續(xù)存儲(chǔ)solve_lu(a,n,bt+i*n,xt+i*n);transpose(xt,m,n,x);/求n維向量v的無(wú)窮范數(shù)double vector_fanshu(double *v,int n)int i;double max=0;for(i=0;i<n;i+)if(max<abs(vi)max=abs(vi);return max;

7、/求解非線性方程組void solve_non_equation(double a,double b,double* x)double a44,f4,detx4;double det=1e-12;/求解精度要求int i,k=0;int m=5000;/最大迭代次數(shù)/設(shè)定迭代初始值x0=;x1=1;x2=1;x3=1;doset_non_b(f,x,a,b);set_non_jacobia(*a,x);solve_lin(*a,4,f,detx,1);if(vector_fanshu(detx,4)/vector_fanshu(x,4)<det)return;for(i=0;i<4

8、;i+)xi+=detxi;k+;while(k<m);printf("newton法在該初值不收斂n");/查找n維向量v中與常數(shù)a最接近的元素的下標(biāo)int near_index(double* v,double a,int n)int i,index;double min;min=abs(v0-a);index=0;for(i=1;i<n;i+)if(min>abs(vi-a)min=abs(vi-a);index=i;return index;/求數(shù)表z(t,u)在點(diǎn)（x，y）處的分片二次代數(shù)差值double chazhi(double a,doub

9、le b)double z66=,;double x6=0,1;double y6=0,2;double l3,l_3,z;int i,j,k,r,t;i=near_index(x,a,6);j=near_index(y,b,6);/插值區(qū)域邊界處插值節(jié)點(diǎn)的選取if(i=0)i=1;if(i=5)i=4;if(j=0)j=1;if(j=5)j=4;for(k=i-1;k<=i+1;k+)lk-i+1=1;for(t=i-1;t<=i+1;t+)if(t!=k)lk-i+1*=(a-xt)/(xk-xt);for(r=j-1;r<=j+1;r+)l_r-j+1=1;for(t=

10、j-1;t<=j+1;t+)if(t!=r)l_r-j+1*=(b-yt)/(yr-yt);z=0;for(k=i-1;k<=i+1;k+)for(r=j-1;r<=j+1;r+)z+=(lk-i+1*l_r-j+1*zkr);return z;/對(duì)m×n數(shù)表u(x,y)進(jìn)行二元多項(xiàng)式擬合void nihe(double*u,double*x,double*y,int m,int n,int p,int q,double*c)/u為擬合數(shù)據(jù)點(diǎn)處的函數(shù)值int i,j;double *b,*bt,*btb,*g,*gt,*gtg,*btug,*temp,*tempt,

11、*ct;b=(double*) calloc(m*p,sizeof(double);bt=(double*) calloc(p*m,sizeof(double);btb=(double*) calloc(p*p,sizeof(double);g=(double*) calloc(n*q,sizeof(double);gt=(double*) calloc(q*n,sizeof(double);gtg=(double*) calloc(q*q,sizeof(double);btug=(double*) calloc(p*q,sizeof(double);temp=(double*) calloc

12、(p*n,sizeof(double);for(i=0;i<m;i+)for(j=0;j<p;j+)bi*p+j=pow(xi,j);for(i=0;i<n;i+)for(j=0;j<q;j+)gi*q+j=pow(yi,j);transpose(b,m,p,bt);transpose(g,n,q,gt);array_mult_array(bt,b,p,m,p,btb);array_mult_array(gt,g,q,n,q,gtg);array_mult_array(bt,u,p,m,n,temp);array_mult_array(temp,g,p,n,q,btug

13、);free(b);free(bt);free(g);free(gt);free(temp);temp=(double*) calloc(p*q,sizeof(double);solve_lin(btb,p,btug,temp,q);free(btb);free(btug);tempt=(double*) calloc(q*p,sizeof(double);ct=(double*) calloc(q*p,sizeof(double);transpose(temp,p,q,tempt);solve_lin(gtg,q,tempt,ct,p);transpose(ct,q,p,c);free(temp);free(ct);free(