北航数值分析大作业3

一、算法设计方案

1.使用牛顿迭代法，对原题中给出的i x i 08.0=，j y j 05.05.0+=，

(010

,020i j ≤≤≤≤)的11*21组j i y x ,分别求出原题中方程组的一组解，于是得到一组和i i y x ,对应的j i t u ,。

2.对于已求出的j i t u ,，使用分片二次代数插值法对原题中关于u t z ,,的数表进行插值得到

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

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

10-<σ时结束计算，输出拟合结果。

4.计算)5,,2,1,8,,2,1)(,(),,(*

***⋅⋅⋅=⋅⋅⋅=j i y x p y x f j i j i 的值并输出结果，以观察),(y x p 逼近),(y x f 的效果。其中j y i x j i 2.05.0,1.0*

*+==。

二、算法实现方案 1、求(,)f x y ：

（1）Newton 法解非线性方程组

0.5cos 2.670.5sin 1.07(1)0.5cos 3.740.5sin 0.79

t u v w x t u v w y t u v w x t u v w y +++-=⎧⎪+++-=⎪

⎨

+++-=⎪⎪+++-=⎩， 其中，t, u, v ,w 为待求的未知量，x, y 为代入的已知量。

设(,,,)T t u v w ξ=，给定精度水平12110ε-=和最大迭代次数M ，则解该线性方程组的迭代格式为：

*(0)(0)(0)(0)(0)(k+1)

()()1()(,,,)()()0,1,T k k k t u v w F F k ξξξ

ξξξ-⎧=⎪'=-⎨⎪=

⎩

在附近选取初值， 迭代终止条件为()(1)

()

1/k k k ξξ

ξε-∞

∞

-≤，若k M >时仍未达到迭代精度，则迭代计算失

败。

其中，雅可比矩阵

0.5*cos(t) + u + v + w - x - 2.67t + 0.5*sin(u) + v + w - y - 1.07()0.5*t + u + cos(v) + w - x - 3.74t + 0.5*u + v + sin(w) - y - 0.79F ξ⎡⎤⎢⎥

⎢

⎥=⎢⎥⎢⎥⎣⎦

，

0.5*sin()

11110.5*cos()11()0.51sin()110.51cos()t u F v w ξ-⎡⎤⎢⎥⎢⎥'=⎢⎥-⎢⎥

⎣⎦

， （2）分片双二次插值：

根据题目给出的表格，

(0,1,,5)(0,1,

,5)

i i t ih i u j j τ= == =（其中，0.2h τ= =0.4，

） 对于给定的(,)t u ,如果(,)t u 满足

,322

,322

i i j j h h

t t t i u u u j τ

τ

-

<≤+ 2≤≤-

<≤+ 2≤≤

则选择(,)(1,,1;1,,1)k r t u k i i i r j j j =-+=-+为插值节点，相应的插值多项式为

1

111

(,)()()(,)j i k

r

k

r

k i r j h t u l t l u g t u ++=-=-=

(2) ∑∑

其中，

1

11

1()(1,,1)()(1,,1)i m

k m i k m

m k

j n

r n j r n

n r

t t l t k i i i t t u u l u r j j j u u +=-≠+=-≠-= =-+--=

=-+-∏∏

如果1422h h t t t t ≤+

>+或,则在式（2）中取i=1或i=4; 如果1422

u u u ττ

≤+>+或u ,则在式（2）中取u=1或u=4。

在区域{(,)|00.8,0.5 1.5}D x y x y =≤≤≤≤上，将

(,)i j x y (0.08,i x i =i=0,1,…,10;j y =0.5+0.2j,j=0,1,…20)代入到非线性方程组（1）中，用

Newton 法解出,,(,)i j i j t u ，再由分片双二次插值得,,,(,)i j i j i j z h t u =，则有,(,)i j i j z f x y =（i=0,1,…,10;j=0,1,…,20），即求出了(,)z f x y =。 2、求(,)p x y ：

乘积型最小二乘拟合曲面：

（1）求系数矩阵C ：

11()()T T T C B B B UG G G --=

其中，

2

00021

112101010111k k k x x x x x x B x x x ⎡⎤⎢

⎥⎢

⎥

=⎢⎥⎢

⎥⎢⎥⎣⎦ 200021

112

20

2020111k k k y y y y y y G y x y ⎡⎤

⎢

⎥⎢

⎥

=⎢⎥⎢⎥⎢⎥⎣⎦

0,0

0,10,201,0

1,11,20,10,0

10,1

10,20,(,)(0,1,

,10;0,1,,20)i j i j

z z z z z z U z f x y i j z z z ⎡⎤⎢⎥⎢

⎥== ==⎢⎥⎢⎥⎣⎦

计算中涉及到对矩阵求逆，接着在后面将会具体说明列主元的高斯消去法求矩阵的逆的方法。

（2）确定最小的k 值，拟合曲面：

,0

(,)k

r s rs

r s p x y c

x y ==

∑

设10

20

200

((,)(,))i j i j i j f x y p x y σ

===-∑∑，给定精度水平7210ε-=和最大迭代次数N ，

则确定最小k 值的迭代格式为：

(),0

1020

()()2

00(,)(0,1,,10;0,1,

,20)

((,)(,))

0,1,2,k

k r s

i j rs i j r s k k i j i j i j p x y c x y i j f x y p x y k σ===⎧= ==⎪⎪

⎪=-⎨⎪

⎪=⎪⎩

∑

∑∑

迭代终止条件为()

2k σ

ε≤，若k N >时仍未达到迭代精度，则迭代计算失败。

待确定满足精度条件的最小k 值后，就可以进行曲面拟合计算了。 3、关于列主元的高斯消去法求矩阵的逆：

设非奇异矩阵n n A C ⨯∈，且1

B A -= ，则

AB I =，

对B 和I 列分块，有