数值分析计算实习题(三).doc

北航《数值分析》计算实习题目三
SY1004114 全昌彪
一、算法设计方案
1、解非线性方程组
首先将x，y当作已知的常数，求解四个未知数t，u，w，u。利用Newton法（简单迭代法不收敛）求解非线性方程组，得到与x，y对应的向量t，u。
求解步骤：
1）、选取初始向量{t,u,v,w}={1，1，1，1}；
2）、计算和；
3）、解关于的线性方程组

（调用Doolittle分解法求解此线性方程组）；
4）、若，则取；否则转5；
5）、计算；
题目中Newton法迭代公式为：

其中
2、分片二次代数插值
解题思路：
由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、子程序说明
子程序名称
功能
subroutine f_fit(t1,t2,c,sigma)
最小二乘法曲面拟合子程序，可给出拟合精度sigma
subroutine f_pxy(c,t1,t2,x,y,p_xy)
以x,y的幂函数为基，得到拟合系数矩阵C
subroutine f_zxy(z)
分片插值子函数，利用已知的（x,y）,得到z(x,y)
subroutine f_zut(u,t,p)
分片插值子函数，利用求取的（u,t）,得到z(u,t)
subroutine DLU(a,b,x)
Doolittle分解求线性方程组子函数
subroutine f_newton_iteration(x,y,u,t)
Newton迭代法解非线性方程组子程序
5、主程序main功能说明
主程序对xi,yi赋值，通过调用子程序对非线性方程组求解，得到相应的数据（t,u,v,w），通过调用插值子程序，得到对应的z=f(x,y),并以文件的形式进行输出。通过调用拟合子程序对拟合系数矩阵及拟合精度的求解，结果以文件形式输出。

二、fortran源程序
!/////曲面拟合子函数，并给出拟合精度/////
subroutine f_fit(t1,t2,c,sigma)
use imsl
implicit none
integer i,j,t1,t2
parameter n1=11
parameter n2=21
dimension b(n1,t1),b_trans(t1,n1),b_trans_b(t1,t1),b_inverse(t1,t1)
dimension g(n2,t2),g_trans(t2,n2),g_trans_g(t2,t2),g_inverse(t2,t2)
double precision b,g,b_trans,g_trans,b_trans_b,g_trans_g,b_inverse,g_inverse
dimension temp_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)
double precision temp_1,temp_2,temp_3,c,u,p_xy,sigma, x,y
!//////初始化x,y//////
do i=1,n1
x(i)=0.08*(i-1)
end do
do j=1,n2
y(j)=0.5+0.05*(j-1)
end do
!//////据题意，求出矩阵b,g//////
do i=1,n1
do j=1,t1
b(i,j)=x(i)**(j-1)
end do
end do

do i=1,n2
do j=1,t2
g(i,j)=y(i)**(j-1)
end do
end do
!//////确定b,g的转置b_trans和g_trans//////
do i=1,n1
do j=1,t1
b_trans(j,i)=b(i,j)
end do
end do
do i=1,n2
do j=1,t2
g_trans(j,i)=g(i,j)
end do
end do
!//////求解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矩阵//////
call f_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
do i=1,n1
do j=1,n2
call f_pxy(c,t1,t2,x(i),y(j),p_xy(i,j))
sigma=sigma+(u(i,j)-p_xy(i,j))**2
end do
end do
write (*,*) sigma
end subroutine f_fit

!/////子函数f_pxy给定拟合曲面的近似表达式/////
subroutine f_pxy(c,t1,t2,x,y,p_xy)
implicit none
integer t1,t2
dimension c(t1,t2)
double precision p_xy,temp_4,x,y,c
integer i,j,k
temp_4=0.0d0
do i=1,t1
do j=1,t2
temp_4=temp_4+c(i,j)*(x**(i-1))*(y**(j-1))
end do
end do
p_xy=temp_4 !拟合曲面
end subroutine f_pxy

