收藏
下载资源下载资源 加入VIP,免费下载

北航研究生数值分析作业第三题.docx

上传人：b****6 文档编号：15770628 上传时间：2023-07-07 格式：DOCX 页数：19 大小：97.74KB
下载 相关 举报
北航研究生数值分析作业第三题.docx_第1页
第1页 / 共19页
北航研究生数值分析作业第三题.docx_第2页
第2页 / 共19页
北航研究生数值分析作业第三题.docx_第3页
第3页 / 共19页
北航研究生数值分析作业第三题.docx_第4页
第4页 / 共19页
北航研究生数值分析作业第三题.docx_第5页
第5页 / 共19页
北航研究生数值分析作业第三题.docx_第6页
第6页 / 共19页
北航研究生数值分析作业第三题.docx_第7页
第7页 / 共19页
北航研究生数值分析作业第三题.docx_第8页
第8页 / 共19页
北航研究生数值分析作业第三题.docx_第9页
第9页 / 共19页
北航研究生数值分析作业第三题.docx_第10页
第10页 / 共19页
北航研究生数值分析作业第三题.docx_第11页
第11页 / 共19页
北航研究生数值分析作业第三题.docx_第12页
第12页 / 共19页
北航研究生数值分析作业第三题.docx_第13页
第13页 / 共19页
北航研究生数值分析作业第三题.docx_第14页
第14页 / 共19页
北航研究生数值分析作业第三题.docx_第15页
第15页 / 共19页
北航研究生数值分析作业第三题.docx_第16页
第16页 / 共19页
北航研究生数值分析作业第三题.docx_第17页
第17页 / 共19页
北航研究生数值分析作业第三题.docx_第18页
第18页 / 共19页
北航研究生数值分析作业第三题.docx_第19页
第19页 / 共19页
亲，该文档总共19页，全部预览完了，如果喜欢就下载吧！
下载资源
资源描述

北航研究生数值分析作业第三题.docx

《北航研究生数值分析作业第三题.docx》由会员分享，可在线阅读，更多相关《北航研究生数值分析作业第三题.docx（19页珍藏版）》请在冰点文库上搜索。

北航研究生数值分析作业第三题.docx

北航研究生数值分析作业第三题

北航研究生数值分析作业第三题

一、算法设计方案：

1.使用牛顿迭代法，对原题中给出的

，

，（

）的11*21组

分别求出原题中方程组的一组解，于是得到一组和

对应的

。

2.对于已求出的

，使用分片二次代数插值法对原题中关于

的数表进行插值得到

。

于是产生了z=f（x,y）的11*21个数值解。

3.从k=1开始逐渐增大k的值，并使用最小二乘法曲面拟合法对z=f（x,y）进行拟合，得到每次的

。

当

时结束计算，输出拟合结果。

4.计算

的值并输出结果，以观察

逼近

的效果。

其中

。

二、源代码：

#include"math.h"

#include"zuoye3.h"

#include"stdio.h"

voidmain（）

