北航数值分析A大作业Word文件下载.docx
《北航数值分析A大作业Word文件下载.docx》由会员分享，可在线阅读，更多相关《北航数值分析A大作业Word文件下载.docx（22页珍藏版）》请在冰点文库上搜索。
voidArray_Mult_Array（double*A,double*B,intm,ints,intn,double*C）;
//矩阵相乘AB＝C
voidTranspose（double*A,intm,intn,double*AT）;
//转置
voidDoolittle（double*A,intn,int*M）;
voidSolve_LU（double*A,intn,double*b,double*x）;
//解方程LUx＝b
voidSolve_lin（double*A,intn,double*B,double*x,intm）;
//解线性方程组Ax＝B
doubleVector_FanShu（double*V,intn）;
voidSolve_non_Equation（doublea,doubleb,double*x）;
//求解非线性方程组
intNear_Index（double*v,doublea,intn）;
//查找n维向量V中与常数a最接近的元素的下标
doubleChaZhi（doublea,doubleb）;
//求数表z（t,u）在点（x，y）处的分片二次代数差值
voidNiHe（double*U,double*x,double*y,intm,intn,intp,intq,double*C）;
//对m×
n数表U（x,y）进行二元多项式拟合
doublePxy（double*C,intp,intq,doublex,doubley）;
//求多项式拟合函数P（x,y）在点（x，y）处的函数值
voidmain（）
{
doublenon[4];
//非线性方程组变量向量（t,u,v,w）
doublef[11][21];
//f（x,y）在拟合节点处的值
doublex[11],y[21];
//拟合节点
double*C;
//拟合系数矩阵
doubledet=1e-7;
//拟合精度要求
doublep,d;
inti,j,k,t;
//设置拟合节点
for（i=0;
i<
11;
i++）x[i]=0.08*i;
for（j=0;
j<
21;
j++）y[j]=0.5+0.05*j;
//求拟合节点处的f（x,y）的值
i++）
{
for（j=0;
j++）
{
Solve_non_Equation（x[i],y[j],non）;
f[i][j]=ChaZhi（non[0],non[1]）;
}
}
printf（"
\n数值分析第三次大作业\n"
）;
//打印数表f（xi,yi）
\nf（xi,yi）\n"
for（t=0;
t<
21/3;
t++）
printf（"
-------------------------------------------------------------------------------------\n"
|f（x,y）|"
for（i=3*t;
3*t+3;
i++）
printf（"
y=%4.2f|"
y[i]）;
\n"
j++）
|x=%4.2f|"
x[j]）;
for（i=3*t;
printf（"
%+.12e|"
f[j][i]）;
k=0;
不同k对应的精度"
\n----------------------------------------------"
do{
C=（double*）calloc（（k+1）*（k+1）,sizeof（double））;
NiHe（*f,x,y,11,21,k+1,k+1,C）;
//求在当前k值拟合的节点处误差
d=0;
for（i=0;
for（j=0;
{
p=Pxy（C,k+1,k+1,x[i],y[j]）;
d+=（（f[i][j]-p）*（f[i][j]-p））;
}
\nk=%d对应的σ=%.12e"
k,d）;
if（d>
=det）
free（C）;
else
\n故可知k=k=%d<
%.0e，满足题设要求"
det）;
break;
}while（++k<
11）;
//打印拟合系数矩阵c
\n\n\n当k=%d时拟合函数p（x,y）的系数矩阵c：
k）;
----------------------------------------------------------------------------\n"
（k+1）/3;
（k+1）;
for（j=3*t;
%+.12e"
C[i*（k+1）+j]）;
//打印（x*,y*）处的拟合逼近情况
\n\n观察以下各点的逼近情况\n"
|（x*,y*）|f（x*,y*）|p（x*,y*）||f-p||\n"
-----------------------------------------------------------------------\n"
for（i=1;
=8;
for（j=1;
=5;
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）;
|（%3.1f,%3.1f）|%+.12e|%+.12e|%f|\n"
0.1*i,0.5+0.2*j,d,p,abs（d-p））;
-----------------------------------------------------------------------\n\n\n"
}
//求题中非线性方程组对应自变量向量x的雅克比矩阵
voidSet_non_JacobiA（double*A,double*x）//A为*4阶系数矩阵，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迭代法的右端式:
voidSet_non_B（double*B,double*x,doublea,doubleb）//a非线性方程的参数，对应题中的x
//b非线性方程的参数，对应题中的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;
m;
n;
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的转置为AT
voidTranspose（double*A,intm,intn,double*AT）
inti,j;
AT[j*m+i]=A[i*n+j];
//对n阶矩阵A进行选主元的Doolittle分解
voidDoolittle（double*A,intn,int*M）//将分解得到的L、U仍然存储在原来的矩阵A中
double*s;
doubleMaxs,temp;
s=（double*）calloc（n,sizeof（double））;
for（k=0;
k++）
{
for（i=k;
s[i]=A[i*n+k];
for（t=0;
k;
t++）s[i]-=A[i*n+t]*A[t*n+k];
Maxs=abs（s[k]）;
M[k]=k;
for（i=k+1;
if（Maxs<
abs（s[i]））
{
Maxs=abs（s[i]）;
M[k]=i;
if（M[k]!
=k）
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;
%.16e方阵奇异\n"
Maxs）;
A[k*n+k]=s[k];
for（j=k+1;
（j<
n）&
&
（k<
n-1）;
t++）
A[k*n+j]-=A[k*n+t]*A[t*n+j];
A[j*n+k]=s[j]/A[k*n+k];
}
voidSolve_LU（double*A,intn,double*b,double*x）
inti,t;
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];
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三角分解
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]=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;
4;
x[i]+=detx[i];
k++;
}while（k<
M）;
Newton法在该初值不收敛\n"
intNear_Index（double*v,doublea,intn）
inti,Index;
doublemin;
min=abs（v[0]-a）;
Index=0;
if（min>
abs（v[i]-a））
min=abs（v[i]-a）;
Index=i;
returnIndex;
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）;
//插值区域边界处插值节点的选取
if（i==0）i=1;
if（i==5）i=4;
if（j==0）j=1;
if（j==5）j=4;
for（k=i-1;
=i+1;
L[k-i+1]=1;
for（t=i-1;
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;
voidNiHe（double*U,double*x,double*y,intm,intn,intp,intq,double*C）
//U为拟合数据点处的函数值
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）;
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（CT）;
free（tempT）;
free（GTG）;
doublePxy（double*C,intp,intq,doublex,doubley）
d
升级会员 