北京航空航天大学2020届研究生
《数值分析》实验作业
第九题
院系：xx学院
学号：
姓名：
2020年11月
Q9:方程组A.4
一、 算法设计方案
（一）总体思路
1.题目要求∑∑===
k i k
j s r rs
y x c
y x p 00
),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。
),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。
2.),(*
*
j i y x f 与1使用相同方法求得，),(*
*
j i y x p 由计算得出的p(x,y)直接带入),(*
*
j i y x 求得。
1. ),(i i y x 与),(i i y x f 的数表的获得
对区域D ={ （x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ，yj=1+0.008j ，得到),(i i y x 共31×21组。将每组带入A4方程组，即可获得五个二元函数组，通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。再将
),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。
2.乘积型最小二乘曲面拟合
2.1使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下：
⎪⎪⎪⎭⎫ ⎝⎛=k i i k x x x x B 0000 ， ⎪⎪⎪⎭⎫ ⎝⎛=k j j
k y y y y G 0000
数表矩阵如下：
⎪⎪⎪⎭
⎫
⎝
⎛=),(),(),(),(0000j i i j y x f y x f y x f y x f U
记C=[rs c ]，则系数rs c 的表达式矩阵为：
11-)(-=G G UG B B B C T T
T ）（
通过求解如下线性方程，即可得到系数矩阵C 。
UG B G G C B B T T T =)()(
2.2计算),(),,(*
***j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值
),(**j i y x f 的计算与),(j i y x f 相同。将),(**j i y x 代入原方程组，求解响应)
,(*
*ij ij u t 进行分片双二次插值求得),(**j i y x f 。),(*
*j
i y x p 的计算则可以直接将),(**j i y x 代入所求p(x,y)。
二、 源程序
********* 第三次数值分析大作业Q9************ integer::i, j, K1, L1, n, m
dimension X(31), Y(21), T(6), U(6), Z(6, 6), UX(11, 21), TY(11, 21), FXY(11, 21), C(6, 6) dimension z1(31, 21), z2(31, 21), z3(31, 21), z4(31, 21), z5(31, 21) dimension X1(8), Y1(5), FXY1(8, 5), PXY1(8, 5), UX1(8, 5), TY1(8, 5)
real(8) X, Y, T, U, Z, FXY, UX, TY, C, Error, X1, Y1, FXY1, PXY1, UX1, TY1, z1, z2, z3, z4, z5
OPEN(1, FILE = '数值分析Q9程序结果.TXT')
do i = 1, 31
X(i) = 1 + 0.008 * (i - 1)
end do
do j = 1, 21
Y(j) = 1 + 0.008 * (j - 1) !生成矩阵
end do
!*****求解非线性方程组，得到数表(X(i), Y(j), f1, f2, f3, f4, f5)*******
do i = 1, 31
do j = 1, 21
call Newton(X(i), Y(j), z1(i, j), z2(i, j), z3(i, j), z4(i, j), z5(i, j)) !用离散牛顿法求解非线性方程组。
if (i == 27.and.j == 21)then
!write(*, *)z1(i, j), z2(i, j), z3(i, j), z4(i, j), z5(i, j)
!pause
end if
end do
end do
write(*, "('数表( X(i),Y(j),f1,f2,f3,f4,f5 ):')")
write(*, "(3X,'X',7X,'Y',10X,'f1(x,y)',13x,'f2(x,y)',13x,'f3(x,y)',13x,'f4(x,y)',13x,'f5(x,y)')")
do i = 1, 31
do j = 1, 21
write(1, '(1X,F5.3,2X,F5.3,2X,5E20.13)') X(i), Y(j), z1(i, j), z2(i, j), z3(i, j), z4(i, j), z5(i, j) write(*, '(3X,I2,1X,I2,3X,F5.3,2X,F5.3,2X,5E20.13)')i, j, X(i), Y(j), z1(i, j), z2(i, j), z3(i, j), z4(i, j), z5(i, j)
end do
end do
!***********最小二乘拟合得到P（x，y）**************
N = 31
M = 21
WRITE(1, '(" ","K和σ分别为：")')
WRITE(*, '(" ","选择过程的K1,L1和σ分别为：")')
DO K1 = 1, 8
DO L1 = 1, 8
