数值分析计算实习题三.docx
《数值分析计算实习题三.docx》由会员分享,可在线阅读,更多相关《数值分析计算实习题三.docx(26页珍藏版)》请在冰豆网上搜索。
数值分析计算实习题三
北航《数值分析》计算实习题目三
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、子程序说明
子程序名称
功能
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源程序
!
/////曲面拟合子函数,并给出拟合精度/////
subroutinef_fit(t1,t2,c,sigma)
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//////
doi=1,n1
doj=1,t1
b(i,j)=x(i)**(j-1)
enddo
enddo
doi=1,n2
doj=1,t2
g(i,j)=y(i)**(j-1)
enddo
enddo
!
//////确定b,g的转置b_trans和g_trans//////
doi=1,n1
doj=1,t1
b_trans(j,i)=b(i,j)
enddo
enddo
doi=1,n2
doj=1,t2
g_trans(j,i)=g(i,j)
enddo
enddo
!
//////求解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
doi=1,n1
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
enddo
enddo
write(*,*)sigma
endsubroutinef_fit
!
/////子函数f_pxy给定拟合曲面的近似表达式/////
subroutinef_pxy(c,t1,t2,x,y,p_xy)
implicitnone
integert1,t2
dimensionc(t1,t2)
doubleprecisionp_xy,temp_4,x,y,c
integeri,j,k
temp_4=0.0d0
doi=1,t1
doj=1,t2
temp_4=temp_4+c(i,j)*(x**(i-1))*(y**(j-1))
enddo
enddo
p_xy=temp_4!
拟合曲面
endsubroutinef_pxy
!
/////分片插值子函数,得出z=z(x,y)/////
subroutinef_zxy(z)
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
x(i)=0.08*(i-1)
y(j)=0.5+0.05*(j-1)
!
调用牛顿迭代法求非线性方程子函数
callf_newton_iteration(x(i),y(j),u(i,j),t(i,j))
callf_zut(u(i,j),t(i,j),z(i,j))
enddo
enddo
endsubroutinef_zxy
!
/////Doolittle分解子函数求解线性方程组/////
subroutineDLU(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
integeri,j,k
dok=1,kkk,1
doj=k,kkk,1
m=0
dot=1,k-1,1
m=l(k,t)*u(t,j)+m
enddo
u(k,j)=a(k,j)-m
enddo
doi=k+1,kkk,1
m=0
dot=1,k-1,1
m=m+l(i,t)*u(t,k)
enddo
l(i,k)=(a(i,k)-m)/u(k,k)
enddo
l(k,k)=1
enddo
!
/////解方程//////
y
(1)=b
(1)
doi=1,kkk,1
m=0
dot=1,i-1,1
m=m+l(i,t)*y(t)
enddo
y(i)=b(i)-m
enddo
x(kkk)=y(kkk)/u(kkk,kkk)
doi=kkk-1,1,-1
m=0
dot=i+1,kkk,1
m=m+u(i,t)*x(t)
enddo
x(i)=(y(i)-m)/u(i,i)
enddo
endsubroutineDLU
!
/////牛顿迭代法子程序/////
subroutinef_newton_iteration(x,y,u,t)
parametern=4
dimensionf(n),aa(n,n),f_delta(n)
doubleprecisionf,aa,f_delta
doubleprecisiont,u,v,w,x,y
doubleprecisionepsion,s1,s2
integeri,j
!
迭代初始值
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=