北航数值分析A大作业.docx
《北航数值分析A大作业.docx》由会员分享，可在线阅读，更多相关《北航数值分析A大作业.docx（22页珍藏版）》请在冰点文库上搜索。
北航数值分析A大作业
一、算法设计方案
1、解非线性方程组
将各拟合节点（xi，yj）分别带入非线性方程组，求出与
相对应的数组te[i][j]，ue[i][j]，求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。
2、二元二次分偏插值
对数表z（t,u）进行分片二次代数插值，求得对应（tij，uij）处的值，即为
的值。
根据给定的数表，可将整个插值区域分成16个小的区域，故先判断tij,uij所在，的区域，再作此区域的插值，计算zij，相应的Lagrange形式的插值多项式为：
其中
（k=m-1,m,m+1）
（r=n-1,n,n+1）
3、曲面拟合
从k=1开始逐渐增大k的值，使用最小二乘法曲面拟合法对z=f（x,y）进行拟合，当
时结束计算。
拟合基函数φr（x）ψs（y）选择为φr（x）=xr，ψs（y）=ys。
拟合系数矩阵c通过连续两次解线性方程组求得。
，
其中
，
4、观察比较
计算
的值并输出结果，以观察
逼近
的效果。
其中
。
二、全部源程序
//hean.cpp:
定义控制台应用程序的入口点。
//
#include"stdafx.h"
#include
#include
#include
voidSet_non_JacobiA（double*A,double*x）;//求题中非线性方程组对应自变量向量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＝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）的值
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数值分析第三次大作业\n"）;
//打印数表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对应的精度"）;
printf（"\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;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对应的σ=%.12e",k,d）;
if（d>=det）
free（C）;
else
{
printf（"\n----------------------------------------------"）;
printf（"\n故可知k=k=%d<%.0e，满足题设要求",det）;
break;
}
}while（++k<11）;
//打印拟合系数矩阵c
printf（"\n\n\n当k=%d时拟合函数p（x,y）的系数矩阵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观察以下各点的逼近情况\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"）;
}
//求题中非线性方程组对应自变量向量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迭代法的右端式:
-F（x）
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;
for（i=0;ifor（j=0;j{
C[i*n+j]=0;
for（k=0;kC[i*n+j]+=A[i*s+k]*B[k*n+j];
}
}
//将m×n阶矩阵A的转置为AT
voidTranspose（double*A,intm,intn,double*AT）
{
inti,j;
for（i=0;ifor（j=0;jAT[j*m+i]=A[i*n+j];
}
//对n阶矩阵A进行选主元的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{
for（i=k;i{
s[i]=A[i*n+k];
for（t=0;t}
Maxs=abs（s[k]）;M[k]=k;
for（i=k+1;i{
if（Maxs{
Maxs=abs（s[i]）;
M[k]=i;
}
}
if（M[k]!
=k）
{
for（t=0;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{
for（t=0;tA[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{
x[i]=b[i];
for（t=0;t}
for（i=n-1;i>-1;i--）
{
for（t=i+1;tx[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三角分解
for（i=0;i{
for（j=0;j{
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;
for（i=0;iif（maxreturnmax;
}
//求解非线性方程组
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）return;
for（i=0;i<4;i++）
x[i]+=detx[i];
k++;
}while（kprintf（"Newton法在该初值不收敛\n"）;
}
//查找n维向量V中与常数a最接近的元素的下标
intNear_Index（double*v,doublea,intn）
{
inti,Index;
doublemin;
min=abs（v[0]-a）;
Index=0;
for（i=1;i{
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]={-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;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数表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））;
for（i=0;i{
for（j=0;j}
for（i=0;i{
for（j=0;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（temp）;
free（CT）;
free（tempT）;
free（GTG）;
}
//求多项式拟合函数P（x,y）在点（x，y）处的函数值
doublePxy（double*C,intp,intq,doublex,doubley）
{
inti,j;
d
升级会员 