{

inti,j,k,r,s;//循环变量,其中k表示的是最小二乘法中x和y的阶数,不做他用

TUVWtuvw[11][21];//对应于x,y的非线性方程的解t,u,v,w

XYxy[11][21];//x,y的值以及f（x,y）,p（x,y）的值

XYstar_xy[8][5];//x*,y*的值以及f（x*,y*）,p（x*,y*）的值

doubleB[11][MAX_K];//关于x的系数矩阵

doubleG[21][MAX_K];//关于y的系数矩阵

doubleU[11][21];//要拟合的值

doubleC[MAX_K][MAX_K];//拟合后得到的系数矩阵

doubleWuCha[MAX_K];//k取不同值时对应的误差

FILE*fp=fopen（"1.txt","w"）;//文件指针

for（i=0;i<=10;i++）

for（j=0;j<=20;j++）

{

xy[i][j].x=0.08*i;//计算x

xy[i][j].y=0.5+0.05*j;//计算y

tuvw[i][j]=Resolve（xy[i][j].x,xy[i][j].y）;//解方程得到t,u

xy[i][j].f_xy=Function（tuvw[i][j].t,tuvw[i][j].u）;//利用插值函数得到f（x,y）

}

for（k=0;k

{

WuCha[k]=0;

for（i=0;i<11;i++）

for（r=0;r

B[i][r]=pow（xy[i][r].x,r）;//得到关于x的矩阵B

for（i=0;i<11;i++）

for（j=0;j<21;j++）

U[i][j]=xy[i][j].f_xy;//得到要拟合的值矩阵U

for（j=0;j<21;j++）

for（s=0;s

G[j][s]=pow（xy[s][j].y,s）;//得到关于y的矩阵G

Calculate（C,B,U,G,k）;//计算系数矩阵C

for（i=0;i<=10;i++）

for（j=0;j<=20;j++）

{

xy[i][j].p_xy=P_nihe（C,xy[i][j].x,xy[i][j].y,k）;//计算p（x,y）

WuCha[k]+=pow（xy[i][j].p_xy-xy[i][j].f_xy,2）;//计算误差

}

if（WuCha[k]<=1e-7）break;//误差是否满足要求

}

//计算x*，y*，f（x*,y*）,p（x*,y*）

for（i=0;i<8;i++）

for（j=0;j<5;j++）

{

star_xy[i][j].x=0.1*（i+1）;

star_xy[i][j].y=0.5+0.2*（j+1）;

tuvw[0][0]=Resolve（star_xy[i][j].x,star_xy[i][j].y）;

star_xy[i][j].f_xy=Function（tuvw[0][0].t,tuvw[0][0].u）;

star_xy[i][j].p_xy=P_nihe（C,star_xy[i][j].x,star_xy[i][j].y,k）;

}

//输出Xi，Yj，f（Xi，Yj）到文件

fprintf（fp,"\n-----|----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|\n"）;

fprintf（fp,"|"）;

for（i=0;i<=5;i++）fprintf（fp,"%.2f|",xy[i][0].x）;

fprintf（fp,"\n-----|----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|\n"）;

for（j=0;j<=20;j++）

for（i=0;i<=5;i++）

{

if（i==0）fprintf（fp,"%.2f|",xy[i][j].y）;

if（xy[i][j].f_xy>0）fprintf（fp,"%.12e|",xy[i][j].f_xy）;

elsefprintf（fp,"%.12e|",xy[i][j].f_xy）;

if（i==5）

{

fprintf（fp,"\n-----|----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|\n"）;

}

}

fprintf（fp,"|"）;

for（i=6;i<=10;i++）fprintf（fp,"%.2f|",xy[i][0].x）;

fprintf（fp,"\n-----|----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|"）;

fprintf（fp,"----------------------|\n"）;

for（j=0;j<=20;j++）

for（i=6;i<=10;i++）

{

if（i==6）fprintf（fp,"%.2f|",xy[i][j].y）;

if（xy[i][j].f_xy>0）fprintf（fp,"%.12e|",xy[i][j].f_xy）;

elsefprintf（fp,"%.12e|",xy[i][j].f_xy）;

if（i==10）

{

fprintf（fp,"\n-----|----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|"）;

fprintf（fp,"----------------------|\n"）;

}

}

//输出选择过程的k值和误差值

for（i=0;i<=k;i++）fprintf（fp,"k=%d时,误差值为:

%.12e\n",i,WuCha[i]）;

//输出达到精度要求的系数矩阵Crs

for（i=0;i<=k;i++）

for（j=0;j<=k;j++）

{

if（C[i][j]>0）fprintf（fp,"%.12e\t",C[i][j]）;

elsefprintf（fp,"%.12e\t",C[i][j]）;

if（j==k）fprintf（fp,"\n"）;

}

//输出x*，y*，f（x*,y*）,p（x*,y*）

fprintf（fp,"\n-----|----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|"）;

fprintf（fp,"----------------------|\n"）;

fprintf（fp,"|"）;

for（j=0;j<5;j++）fprintf（fp,"%.2f|",star_xy[0][j].y）;

fprintf（fp,"\n-----|----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|"）;

fprintf（fp,"----------------------|\n"）;

for（i=0;i<8;i++）

for（j=0;j<5;j++）

{

if（j==0）fprintf（fp,"%.2f|",star_xy[i][j].x）;

if（star_xy[i][j].f_xy>0）fprintf（fp,"%.12e|",star_xy[i][j].f_xy）;

elsefprintf（fp,"%.12e|",star_xy[i][j].f_xy）;

if（j==4）

{

fprintf（fp,"\n|"）;

for（r=0;r<=4;r++）

{

if（star_xy[i][r].p_xy>0）fprintf（fp,"%.12e|",star_xy[i][r].p_xy）;

elsefprintf（fp,"%.12e|",star_xy[i][r].p_xy）;

}

fprintf（fp,"\n-----|----------------------|----------------------|"）;

fprintf（fp,"----------------------|----------------------|"）;

fprintf（fp,"----------------------|\n"）;

}

}

fclose（fp）;

}

