数值分析B第二次大作业Word文档下载推荐.docx
- 文档编号:4085560
- 上传时间:2023-05-02
- 格式:DOCX
- 页数:31
- 大小:42.13KB
数值分析B第二次大作业Word文档下载推荐.docx
《数值分析B第二次大作业Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《数值分析B第二次大作业Word文档下载推荐.docx(31页珍藏版)》请在冰点文库上搜索。
#definesss1.0e-219
#defineL500//迭代最大次数;
#defineK10//确定最小二乘法拟合时的最多次数;
doubleaaa_u[6]={0,0.4,0.8,1.2,1.6,2.0},//输入矩阵数据表在z(u,t)
aaa_t[6]={0,0.2,0.4,0.6,0.8,1.0},
aaa_ztu[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},
},
aaa_val_z[11][21]={0},
init_tuvw[4]={1,2,1,2},
aaa_crs[K][K]={0};
intvaluez(doublex,doubley,inti,intj)//用牛顿迭代法计算x(i),y(i)对应的函数值z
{
doublevtuvw[4],det[4],ff[4][4],f[4],mo_k,mo_delta_k,shang,t,u;
intn=0,k,flag_max,flag_stat;
voidsolve(doubleff[4][4],doubledet[4],doublef[4]);
doubleSharding(double,double);
flag_stat=1;
for(k=0;
k<
M;
k++)vtuvw[k]=init_tuvw[k];
do
{
f[0]=-1.0*(0.5*cos(vtuvw[0])+vtuvw[1]+vtuvw[2]+vtuvw[3]-x-2.67);
//输入非线性方程组
f[1]=-1.0*(vtuvw[0]+0.5*sin(vtuvw[1])+vtuvw[2]+vtuvw[3]-y-1.07);
f[2]=-1.0*(0.5*vtuvw[0]+vtuvw[1]+cos(vtuvw[2])+vtuvw[3]-x-3.74);
f[3]=-1.0*(vtuvw[0]+0.5*vtuvw[1]+vtuvw[2]+sin(vtuvw[3])-y-0.79);
ff[0][0]=-0.5*sin(vtuvw[0]);
ff[0][1]=1;
ff[0][2]=1;
ff[0][3]=1;
//对非线性方程组计算关于t,u,v,w的偏导数
ff[1][0]=1;
ff[1][1]=0.5*cos(vtuvw[1]);
ff[1][2]=1;
ff[1][3]=1;
ff[2][0]=0.5;
ff[2][1]=1;
ff[2][2]=-1*sin(vtuvw[2]);
ff[2][3]=1;
ff[3][0]=1;
ff[3][1]=0.5;
ff[3][2]=1;
ff[3][3]=cos(vtuvw[3]);
solve(ff,det,f);
flag_max=0;
for(k=1;
k++)if(fabs(det[k])>
fabs(det[flag_max]))flag_max=k;
mo_delta_k=det[flag_max];
k++)if(fabs(vtuvw[k])>
fabs(vtuvw[flag_max]))flag_max=k;
mo_k=vtuvw[flag_max];
shang=fabs(mo_delta_k/mo_k);
if(shang<
ss)
t=vtuvw[0]+det[0],u=vtuvw[1]+det[1];
flag_stat=0;
break;
}
else
for(k=0;
4;
k++)vtuvw[k]+=det[k];
while(n++<
=L);
if(!
flag_stat)
printf("
x%d=%f,y%d=%f:
f(x%d,y%d)=%.12le\n"
i,0.08*i,j,0.5+0.05*j,i,j,aaa_val_z[i][j]=Sharding(u,t));
}
returnflag_stat;
}
voidsolve(doubleff[4][4],doubledet[4],doublef[4])//牛顿消元法解线性方程组
inti,j,k,ik;
doublehe,mik;
for(k=0;
k<
4-1;
k++)
ik=k;
for(i=k;
i<
i++)
if(fabs(ff[ik][k])<
fabs(ff[i][k]))ik=i;
for(j=k;
j<
j++)
he=ff[k][j];
ff[k][j]=ff[ik][j];
ff[ik][j]=he;
he=f[k];
f[k]=f[ik];
f[ik]=he;
for(i=k+1;
i<
i++)
mik=ff[i][k]/ff[k][k];
j++)
ff[i][j]=ff[i][j]-mik*ff[k][j];
f[i]-=mik*f[k];
det[4-1]=f[4-1]/ff[4-1][4-1];
for(k=4-2;
k>
=0;
k--)
he=0;
for(j=k+1;
j<
j++)he+=ff[k][j]*det[j];
det[k]=(f[k]-he)/ff[k][k];
doubleSharding(doubleu,doublet)//分片法进行代数二次插值
intr,i,j,k;
doubleduiyzh_u[3],duiyzh_t[3],z;
z=0;
r=int(fabs((t/0.2)+0.5));
k=int(fabs((u/0.4)+0.5));
if(r=0)r=1;
if(r=5)r=4;
if(k=0)k=1;
if(k=5)k=4;
duiyzh_u[0]=(u-aaa_u[k])*(u-aaa_u[k+1])/(aaa_u[k-1]-aaa_u[k])/(aaa_u[k-1]-aaa_u[k+1]);
duiyzh_u[1]=(u-aaa_u[k-1])*(u-aaa_u[k+1])/(aaa_u[k]-aaa_u[k-1])/(aaa_u[k]-aaa_u[k+1]);
duiyzh_u[2]=(u-aaa_u[k-1])*(u-aaa_u[k])/(aaa_u[k+1]-aaa_u[k-1])/(aaa_u[k+1]-aaa_u[k]);
duiyzh_t[0]=(t-aaa_t[r])*(t-aaa_t[r+1])/(aaa_t[r-1]-aaa_t[r])/(aaa_t[r-1]-aaa_t[r+1]);
duiyzh_t[1]=(t-aaa_t[r-1])*(t-aaa_t[r+1])/(aaa_t[r]-aaa_t[r-1])/(aaa_t[r]-aaa_t[r+1]);
duiyzh_t[2]=(t-aaa_t[r-1])*(t-aaa_t[r])/(aaa_t[r+1]-aaa_t[r-1])/(aaa_t[r+1]-aaa_t[r]);
for(i=r-1;
i<
=r+1;
for(j=k-1;
=k+1;
j++)z+=duiyzh_t[i-r+1]*duiyzh_u[j-k+1]*aaa_ztu[i][j];
returnz;
intkchoose()//计算k的值,以及相应的误差
intx,y,i,j,k,kvalue;
doublevec_x[11],vec_y[21],aaa_btb[K][K]={0},aaa_gtg[K][K]={0},aaa_btbbt[K][11],aaa_btbbtu[K][21],he,wucha;
voidinv(doubleaaarix[K][K],int);
kvalue=1;
//赋初值
选择过程中的k与截断误差:
\n"
);
for(i=0;
11;
i++)vec_x[i]=0.08*i;
for(j=0;
21;
j++)vec_y[j]=0.5+0.05*j;
do
kvalue;
j++)
{
k++)he+=pow(vec_x[k],i)*pow(vec_x[k],j);
aaa_btb[i][j]=he;
inv(aaa_btb,kvalue);
k++)he+=pow(vec_y[k],i)*pow(vec_y[k],j);
aaa_gtg[i][j]=he;
inv(aaa_gtg,kvalue);
k++)he+=aaa_btb[i][k]*pow(vec_x[j],k);
aaa_btbbt[i][j]=he;
k++)he+=aaa_btbbt[i][k]*aaa_val_z[k][j];
aaa_btbbtu[i][j]=he;
k++)he+=aaa_btbbtu[i][k]*pow(vec_y[k],j);
k++)he+=aaa_btb[i][k]*aaa_gtg[k][j];
aaa_crs[i][j]=he;
wucha=0;
for(x=0;
x<
x++)
for(y=0;
y<
y++)
he+=aaa_crs[i][j]*pow(0.08*x,i)*pow(0.5+0.05*y,j);
wucha+=(he-aaa_val_z[x][y])*(he-aaa_val_z[x][y]);
%d%.12le\n\n"
kvalue,wucha);
if(wucha<
1.0e-7)
Crs的系数矩阵:
Crs[%d][%d]=\n{\n"
kvalue,kvalue);
{"
kvalue-1;
j++)printf("
%.12le,"
aaa_crs[i][j]);
%.12le"
aaa_crs[i][kvalue-1]);
},\n"
}\n\n"
returnkvalue;
while(kvalue++<
K);
return0;
}
voidinv(doubleaaarix[K][K],intrank)
inti,j,t,n;
doubleaaa_he[K][K]={0},aaarix_bak[K][K],vec_he[K],vec_x[K]={0},vec_dx[K]={0},vec_r[K]={0},he,max_x,max_dx;
voidpreprocess(double[K][K],double[K],int);
//对矩阵进行处理,使其条件数变小。
voidLUbreakdown(double[K][K],int);
//LU分解法
voidLUsolving(double[K][K],double[K],double[K],int);
//LU分解法求矩阵的逆
n=0;
rank;
j++)aaarix_bak[i][j]=aaa_he[i][j]=aaarix[i][j];
vec_he[i]=0;
LUbreakdown(aaa_he,rank);
if(i==0)vec_he[i]=1;
else
vec_he[i-1]=0;
vec_he[i]=1;
LUsolving(aaa_he,vec_x,vec_he,rank);
while(n++<
=3)
for(t=0;
t<
t++)he+=aaarix_bak[j][t]*vec_x[t];
vec_r[j]=vec_he[j]-he;
LUsolving(aaa_he,vec_dx,vec_r,rank);
max_x=vec_x[0],max_dx=vec_dx[0];
if(vec_x[j]>
max_x)max_x=vec_x[j];
if(vec_dx[j]>
max_dx)max_dx=vec_dx[j];
j++)vec_x[j]+=vec_dx[j];
j++)aaarix[j][i]=vec_x[j];
voidpreprocess(doubleaaarix[K][K],doublevec_he[K],intrank)/对矩阵进行处理,使其条件数变小。
inti,k;
doublehe;
for(i=0;
rank;
i++)
he=aaarix[i][0];
k<
k++)if(aaarix[i][k]>
he)he=aaarix[i][k];
k++)aaarix[i][k]/=he;
vec_he[i]/=he;
voidLUbreakdown(doubleaaarix[K][K],intrank)//LU分解法
inti,j,k,t;
k++)
he=0;
for(t=0;
t<
k;
t++)
he+=aaarix[i][t]*aaarix[t][k];
aaarix[i][k]-=he;
if(k<
rank-1){
for(j=k+1;
j++){
he+=aaarix[k][t]*aaarix[t][j];
aaarix[k][j]=(aaarix[k][j]-he)/aaarix[k][k];
voidLUsolving(doubleaaarix[K][K],doublevec_x[K],doublevec_he[K],intrank)//LU分解法求矩阵的逆
inti,t;
vec_x[0]=vec_he[0]/aaarix[0][0];
for(i=1;
i;
t++)he+=aaarix[i][t]*vec_x[t];
vec_x[i]=(vec_he[i]-he)/aaarix[i][i];
for(i=rank-2;
i>
i--)
for(t=i+1;
vec_x[i]-=he;
voidfvalue(inthe){
inti,j;
近似表达式:
8;
for(j=0;
5;
he+=valuez(0.1*(i+1),0.5+0.2*(j+1),i,j);
voidshuchu(){//输出数表
inti,j,x,y,k;
doublet;
数表:
x*,y*,f(x*,y*),p(x*,y*):
x++){
for(y=0;
y++){
t=0;
for(i=0;
k;
i++){
t=t+aaa_crs[i][j]*pow(0.1*(x+1),i)*pow(0.5+0.2*(y+1),j);
x[%d]=%f,y[%d]=%f:
f(x*[%d],y*[%d])=%.12le,p(x*[%d],y*[%d])=%.12le\n"
x+1,0.1*(x+1),y+1,0.5+0.2*(y+1),x+1,y+1,aaa_val_z[x][y],x+1,y+1,t);
voidmain()//主函数
inti,j,k,he=0;
=10;
i++)//计算xi,yi对应的值
for(j=0;
j<
=20;
he+=valuez(0.08*i,0.5+0.05*j,i,j);
k=kchoose();
//计算选择过程中的k值,以及相应的误差值
fvalue(he);
//输出近似表达式的每个值
he)
shuchu();
//输出数表
三、输出结果:
xi,yi,f(xi,yi);
x0=0.000000,y0=0.500000:
f(x0,y0)=4.465040184799e-001
x0=0.000000,y1=0.550000:
f(x0,y1)=3.246832629274e-001
x0=
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 第二次 作业
![提示](https://static.bingdoc.com/images/bang_tan.gif)