!/////分片插值子函数，得出z=z(x,y)/////
subroutine f_zxy(z)
dimension x(11),y(21),u(11,21),t(11,21),z(11,21)
double precision x,y,u,t,z
integer i,j
do i=1,11
do j=1,21
x(i)=0.08*(i-1)
y(j)=0.5+0.05*(j-1)
!调用牛顿迭代法求非线性方程子函数
call f_newton_iteration(x(i),y(j),u(i,j),t(i,j))
call f_zut(u(i,j),t(i,j),z(i,j))
end do
end do
end subroutine f_zxy

!///// Doolittle分解子函数求解线性方程组/////
subroutine DLU(a,b,x)
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
integer i,j,k
do k=1,kkk,1
do j=k,kkk,1
m=0
do t=1,k-1,1
m=l(k,t)*u(t,j)+m
end do
u(k,j)=a(k,j)-m
end do
do i=k+1,kkk,1
m=0
do t=1,k-1,1
m=m+l(i,t)*u(t,k)
end do
l(i,k)=(a(i,k)-m)/u(k,k)
end do
l(k,k)=1
end do
!/////解方程//////
y(1)=b(1)
do i=1,kkk,1
m=0
do t=1,i-1,1
m=m+l(i,t)*y(t)
end do
y(i)=b(i)-m
end do
x(kkk)=y(kkk)/u(kkk,kkk)
do i=kkk-1,1,-1
m=0
do t=i+1,kkk,1
m=m+u(i,t)*x(t)
end do
x(i)=(y(i)-m)/u(i,i)
end do
end subroutine DLU

!/////牛顿迭代法子程序/////
subroutine f_newton_iteration(x,y,u,t)
parameter n=4
dimension f(n),aa(n,n),f_delta(n)
double precision f,aa,f_delta
double precision t,u,v,w,x,y
double precision epsion,s1,s2
integer i,j
!迭代初始值
t=1.0
u=1.0
v=1.0
w=1.0
epsion=1.0
do while(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)
call DLU(aa,f,f_delta)
s1=f_delta(1)*f_delta(1)+f_delta(2)*f_delta(2)+f_delta(3)*f_delta(3)+f_delta(4)*f_delta(4)
s1=sqrt(s1)
s2=t*t+u*u+v*v+w*w
s2=sqrt(s2)
epsion=s1/s2
t=t+f_delta(1)
u=u+f_delta(2)
v=v+f_delta(3)
w=w+f_delta(4)
end do
end subroutine f_newton_iteration

!!/////分片插值子函数，得出z=z(u,t)/////
subroutine f_zut(u,t,p)
implicit none
parameter n=6
dimension x(n),y(n)
double precision x,y,h,h1
double precision u,t
integer i,j,ii,jj !ii,jj 插值节点
integer k,r,ti
dimension lx(n),ly(n),z(n,n)
double precision p,lx,ly,z !p,lx,ly,z 插值多项式中的项
x(1)=0
y(1)=0
h=0.4
h1=0.2
do i=2,n
x(i)=x(1)+(i-1)*h
y(i)=y(1)+(i-1)*h1
end do
ii=0
jj=0
p=0
!/////给定z(n,n)/////
z(:,1)=(/-0.5,-0.34,0.14,0.94,2.06,3.5/)
z(:,2)=(/-0.42,-0.5,-0.26,0.3,1.18,2.38/)
z(:,3)=(/-0.18,-0.5,-0.5,-0.18,0.46,1.42/)
z(:,4)=(/0.22,-0.34,-0.58,-0.5,-0.1,0.62/)
z(:,5)=(/0.78,-0.02,-0.5,-0.66,-0.5,-0.02/)
z(:,6)=(/1.5,0.46,-0.26,-0.66,-0.74,-0.5/)
do i=3,n-2
if ((u.gt.(x(i)-h/2)).and.(u.le.(x(i)+h/2))) then
ii=i
goto 10
else if (u.le.(x(3)-h/2)) then
ii=2
else
ii=n-1
end if
end do
10 do j=3,n-2
if ((t.gt.(y(j)-h1/2)).and.(t.le.(y(j)+h1/2))) then
jj=j
goto 20
else if (t.le.(y(3)-h1/2)) then
jj=2
else
jj=n-1
end if
end do
20 do k=ii-1,ii+1
do r=jj-1,jj+1
lx(k)=1
ly(r)=1
do ti=ii-1,ii+1
if (ti.ne.k) then
lx(k)=lx(k)*(u-x(ti))/(x(k)-x(ti))
end if
end do
do ti=jj-1,jj+1
if (ti.ne.r) then
ly(r)=ly(r)*(t-y(ti))/(y(r)-y(ti))
end if
end do
p=p+lx(k)*ly(r)*z(k,r)
end do
end do
end subroutine f_zut

