数值分析计算实习题三Word格式文档下载.docx
《数值分析计算实习题三Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《数值分析计算实习题三Word格式文档下载.docx（26页珍藏版）》请在冰豆网上搜索。
解题思路：
由1得到的x，y和t，u的映射表，f（t（x，y），u（x，y）），即求得f（x，y）。
但由于得到的t，u不可能正好是题目提供的二维数表中的值，需要用相关规则对插值节点加以规范。
利用（x，y）以及对应的f（x，y），就可能通过二元拉格朗日插值多项式得到f（x，y）的表达式。
插值节点：
1）、根据计算得到的t、u值，选取插值节点；
选择标准如下：
假设对（t，u），这里用（x，y）代替：
设：
a）、若满足：
则应选择
为插值节点
b）、若满足：
或
2）、双元二次插值子程序
相应的插值多项式为：
3、最小二乘法曲面拟合
设在三维直角坐标系
中给定（m+1）*（n+1）个点（即三维坐标）
在本题中即为
。
选定M+1个x的函数
以及N+1个y的函数
本题中
，
，于是得到乘积型基函数构成的曲面，
随着k值的不断增大，精度
会越来越大，题目要求精度为
，此时的k即为要求的最小值。
1）、求解矩阵A
固定
，以
为基函数对数据
作最小二乘拟合，得到n+1条拟合曲线
是法方程
的解，而
，求解n+1线性方程组，得到矩阵A。
2）、求解矩阵G
3）、系数矩阵C
4、子程序说明
子程序名称
功能
subroutinef_fit（t1,t2,c,sigma）
最小二乘法曲面拟合子程序，可给出拟合精度sigma
subroutinef_pxy（c,t1,t2,x,y,p_xy）
以x,y的幂函数为基，得到拟合系数矩阵C
subroutinef_zxy（z）
分片插值子函数，利用已知的（x,y）,得到z（x,y）
subroutinef_zut（u,t,p）
分片插值子函数，利用求取的（u,t）,得到z（u,t）
subroutineDLU（a,b,x）
Doolittle分解求线性方程组子函数
subroutinef_newton_iteration（x,y,u,t）
Newton迭代法解非线性方程组子程序
5、主程序main功能说明
主程序对xi,yi赋值，通过调用子程序对非线性方程组求解，得到相应的数据（t,u,v,w），通过调用插值子程序，得到对应的z=f（x,y）,并以文件的形式进行输出。
通过调用拟合子程序对拟合系数矩阵及拟合精度的求解，结果以文件形式输出。
二、fortran源程序
!
/////曲面拟合子函数，并给出拟合精度/////
useimsl
implicitnone
integeri,j,t1,t2
parametern1=11
parametern2=21
dimensionb（n1,t1）,b_trans（t1,n1）,b_trans_b（t1,t1）,b_inverse（t1,t1）
dimensiong（n2,t2）,g_trans（t2,n2）,g_trans_g（t2,t2）,g_inverse（t2,t2）
doubleprecisionb,g,b_trans,g_trans,b_trans_b,g_trans_g,b_inverse,g_inverse
dimensiontemp_1（t1,n1）,temp_2（t1,n2）,temp_3（t1,t2）,&
c（t1,t2）,u（n1,n2）,p_xy（n1,n2）,x（n1）,y（n2）
doubleprecisiontemp_1,temp_2,temp_3,c,u,p_xy,sigma,x,y
//////初始化x,y//////
doi=1,n1
x（i）=0.08*（i-1）
enddo
doj=1,n2
y（j）=0.5+0.05*（j-1）
enddo
//////据题意，求出矩阵b,g//////
doj=1,t1
b（i,j）=x（i）**（j-1）
enddo
doi=1,n2
doj=1,t2
g（i,j）=y（i）**（j-1）
//////确定b,g的转置b_trans和g_trans//////
b_trans（j,i）=b（i,j）
g_trans（j,i）=g（i,j）
//////求解b_trans_b和g_trans_g的逆矩阵b_inverse和g_inverse//////
b_trans_b=matmul（b_trans,b）!
matmul为矩阵相乘函数，库函数
g_trans_g=matmul（g_trans,g）
b_inverse=.i.b_trans_b
g_inverse=.i.g_trans_g
//////求解c矩阵//////
callf_zxy（u）
temp_1=matmul（b_inverse,b_trans）
temp_2=matmul（temp_1,u）
temp_3=matmul（temp_2,g）
c=matmul（temp_3,g_inverse）
/////求拟合精度误差/////
sigma=0
doj=1,n2
callf_pxy（c,t1,t2,x（i）,y（j）,p_xy（i,j））
sigma=sigma+（u（i,j）-p_xy（i,j））**2
write（*,*）sigma
endsubroutinef_fit
/////子函数f_pxy给定拟合曲面的近似表达式/////
integert1,t2
dimensionc（t1,t2）
doubleprecisionp_xy,temp_4,x,y,c
integeri,j,k
temp_4=0.0d0
doi=1,t1
temp_4=temp_4+c（i,j）*（x**（i-1））*（y**（j-1））
p_xy=temp_4!
拟合曲面
endsubroutinef_pxy
/////分片插值子函数，得出z=z（x,y）/////
dimensionx（11）,y（21）,u（11,21）,t（11,21）,z（11,21）
doubleprecisionx,y,u,t,z
integeri,j
doi=1,11
doj=1,21
调用牛顿迭代法求非线性方程子函数
callf_newton_iteration（x（i）,y（j）,u（i,j）,t（i,j））
callf_zut（u（i,j）,t（i,j）,z（i,j））
endsubroutinef_zxy
/////Doolittle分解子函数求解线性方程组/////
integer,parameter:
:
kkk=4
real（kind=8）m
real（kind=8）,dimension（kkk,kkk）,intent（in）:
a
real（kind=8）,dimension（kkk,kkk）:
l=0,u=0
real（kind=8）,dimension（kkk）:
y
real（kind=8）,dimension（kkk）,intent（out）:
x
real（kind=8）,dimension（kkk）,intent（in）:
b
dok=1,kkk,1
doj=k,kkk,1
m=0
dot=1,k-1,1
m=l（k,t）*u（t,j）+m
u（k,j）=a（k,j）-m
doi=k+1,kkk,1
m=m+l（i,t）*u（t,k）
l（i,k）=（a（i,k）-m）/u（k,k）
l（k,k）=1
/////解方程//////
y
（1）=b
（1）
doi=1,kkk,1
dot=1,i-1,1
m=m+l（i,t）*y（t）
y（i）=b（i）-m
x（kkk）=y（kkk）/u（kkk,kkk）
doi=kkk-1,1,-1
dot=i+1,kkk,1
m=m+u（i,t）*x（t）
x（i）=（y（i）-m）/u（i,i）
endsubroutineDLU
/////牛顿迭代法子程序/////
parametern=4
dimensionf（n）,aa（n,n）,f_delta（n）
doubleprecisionf,aa,f_delta
doubleprecisiont,u,v,w,x,y
doubleprecisionepsion,s1,s2
迭代初始值
t=1.0
u=1.0
v=1.0
w=1.0
epsion=1.0
dowhile（epsion.ge.1e-12）
f
（1）=-（0.5*cos（t）+u+v+w-x-2.67）
f
（2）=-（t+0.5*sin（u）+v+w-y-1.07）
f（3）=-（0.5*t+u+cos（v）+w-x-3.74）
f（4）=-（t+0.5*u+v+sin（w）-y-0.79）
aa（1,1）=-0.5*sin（t）
aa（1,2）=1
aa（1,3）=1
aa（1,4）=1
aa（2,1）=1
aa（2,2）=0.5*cos（u）
aa（2,3）=1
aa（2,4）=1
aa（3,1）=0.5
aa（3,2）=1
aa（3,3）=-sin（v）
aa（3,4）=1
aa（4,1）=1
aa（4,2）=0.5
aa（4,3）=1
aa（4,4）=cos（w）
callDLU（aa,f,f_delta）
s1=
升级会员 