TUVWResolve（doublex,doubley）

{

doublec1=2.67+x;

doublec2=1.07+y;

doublec3=3.74+x;

doublec4=0.79+y;//计算4个方程的常数项

doubletuvw[N]={0,1,1,1,1};//初始迭代向量

doubleF[N],derta[N];//非线性方程对应的泛函值,迭代差值

doubleDF[N][N]={{0,0,0,0,0},{0,-0,1,1,1},{0,1,-0,1,1},{0,0.5,1,-0,1},{0,1,0.5,1,-0}};

//F'（x）系数矩阵

inti;

for（;;）//开始迭代

{

F[1]=0.5*cos（tuvw[1]）+tuvw[2]+tuvw[3]+tuvw[4]-c1;

F[2]=tuvw[1]+0.5*sin（tuvw[2]）+tuvw[3]+tuvw[4]-c2;

F[3]=0.5*tuvw[1]+tuvw[2]+cos（tuvw[3]）+tuvw[4]-c3;

F[4]=tuvw[1]+0.5*tuvw[2]+tuvw[3]+sin（tuvw[4]）-c4;

DF[1][1]=-sin（tuvw[1]）*0.5;

DF[2][2]=0.5*cos（tuvw[2]）;

DF[3][3]=-sin（tuvw[3]）;

DF[4][4]=cos（tuvw[4]）;

Gauss（&DF[0][0],F,derta,5）;//GAUSS法求迭代差值

for（i=1;i

if（FindMax（derta）/FindMax（tuvw）<=MIS）break;//判断是满足精度要求

}

TUVWtemp={tuvw[1],tuvw[2],tuvw[3],tuvw[4]};

returntemp;//返回求解结果

}

voidGauss（double*A,constdouble*b,double*result,constintM）

{

intk,i,j;//GAUSS法会改变系数矩阵和等号右边的向量，使用时要格外注意

intMaxLineNum;

doublemax;

doubletemp;

doubleMik;

double*b1=newdouble[M];//b的备份

double*A1=newdouble[M*M];//A的备份

for（i=1;i

b1[i]=b[i];//创建b的备份，以防止解方程组时修改b

for（i=1;i

for（j=1;j

A1（i,j）=A（i,j）;//创建A的备份，以防止解方程组时修改A

for（k=1;k

{

max=fabs（A1（k,k））;

MaxLineNum=k;

for（i=k;i

{

if（fabs（A1（i,k））>max）

{

max=fabs（A1（i,k））;

MaxLineNum=i;

}

}

if（k!

=MaxLineNum）

{

for（j=k;j

{

temp=A1（k,j）;

A1（k,j）=A1（MaxLineNum,j）;

A1（MaxLineNum,j）=temp;

}

temp=b1[k];

b1[k]=b1[MaxLineNum];

b1[MaxLineNum]=temp;

}

for（i=k+1;i

{

Mik=A1（i,k）/A1（k,k）;

for（j=k+1;j

A1（i,j）-=Mik*A1（k,j）;

b1[i]-=Mik*b1[k];

}

}

for（k=M-1;k>=1;k--）

{

if（k==M-1）result[k]=b1[k]/A1（k,k）;

else

{

for（j=k+1;j

b1[k]-=A1（k,j）*result[j];

result[k]=b1[k]/A1（k,k）;

}

}

delete[]b1;

delete[]A1;

}

doubleFindMax（constdoublea[N]）//求无穷范数

{

doublemax=fabs（a[1]）;

for（inti=1;i

if（max

returnmax;

}

doubleFunction（doublet,doubleu）//分片二次插值函数

{

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}};

inti,j;

switch（（int）（10*t））//判断插值的中心点x坐标