!/////主函数main/////
program main
implicit none
double precision vector_x(0:10),vector_y(0:20)
double precision xx(0:7),yy(0:4)
double precision tu(0:230,0:1),ut(0:1),zt(0:230),c(0:10,0:10),tu2(45,0:1),zt2(0:39),tt(45)
double precision z,sigma
integer i,j,m,n,l,t2,n2
n=10
m=20
do i=0,10
vector_x(i)=0.08*i
end do
do j=0,20
vector_y(j)=0.5+0.05*j
end do
open(11,file='xytuz.txt')
write(11,100)'x','y','t','u','z'
do i=0,10
do j=0,20
call f_newton_iteration(vector_x(i),vector_y(j),&
tu(i*21+j,1),tu(i*21+j,0))
call f_zut(tu(i*21+j,1),tu(i*21+j,0),zt(i*21+j))
write(11,200) vector_x(i),vector_y(j),&
tu(i*21+j,1),tu(i*21+j,0),zt(i*21+j)
end do
end do
close(11)
open(22,file='f(x,y).txt')
write(22,400)'xi','yi','f(xi,yi)'
do i=0,n
do j=0,m
write(22,300) vector_x(i),vector_y(j),zt(i*21+j)
end do
end do
close(22)
do l=1,6
call f_fit(l,l,c,sigma)
end do
open(33,file='matrix_c.txt')
do i=0,5
write(33,'(20e20.12)')(c(i,j),j=0,5)
end do
close(33)
t2=8
n2=5
do i=1,t2
xx(i-1)=0.1*i
end do
do j=1,n2
yy(j-1)=0.5+0.2*j
end do
do i=1,t2
do j=1,n2
call f_newton_iteration(xx(i-1),yy(j-1),tu2((i-1)*5+j,1), tu2((i-1)*5+j,0))
end do
end do
do i=1,40
call f_zut(tu2(i,1),tu2(i,0),zt2(i-1))
end do
do i=1,t2
do j=1,n2
call f_pxy(c,6,6,xx(i-1),yy(j-1),tt((i-1)*5+j))
end do
end do
open(44,file='f2(x,y).txt')
write(44,600)'x*i','y*i','f(xi,yi)','p(xi,yi)'
do i=1,t2
do j=1,n2
write(44,500) xx(i-1),yy(j-1),zt2((i-1)*5+j-1),tt((i-1)*5+j)
end do
end do
close(44)
! /////输出格式控制/////
100 format(2a10,3a15)
200 format(2e12.4,3e20.12)
300 format(2e12.4,e20.12)
400 format(2a10,a18)
500 format(2e20.7,2e20.12)
600 format(2a18,2a18)
end program main
三、结果汇总
1、数表：xi,yi、f(xi,yi)
xi yi f(xi,yi)
0.0000E+00 0.5000E+00 0.446504069241E+00
0.0000E+00 0.5500E+00 0.324683310597E+00
0.0000E+00 0.6000E+00 0.210159730631E+00
0.0000E+00 0.6500E+00 0.103043643055E+00
0.0000E+00 0.7000E+00 0.340192524000E-02
0.0000E+00 0.7500E+00 -0.887357886680E-01
0.0000E+00 0.8000E+00 -0.173371611924E+00
0.0000E+00 0.8500E+00 -0.250534594048E+00
0.0000E+00 0.9000E+00 -0.320276491493E+00
0.0000E+00 0.9500E+00 -0.382668052743E+00
0.0000E+00 0.1000E+01 -0.437795745638E+00
0.0000E+00 0.1050E+01 -0.485758921243E+00
0.0000E+00 0.1100E+01 -0.526667236102E+00
0.0000E+00 0.1150E+01 -0.560638463024E+00
0.0000E+00 0.1200E+01 -0.587796524654E+00
0.0000E+00 0.1250E+01 -0.608269768328E+00
0.0000E+00 0.1300E+01 -0.622189446194E+00
0.0000E+00 0.1350E+01 -0.629688384281E+00
0.0000E+00 0.1400E+01 -0.630899769389E+00
0.0000E+00 0.1450E+01 -0.625956164380E+00
0.0000E+00 0.1500E+01 -0.614988550047E+00
0.8000E-01 0.5000E+00 0.638015234587E+00
0.8000E-01 0.5500E+00 0.506611795876E+00
0.8000E-01 0.6000E+00 0.382176408075E+00
0.8000E-01 0.6500E+00 0.264863525560E+00
0.8000E-01 0.7000E+00 0.154780230633E+00
0.8000E-01 0.7500E+00 0.519927097057E-01
0.8000E-01 0.8000E+00 -0.434680179481E-01
0.8000E-01 0.8500E+00 -0.131601038093E+00
0.8000E-01 0.9000E+00 -0.212431072587E+00
0.8000E-01 0.9500E+00 -0.286004537561E+00
0.8000E-01 0.1000E+01 -0.352386059629E+00
0.8000E-01 0.1050E+01 -0.411655437858E+00
0.8000E-01 0.1100E+01 -0.463904893928E+00
0.8000E-01 0.1150E+01 -0.509236708669E+00
0.8000E-01 0.1200E+01 -0.547761104059E+00
0.8000E-01 0.1250E+01 -0.579594377202E+00
0.8000E-01 0.1300E+01 -0.604857251213E+00
0.8000E-01 0.1350E+01 -0.623673426166E+00
0.8000E-01 0.1400E+01 -0.636168257103E+00
0.8000E-01 0.1450E+01 -0.642467668455E+00
0.8000E-01 0.1500E+01 -0.642697116743E+00
0.1600E+00 0.5000E+00 0.840081396883E+00
0.1600E+00 0.5500E+00 0.699764166271E+00
0.1600E+00 0.6000E+00 0.566061444836E+00
0.1600E+00 0.6500E+00 0.439171614444E+00
0.1600E+00 0.7000E+00 0.319242167288E+00
0.1600E+00 0.7500E+00 0.206376218657E+00
0.1600E+00 0.8000E+00 0.100638546919E+00
0.1600E+00 0.8500E+00 0.206075984654E-02
0.1600E+00 0.9000E+00 -0.893540080136E-01
0.1600E+00 0.9500E+00 -0.173626954709E+00
0.1600E+00 0.1000E+01 -0.250799943998E+00
0.1600E+00 0.1050E+01 -0.320932252750E+00
0.1600E+00 0.1100E+01 -0.384097719188E+00
0.1600E+00 0.1150E+01 -0.440382160815E+00
0.1600E+00 0.1200E+01 -0.489881139368E+00
0.1600E+00 0.1250E+01 -0.532697954784E+00
0.1600E+00 0.1300E+01 -0.568941871351E+00
0.1600E+00 0.1350E+01 -0.598726545062E+00
0.1600E+00 0.1400E+01 -0.622168637157E+00
0.1600E+00 0.1450E+01 -0.639386546614E+00
0.1600E+00 0.1500E+01 -0.650499364497E+00
0.2400E+00 0.5000E+00 0.105151509248E+01
0.2400E+00 0.5500E+00 0.902927427685E+00
0.2400E+00 0.6000E+00 0.760580262813E+00
0.2400E+00 0.6500E+00 0.624715195878E+00
0.2400E+00 0.7000E+00 0.495519776775E+00
0.2400E+00 0.7500E+00 0.373134067685E+00
0.2400E+00 0.8000E+00 0.257656776698E+00
0.2400E+00 0.8500E+00 0.149150588766E+00
0.2400E+00 0.9000E+00 0.476470043700E-01
0.2400E+00 0.9500E+00 -0.468493081747E-01
0.2400E+00 0.1000E+01 -0.134356747670E+00
0.2400E+00 0.1050E+01 -0.214913334079E+00
0.2400E+00 0.1100E+01 -0.288573687009E+00
0.2400E+00 0.1150E+01 -0.355406352189E+00
0.2400E+00 0.1200E+01 -0.415491385199E+00
0.2400E+00 0.1250E+01 -0.468918240380E+00
0.2400E+00 0.1300E+01 -0.515783875727E+00
0.2400E+00 0.1350E+01 -0.556191070573E+00
0.2400E+00 0.1400E+01 -0.590246929357E+00
0.2400E+00 0.1450E+01 -0.618061557834E+00
0.2400E+00 0.1500E+01 -0.639746851961E+00
0.3200E+00 0.5000E+00 0.127124675843E+01
0.3200E+00 0.5500E+00 0.111500201761E+01
0.3200E+00 0.6000E+00 0.964607722308E+00
0.3200E+00 0.6500E+00 0.820347363242E+00
0.3200E+00 0.7000E+00 0.682447673529E+00
0.3200E+00 0.7500E+00 0.551085227171E+00
0.3200E+00 0.8000E+00 0.426392408546E+00
0.3200E+00 0.8500E+00 0.308463021403E+00
0.3200E+00 0.9000E+00 0.197357157387E+00
0.3200E+00 0.9500E+00 0.931056491500E-01
0.3200E+00 0.1000E+01 -0.428596469251E-02
0.3200E+00 0.1050E+01 -0.948339139383E-01
0.3200E+00 0.1100E+01 -0.178572979137E+00
0.3200E+00 0.1150E+01 -0.255553768875E+00
0.3200E+00 0.1200E+01 -0.325840141101E+00
0.3200E+00 0.1250E+01 -0.389506980643E+00
0.3200E+00 0.1300E+01 -0.446638198547E+00
0.3200E+00 0.1350E+01 -0.497324947822E+00
0.3200E+00 0.1400E+01 -0.541664031390E+00
0.3200E+00 0.1450E+01 -0.579756487909E+00
0.3200E+00 0.1500E+01 -0.611706299458E+00
0.4000E+00 0.5000E+00 0.149832107231E+01
0.4000E+00 0.5500E+00 0.133499864071E+01
0.4000E+00 0.6000E+00 0.117712512422E+01
0.4000E+00 0.6500E+00 0.102502405047E+01
0.4000E+00 0.7000E+00 0.878960016740E+00
0.4000E+00 0.7500E+00 0.739145103470E+00
0.4000E+00 0.8000E+00 0.605744887693E+00
0.4000E+00 0.8500E+00 0.478883880984E+00
0.4000E+00 0.9000E+00 0.358650648829E+00
0.4000E+00 0.9500E+00 0.245102261159E+00
0.4000E+00 0.1000E+01 0.138268376673E+00
0.4000E+00 0.1050E+01 0.381548905914E-01
0.4000E+00 0.1100E+01 -0.552527979027E-01
0.4000E+00 0.1150E+01 -0.141986870438E+00
0.4000E+00 0.1200E+01 -0.222094430862E+00
0.4000E+00 0.1250E+01 -0.295635227198E+00
0.4000E+00 0.1300E+01 -0.362679507525E+00
0.4000E+00 0.1350E+01 -0.423306161793E+00
0.4000E+00 0.1400E+01 -0.477601035622E+00
0.4000E+00 0.1450E+01 -0.525655428548E+00
0.4000E+00 0.1500E+01 -0.567564753485E+00
0.4800E+00 0.5000E+00 0.173189277935E+01
0.4800E+00 0.5500E+00 0.156203460147E+01
0.4800E+00 0.6000E+00 0.139721693052E+01
0.4800E+00 0.6500E+00 0.123780101007E+01
0.4800E+00 0.7000E+00 0.108408753012E+01
0.4800E+00 0.7500E+00 0.936322767079E+00
0.4800E+00 0.8000E+00 0.794704459103E+00
0.4800E+00 0.8500E+00 0.659387211130E+00
0.4800E+00 0.9000E+00 0.530487603091E+00
0.4800E+00 0.9500E+00 0.408088704433E+00
0.4800E+00 0.1000E+01 0.292244222111E+00
0.4800E+00 0.1050E+01 0.182982228522E+00
0.4800E+00 0.11