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

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

2、k=1開始逐漸增大k的值，使用最小二乘法曲面擬合法對z=f(x,y)進行擬合，當時結束計算。擬合基函數(shù)r(x)s(y)選擇為r(x)=xr，s(y)=ys。擬合系數(shù)矩陣c通過連續(xù)兩次解線性方程組求得。， 其中，4、觀察比較計算的值并輸出結果，以觀察逼近的效果。其中。二、全部源程序/ hean.cpp : 定義控制臺應用程序的入口點。/#include stdafx.h#include #include #include void Set_non_JacobiA(double* A,double* x);/求題中非線性方程組對應自變量向量x的雅克比/矩陣void Set_non_B(double

3、* B,double* x,double a,double b);/求非線性方程組Newton迭代法的右/端式:-F(x)void Array_Mult_Array(double* A,double* B,int m,int s,int n,double* C);/矩陣相乘ABCvoid Transpose(double *A,int m,int n,double* AT);/轉置void Doolittle(double *A,int n,int *M);void Solve_LU(double* A,int n,double* b,double* x);/解方程LUxbvoid Solve

4、_lin(double* A,int n,double* B,double* x,int m);/解線性方程組AxBdouble Vector_FanShu(double *V,int n);void Solve_non_Equation(double a,double b,double* x);/求解非線性方程組int Near_Index(double* v,double a,int n);/查找n維向量V中與常數(shù)a最接近的元素的下標double ChaZhi(double a,double b);/求數(shù)表z(t,u)在點（x，y）處的分片二次代數(shù)差值void NiHe(double*U,

5、double*x,double*y,int m,int n,int p,int q,double*C);/對mn數(shù)表U(x,y)進行二元多項式擬合double Pxy(double *C,int p,int q,double x,double y);/求多項式擬合函數(shù)P(x,y)在點（x，y）處的函數(shù)值void main()double non4; /非線性方程組變量向量（t,u,v,w）double f1121;/f(x,y)在擬合節(jié)點處的值double x11,y21;/擬合節(jié)點double *C;/擬合系數(shù)矩陣double det=1e-7;/擬合精度要求double p,d;int i

6、,j,k,t;/設置擬合節(jié)點for(i=0;i11;i+)xi=0.08*i;for(j=0;j21;j+)yj=0.5+0.05*j;/求擬合節(jié)點處的f(x,y)的值for(i=0;i11;i+)for(j=0;j21;j+)Solve_non_Equation(xi,yj,non);fij=ChaZhi(non0,non1);printf(n 數(shù)值分析第三次大作業(yè)n);/打印數(shù)表f(xi,yi)printf(n f(xi,yi) n);for(t=0;t21/3;t+)printf(-n);printf(| f(x,y) |);for(i=3*t;i3*t+3;i+)printf( y=%

7、4.2f |,yi);printf(n);for(j=0;j11;j+)printf(| x=%4.2f |,xj);for(i=3*t;i3*t+3;i+)printf( %+.12e |,fji);printf(n);printf(-n);printf(n);k=0;printf( 不同k對應的精度 );printf(n-);doC=(double*) calloc(k+1)*(k+1),sizeof(double);NiHe(*f,x,y,11,21,k+1,k+1,C);/求在當前k值擬合的節(jié)點處誤差d=0;for(i=0;i11;i+)for(j=0;j=det)free(C);el

