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

一、算法設(shè)計(jì)方案1、解非線性方程組將各擬合節(jié)點(diǎn)（xi，yj）分別帶入非線性方程組，求出與相對應(yīng)的數(shù)組te[i][j]，ue[i][j]，求解非線性方程組選擇Newton迭代法，迭代過程中需要求解線性方程組，選擇選主元的Doolittle分解法。2、二元二次分偏插值對數(shù)表z(t,u)進(jìn)行分片二次代數(shù)插值，求得對應(yīng)（tij，uij）處的值，即為的值。根據(jù)給定的數(shù)表可將整個插值區(qū)域分成16個小的區(qū)域，故先判斷tij,uij所在，的區(qū)域，再作此區(qū)域的插值，計(jì)算zij，相應(yīng)的Lagrange形式的插值多項(xiàng)式為：其中(k=m-1,m,m+1)(r=n-1,n,n+1)3、曲面擬合從k=1開始逐漸增大k的值，使用最小二乘法曲面擬合法對z=f(x,y)進(jìn)行擬合，當(dāng)時結(jié)束計(jì)算。擬合基函數(shù)φr(x)ψs(y)選擇為φr(x)=xr，ψs(y)=ys。擬合系數(shù)矩陣c通過連續(xù)兩次解線性方程組求得。，其中，4、觀察比較計(jì)算的值并輸出結(jié)果，以觀察逼近的效果。其中。二、全部源程序//hean.cpp:定義控制臺應(yīng)用程序的入口點(diǎn)。//#include"stdafx.h"#include<stdio.h>#include<stdlib.h>#include<math.h>voidSet_non_JacobiA(double*A,double*x);//求題中非線性方程組對應(yīng)自變量向量x的雅克比//矩陣voidSet_non_B(double*B,double*x,doublea,doubleb);//求非線性方程組Newton迭代法的右//端式:-F(x)voidArray_Mult_Array(double*A,double*B,intm,ints,intn,double*C);//矩陣相乘AB＝CvoidTranspose(double*A,intm,intn,double*AT);//轉(zhuǎn)置voidDoolittle(double*A,intn,int*M);voidSolve_LU(double*A,intn,double*b,double*x);//解方程LUx＝bvoidSolve_lin(double*A,intn,double*B,double*x,intm);//解線性方程組Ax＝BdoubleVector_FanShu(double*V,intn);voidSolve_non_Equation(doublea,doubleb,double*x);//求解非線性方程組intNear_Index(double*v,doublea,intn);//查找n維向量V中與常數(shù)a最接近的元素的下標(biāo)doubleChaZhi(doublea,doubleb);//求數(shù)表z(t,u)在點(diǎn)（x，y）處的分片二次代數(shù)差值voidNiHe(double*U,double*x,double*y,intm,intn,intp,intq,double*C);//對m×n數(shù)表U(x,y)進(jìn)行二元多項(xiàng)式擬合doublePxy(double*C,intp,intq,doublex,doubley);//求多項(xiàng)式擬合函數(shù)P(x,y)在點(diǎn)（x，y）處的函數(shù)值voidmain(){ doublenon[4]; //非線性方程組變量向量（t,u,v,w） doublef[11][21]; //f(x,y)在擬合節(jié)點(diǎn)處的值 doublex[11],y[21]; //擬合節(jié)點(diǎn) double*C; //擬合系數(shù)矩陣 doubledet=1e-7; //擬合精度要求 doublep,d; inti,j,k,t; //設(shè)置擬合節(jié)點(diǎn) for(i=0;i<11;i++) x[i]=0.08*i; for(j=0;j<21;j++) y[j]=0.5+0.05*j; //求擬合節(jié)點(diǎn)處的f(x,y)的值 for(i=0;i<11;i++) { for(j=0;j<21;j++) { Solve_non_Equation(x[i],y[j],non); f[i][j]=ChaZhi(non[0],non[1]); } } printf("\n數(shù)值分析第三次大作業(yè)\n"); //打印數(shù)表f(xi,yi) printf("\nf(xi,yi)\n"); for(t=0;t<21/3;t++) { printf("-------------------------------------------------------------------------------------\n"); printf("|f(x,y)|"); for(i=3*t;i<3*t+3;i++) { printf("y=%4.2f|",y[i]); } printf("\n"); for(j=0;j<11;j++) { printf("|x=%4.2f|",x[j]); for(i=3*t;i<3*t+3;i++) printf("%+.12e|",f[j][i]); printf("\n"); } printf("-------------------------------------------------------------------------------------\n"); printf("\n"); } k=0; printf("不同k對應(yīng)的精度"); printf("\n----------------------------------------------"); do{ C=(double*)calloc((k+1)*(k+1),sizeof(double)); NiHe(*f,x,y,11,21,k+1,k+1,C); //求在當(dāng)前k值擬合的節(jié)點(diǎn)處誤差 d=0; for(i=0;i<11;i++) { for(j=0;j<21;j++) { p=Pxy(C,k+1,k+1,x[i],y[j]); d+=((f[i][j]-p)*(f[i][j]-p)); } } printf("\nk=%d對應(yīng)的σ=%.12e",k,d); if(d>=det) free(C); else { printf("\n----------------------------------------------"); printf("\n故可知k=k=%d<%.0e，滿足題設(shè)要求",det); break; } }while(++k<11); //打印擬合系數(shù)矩陣c printf("\n\n\n當(dāng)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;j<3*t+3;j++) printf("%+.12e",C[i*(k+1)+j]); printf("\n"); } printf("----------------------------------------------------------------------------\n"); } //打印(x*,y*)處的擬合逼近情況 printf("\n\n觀察以下各點(diǎn)的逼近情況\n"); printf("|(x*,y*)|f(x*,y*)|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(non[0],non[1]); 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("-----------------------------------------------------------------------\n\n\n");}//求題中非線性方程組對應(yīng)自變量向量x的雅克比矩陣voidSet_non_JacobiA(double*A,double*x)//A為*4階系數(shù)矩陣，x為自變量(t,u,v,w){ A[0]=-0.5*sin(x[0]); A[1]=1; A[2]=1; A[3]=1; A[4]=1; A[5]=0.5*cos(x[1]); A[6]=1; A[7]=1; A[8]=0.5; A[9]=1; A[10]=-sin(x[2]); A[11]=1; A[12]=1; A[13]=0.5; A[14]=1; A[15]=cos(x[3]);}//求非線性方程組Newton迭代法的右端式:-F(x)voidSet_non_B(double*B,double*x,doublea,doubleb)//a非線性方程的參數(shù)，對應(yīng)題中的x //b非線性方程的參數(shù)，對應(yīng)題中的y{ B[0]=-(0.5*cos(x[0])+x[1]+x[2]+x[3]-a-2.67); B[1]=-(x[0]+0.5*sin(x[1])+x[2]+x[3]-b-1.07); B[2]=-(0.5*x[0]+x[1]+cos(x[2])+x[3]-a-3.74); B[3]=-(x[0]+0.5*x[1]+x[2]+sin(x[3])-b-0.79);}//矩陣相乘AB＝C，其中A為m×s階，B為s×n階。voidArray_Mult_Array(double*A,double*B,intm,ints,intn,double*C){ inti,j,k; for(i=0;i<m;i++) for(j=0;j<n;j++) { C[i*n+j]=0; for(k=0;k<s;k++) C[i*n+j]+=A[i*s+k]*B[k*n+j]; }}//將m×n階矩陣A的轉(zhuǎn)置為ATvoidTranspose(double*A,intm,intn,double*AT){ inti,j; for(i=0;i<m;i++) for(j=0;j<n;j++) AT[j*m+i]=A[i*n+j];}//對n階矩陣A進(jìn)行選主元的Doolittle分解voidDoolittle(double*A,intn,int*M)//將分解得到的L、U仍然存儲在原來的矩陣A中{ inti,j,k,t; double*s; doubleMaxs,temp; s=(double*)calloc(n,sizeof(double)); for(k=0;k<n;k++) { for(i=k;i<n;i++) { s[i]=A[i*n+k]; for(t=0;t<k;t++) s[i]-=A[i*n+t]*A[t*n+k]; } Maxs=abs(s[k]);M[k]=k; for(i=k+1;i<n;i++) { if(Maxs<abs(s[i])) { Maxs=abs(s[i]); M[k]=i; } } if(M[k]!=k) { for(t=0;t<n;t++) { temp=A[k*n+t]; A[k*n+t]=A[M[k]*n+t]; A[M[k]*n+t]=temp; } temp=s[k]; s[k]=s[M[k]]; s[M[k]]=temp; } if(Maxs<(1e-14)) { s[k]=1e-14; printf("%.16e方陣奇異\n",Maxs); } A[k*n+k]=s[k]; for(j=k+1;(j<n)&&(k<n-1);j++) { for(t=0;t<k;t++) A[k*n+j]-=A[k*n+t]*A[t*n+j]; A[j*n+k]=s[j]/A[k*n+k]; } }}//解方程LUx＝bvoidSolve_LU(double*A,intn,double*b,double*x){ inti,t; for(i=0;i<n;i++) { x[i]=b[i]; for(t=0;t<i;t++) x[i]-=A[i*n+t]*x[t]; } for(i=n-1;i>-1;i--) { for(t=i+1;t<n;t++) x[i]-=A[i*n+t]*x[t]; x[i]/=A[i*n+i]; }}//解線性方程組Ax＝BvoidSolve_lin(double*A,intn,double*B,double*x,intm)//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); Doolittle(A,n,M);//將A三角分解 for(i=0;i<m;i++) { for(j=0;j<n-1;j++) { temp=BT[i*n+j]; BT[i*n+j]=BT[i*n+M[j]]; BT[i*n+M[j]]=temp; }//將B轉(zhuǎn)置，使得同一方程組對應(yīng)的系數(shù)連續(xù)存儲 Solve_LU(A,n,BT+i*n,xT+i*n); } Transpose(xT,m,n,x);}//求n維向量V的無窮范數(shù)doubleVector_FanShu(double*V,intn){ inti; doublemax=0; for(i=0;i<n;i++) if(max<abs(V[i])) max=abs(V[i]); returnmax;}//求解非線性方程組voidSolve_non_Equation(doublea,doubleb,double*x){ doubleA[4][4],F[4],detx[4]; doubledet=1e-12; //求解精度要求 inti,k=0; intM=5000; //最大迭代次數(shù) //設(shè)定迭代初始值 x[0]=0.5; x[1]=1; x[2]=1; x[3]=1; do { Set_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;i++) x[i]+=detx[i]; k++; }while(k<M); printf("Newton法在該初值不收斂\n");}//查找n維向量V中與常數(shù)a最接近的元素的下標(biāo)intNear_Index(double*v,doublea,intn){ inti,Index; doublemin; min=abs(v[0]-a); Index=0; for(i=1;i<n;i++) { if(min>abs(v[i]-a)) { min=abs(v[i]-a); Index=i; } } returnIndex;}//求數(shù)表z(t,u)在點(diǎn)（x，y）處的分片二次代數(shù)差值doubleChaZhi(doublea,doubleb){ doublez[6][6]={-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}; doublex[6]={ 0, 0.2, 0.4, 0.6, 0.8, 1 }; doubley[6]={ 0, 0.4, 0.8, 1.2, 1.6, 2 }; doubleL[3],L_[3],Z; inti,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++) { L[k-i+1]=1; for(t=i-1;t<=i+1;t++) { if(t!=k) L[k-i+1]*=((a-x[t])/(x[k]-x[t])); } } 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-y[t])/(y[r]-y[t])); } } Z=0; for(k=i-1;k<=i+1;k++) { for(r=j-1;r<=j+1;r++) { Z+=(L[k-i+1]*L_[r-j+1]*z[k][r]); } } returnZ;}//對m×n數(shù)表U(x,y)進(jìn)行二元多項(xiàng)式擬合voidNiHe(double*U,double*x,double*y,intm,intn,intp,intq,double*C)//U為擬合數(shù)據(jù)點(diǎn)處的函數(shù)值{ inti,j; double*B,*BT,*BTB,*G,*GT,*GTG,*BTUG,*temp,*tempT,*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*)callo