{

case0:

case1:

case2:

i=1;break;

case3:

case4:

i=2;break;

case5:

case6:

i=3;break;

case7:

case8:

case9:

i=4;break;

default:

printf（"error"）;break;

}

switch（（int）（5*u））//判断插值中心点y坐标

{

case0:

case1:

case2:

j=1;break;

case3:

case4:

j=2;break;

case5:

case6:

j=3;break;

case7:

case8:

case9:

j=4;break;

default:

printf（"error"）;break;

}

doublevalue=0;

intk,r;

for（k=i-1;k<=i+1;k++）

for（r=j-1;r<=j+1;r++）

value+=z[k][r]*Lagrange_t（t,i,k）*Lagrange_u（u,j,r）;

//分片二次插值

returnvalue;

}

doubleLagrange_t（constdoublex,constinti,constintk）

{

intt;

doublevalue=1;

for（t=i-1;t<=i+1;t++）

if（t!

=k）value*=（x-t*0.2）/（k*0.2-t*0.2）;

returnvalue;

}

doubleLagrange_u（constdoubley,constintj,constintr）

{

intt;

doublevalue=1;

for（t=j-1;t<=j+1;t++）

if（t!

=r）value*=（y-t*0.4）/（r*0.4-t*0.4）;

returnvalue;

}

voidCalculate（double（&C）[MAX_K][MAX_K],double（&B）[11][MAX_K],

double（&U）[11][21],double（&G）[21][MAX_K],constintk）

{//系数矩阵C=（BT*B）-1*BT*U*G*（GT*G）-1

inti,j,t;

doubleBTB[MAX_K][MAX_K];

doubleBTB_1[MAX_K][MAX_K];

doubleGTG[MAX_K][MAX_K];

doubleGTG_1[MAX_K][MAX_K];

doubleBTU[MAX_K][21];

doubleBTUG[MAX_K][MAX_K];

doubleBTB_1_BTUG[MAX_K][MAX_K];

for（i=0;i

for（j=0;j

BTB[i][j]=GTG[i][j]=（i==j?

1:

0）;

for（i=0;i<=k;i++）

for（j=0;j<=k;j++）

{

BTB[i][j]=GTG[i][j]=0;

for（t=0;t<11;t++）BTB[i][j]+=B[t][i]*B[t][j];

for（t=0;t<21;t++）GTG[i][j]+=G[t][i]*G[t][j];

}

Mat_Inverse（BTB,BTB_1）;//求逆

Mat_Inverse（GTG,GTG_1）;//求逆

//将各个矩阵相乘

for（i=0;i<=k;i++）

for（j=0;j<21;j++）

{

BTU[i][j]=0;

for（t=0;t<11;t++）BTU[i][j]+=B[t][i]*U[t][j];

}

for（i=0;i<=k;i++）

for（j=0;j<=k;j++）

{

BTUG[i][j]=0;

for（t=0;t<21;t++）BTUG[i][j]+=BTU[i][t]*G[t][j];

}

for（i=0;i<=k;i++）

for（j=0;j<=k;j++）

{

BTB_1_BTUG[i][j]=0;

for（t=0;t<=k;t++）BTB_1_BTUG[i][j]+=BTB_1[i][t]*BTUG[t][j];

}

for（i=0;i<=k;i++）

for（j=0;j<=k;j++）

{

C[i][j]=0;

for（t=0;t<=k;t++）C[i][j]+=BTB_1_BTUG[i][t]*GTG_1[t][j];

}

}

voidMat_Inverse（double（&Mat）[MAX_K][MAX_K],double（&Mat_1）[MAX_K][MAX_K]）

{

inti,j;

doubleA[MAX_K+1][MAX_K+1];

doubleb[MAX_K+1];

doubleresult[MAX_K+1];

for（i=1;i

for（j=1;j

A[i][j]=Mat[i-1][j-1];

for（i=1;i

{

for（j=1;j

1:

0）;

Gauss（&A[0][0],b,result,MAX_K+1）;

for（j=1;j

展开阅读全文
相关资源
猜你喜欢
相关搜索
资源标签

当前位置：首页 > 幼儿教育 > 幼儿读物

copyright@ 2008-2023 冰点文库 网站版权所有

经营许可证编号:鄂ICP备19020893号-2