8、seprintf(n-);printf(n 故可知k=k=%d%.0e，滿足題設要求,det);break;while(+k11);/打印擬合系數(shù)矩陣cprintf(nnn 當k=%d時擬合函數(shù)p(x,y)的系數(shù)矩陣c：n,k);printf(-n);for(t=0;t(k+1)/3;t+)for(i=0;i(k+1);i+)for(j=3*t;j3*t+3;j+)printf( %+.12e ,Ci*(k+1)+j);printf(n);printf(-n);/打印(x*,y*)處的擬合逼近情況printf(nn觀察以下各點的逼近情況n);printf( | (x*,y*) | f(x*,y

9、*) | p(x*,y*) | |f-p| |n);printf(-n);for(i=1;i=8;i+)for(j=1;j=5;j+)Solve_non_Equation(0.1*i,0.5+0.2*j,non);d=ChaZhi(non0,non1);p=Pxy(C,k+1,k+1,0.1*i,0.5+0.2*j);printf( | (%3.1f,%3.1f)| %+.12e | %+.12e | %f |n,0.1*i,0.5+0.2*j,d,p,abs(d-p);printf(-nnn);/求題中非線性方程組對應自變量向量x的雅克比矩陣void Set_non_JacobiA(doub

10、le* A,double* x)/A為*4階系數(shù)矩陣，x為自變量(t,u,v,w)A0=-0.5*sin(x0);A1=1;A2=1;A3=1;A4=1;A5=0.5*cos(x1);A6=1;A7=1;A8=0.5;A9=1;A10=-sin(x2);A11=1;A12=1;A13=0.5;A14=1;A15=cos(x3);/求非線性方程組Newton迭代法的右端式:-F(x)void Set_non_B(double* B,double* x,double a,double b)/a非線性方程的參數(shù)，對應題中的x/b非線性方程的參數(shù)，對應題中的yB0=-(0.5*cos(x0)+x1+x

11、2+x3-a-2.67);B1=-(x0+0.5*sin(x1)+x2+x3-b-1.07);B2=-(0.5*x0+x1+cos(x2)+x3-a-3.74);B3=-(x0+0.5*x1+x2+sin(x3)-b-0.79);/矩陣相乘ABC，其中A為ms階，B為sn階。void Array_Mult_Array(double* A,double* B,int m,int s,int n,double* C)int i,j,k;for(i=0;im;i+)for(j=0;jn;j+)Ci*n+j=0;for(k=0;ks;k+)Ci*n+j+=Ai*s+k*Bk*n+j;/將mn階矩陣A的

12、轉置為ATvoid Transpose(double *A,int m,int n,double* AT)int i,j;for(i=0;im;i+)for(j=0;jn;j+)ATj*m+i=Ai*n+j;/對n階矩陣A進行選主元的Doolittle分解void Doolittle(double *A,int n,int *M)/將分解得到的L、U仍然存儲在原來的矩陣A中int i,j,k,t;double *s;double Maxs,temp;s=(double*) calloc(n,sizeof(double);for(k=0;kn;k+) for(i=k;in;i+)si=Ai*n+

13、k;for(t=0;tk;t+)si-= Ai*n+t * At*n+k;Maxs=abs(sk); Mk=k;for(i=k+1;in;i+)if(Maxsabs(si) Maxs=abs(si);Mk=i;if(Mk!=k)for(t=0;tn;t+)temp=Ak*n+t;Ak*n+t=AMk*n+t;AMk*n+t=temp;temp=sk;sk=sMk;sMk=temp;if(Maxs(1e-14) sk=1e-14;printf(%.16e方陣奇異n,Maxs);Ak*n+k=sk;for(j=k+1;(jn)&(kn-1);j+)for(t=0;tk;t+)Ak*n+j-=Ak*

14、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;in;i+)xi=bi;for(t=0;t-1;i-)for(t=i+1;tn;t+)xi-=Ai*n+t*xt;xi/=Ai*n+i;/解線性方程組AxBvoid Solve_lin(double* A,int n,double* B,double* x,int m)/B為nm矩陣int* M,i,j;M=(int*) calloc(n,sizeof(int);double *BT,*xT,

15、temp;BT=(double*) calloc(n*m,sizeof(double);xT=(double*) calloc(n*m,sizeof(double);Transpose(B,n,m,BT); Doolittle(A,n,M);/將A三角分解for(i=0;im;i+)for(j=0;jn-1;j+)temp=BTi*n+j;BTi*n+j=BTi*n+Mj;BTi*n+Mj=temp;/將B轉置，使得同一方程組對應的系數(shù)連續(xù)存儲Solve_LU(A,n,BT+i*n,xT+i*n);Transpose(xT,m,n,x);/求n維向量V的無窮范數(shù)double Vector_Fa

16、nShu(double *V,int n)int i;double max=0;for(i=0;in;i+)if(maxabs(Vi)max=abs(Vi);return max;/求解非線性方程組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ù)/設定迭代初始值x0=0.5;x1=1;x2=1;x3=1;doSet_non_B(F,x,a,b);Set_non_JacobiA(*A,x);Solv

17、e_lin(*A,4,F,detx,1);if(Vector_FanShu(detx,4)/Vector_FanShu(x,4)det)return;for(i=0;i4;i+)xi+=detxi;k+;while(kM);printf(Newton法在該初值不收斂n);/查找n維向量V中與常數(shù)a最接近的元素的下標int Near_Index(double* v,double a,int n)int i,Index;double min;min=abs(v0-a);Index=0;for(i=1;iabs(vi-a)min=abs(vi-a);Index=i;return Index;/求數(shù)表

18、z(t,u)在點（x，y）處的分片二次代數(shù)差值double ChaZhi(double a,double b)double z66=-0.5,-0.34,0.14,0.94,2.06,3.5,-0.42,-0.5,-0.26,0.3,1.18,2.38,-0.18,-0.5,-0.5,-0.18,0.46,1.42,0.22,-0.34,-0.58,-0.5,-0.1,0.62,0.78,-0.02,-0.5,-0.66,-0.5,-0.02,1.5,0.46,-0.26,-0.66,-0.74,-0.5;double x6=0,0.2,0.4,0.6,0.8,1;double y6=0,0.

19、4,0.8,1.2,1.6,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é)點的選取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=j-1;t=j+1;t+)if(t!=r)L_r-j+1*=(b-yt)/(yr-y

20、t);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;/對mn數(shù)表U(x,y)進行二元多項式擬合void NiHe(double*U,double*x,double*y,int m,int n,int p,int q,double*C)/U為擬合數(shù)據(jù)點處的函數(shù)值int i,j;double *B,*BT,*BTB,*G,*GT,*GTG,*BTUG,*temp,*tempT,*CT;B=(double*) calloc(m*p,sizeof(double);BT=(double*) callo

21、c(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(p*n,sizeof(double);for(i=0;im;i+)for(j=0;jp;j+)Bi*p+j=pow(xi,j);for(i=0;in;i+)for(j=0;jq;j