北航数值分析A大作业Word下载.docx
《北航数值分析A大作业Word下载.docx》由会员分享，可在线阅读，更多相关《北航数值分析A大作业Word下载.docx（10页珍藏版）》请在冰豆网上搜索。
的值并输出结果，以观察
逼近
的效果。
。
二、全部源程序
2e|"
f[j][i]）;
printf（"
\n"
）;
}
printf（"
-------------------------------------------------------------------------------------\n"
}
k=0;
printf（"
不同k对应的精度"
\n----------------------------------------------"
do{
C=（double*）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）;
else
{
\n故可知k=k=%d<
%.0e，满足题设要求"
det）;
break;
}while（++k<
11）;
2e"
C[i*（k+1）+j]）;
----------------------------------------------------------------------------\n"
2e|%+.12e|%f|\n"
*i,+*j,d,p,abs（d-p））;
-----------------------------------------------------------------------\n\n\n"
}
6e方阵奇异\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＝b
voidSolve_LU（double*A,intn,double*b,double*x）
{
inti,t;
for（i=0;
i<
n;
i++）
{
x[i]=b[i];
for（t=0;
i;
t++）x[i]-=A[i*n+t]*x[t];
for（i=n-1;
i>
-1;
i--）
for（t=i+1;
x[i]-=A[i*n+t]*x[t];
x[i]/=A[i*n+i];
//解线性方程组Ax＝B
voidSolve_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三角分解
m;
for（j=0;
j<
n-1;
temp=BT[i*n+j];
BT[i*n+j]=BT[i*n+M[j]];
BT[i*n+M[j]]=temp;
}//将B转置，使得同一方程组对应的系数连续存储
Solve_LU（A,n,BT+i*n,xT+i*n）;
Transpose（xT,m,n,x）;
//求n维向量V的无穷范数
doubleVector_FanShu（double*V,intn）
inti;
doublemax=0;
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;
//最大迭代次数
//设定迭代初始值
x[0]=;
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;
4;
i++）
x[i]+=detx[i];
k++;
}while（k<
M）;
Newton法在该初值不收敛\n"
//查找n维向量V中与常数a最接近的元素的下标
intNear_Index（double*v,doublea,intn）
inti,Index;
doublemin;
min=abs（v[0]-a）;
Index=0;
for（i=1;
if（min>
abs（v[i]-a））
min=abs（v[i]-a）;
Index=i;
returnIndex;
//求数表z（t,u）在点（x，y）处的分片二次代数差值
doubleChaZhi（doublea,doubleb）
doublez[6][6]={,,,,,,
,,,,,
,,,,};
doublex[6]={0,,,,,1};
doubley[6]={0,,,,,2};
doubleL[3],L_[3],Z;
inti,j,k,r,t;
i=Near_Index（x,a,6）;
j=Near_Index（y,b,6）;
//插值区域边界处插值节点的选取
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++）
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;
=r）L_[r-j+1]*=（（b-y[t]）/（y[r]-y[t]））;
Z=0;
for（r=j-1;
Z+=（L[k-i+1]*L_[r-j+1]*z[k][r]）;
returnZ;
//对m×
n数表U（x,y）进行二元多项式拟合
voidNiHe（double*U,double*x,double*y,intm,intn,intp,intq,double*C）
//U为拟合数据点处的函数值
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*）calloc（p*q,sizeof（double））;
temp=（double*）calloc（p*n,sizeof（double））;
p;
j++）B[i*p+j]=pow（x[i],j）;
q;
j++）G[i*q+j]=pow（y[i],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）;
free（B）;
free（BT）;
free（G）;
free（GT）;
free（temp）;
temp=（double*）calloc（p*q,sizeof（double））;
Solve_lin（BTB,p,BTUG,temp,q）
升级会员 