数值分析第3次上机作业

《数值分析》

第3次上机作业

姓名：

学号：

2015.12.9

1 算法设计思路和方案

1.1 总体方案设计

（1）解非线性方程组。将给定的(,)i i x y 当作已知量代入题目给定的非线性方程组，求得与(,)i i x y 相对应的数组t[i][j],u[i][j]。

（2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=(,)i i f x y 。

（3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。

（4）观察和(,)i i p x y 的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,)i i f x y ，再与对应的(,)i i p x y 比较即可，这里求解(,)i i p x y 可以直接使用（3）中的C[r][s]和k 。

1.2 具体算法设计：

1.2.1 解非线性方程组

牛顿法解方程组()0F x =的解*x ，可采用如下算法：

（1）在*x 附近选取(0)x D ∈，给定精度水平0ε>和最大迭代次数M 。 （2）对于0,1,

k M =执行

①计算()()k F x 和()()k F x '。 ②求解关于()k x ∆的线性方程组

()()()()()k k k F x x F x '∆=-

③若()

()

k k x x ε∞

∞

∆≤，则取*()k x x ≈，并停止计算；否则转④。

④计算(1)()()k k k x x x +=+∆。

⑤若k M <，则继续，否则，输出M 次迭代不成功的信息，并停止计算。

1.2.2 分片双二次插值

给定已知数表以及需要插值的节点，进行分片二次插值的算法：

设已知数表中的点为：00(0,1,,)

(0,1,,)

i j x x ih i n y y j j m τ=+=⎧⎪⎨

=+=⎪⎩ ，需要插值的节点为(,)x y 。 （1）根据(,)x y 选择插值节点(,)i j x y ： 若12h x x ≤+

或12

n h

x x ->-，插值节点对应取1i =或1i n =-， 若12y y τ≤+或12

n y y τ

->-，插值节点对应取1j =或1i m =-。

若,2222,2222

i i j j h h x x x i n y y y j m ττ⎧-<≤+≤≤-⎪⎪⎨⎪-<≤+≤≤-⎪⎩

则选择(,)(1,,1;1,,1)k r x y k i i i r j j j =-+=-+为插值节点。

（2）计算

1

11

1()(1,,1)()(1,,1)

i t

k t i k t

t k j t

r t j r t

t r

x x l x k i i i x x y y l y r j j j y y +=-≠+=-≠-==-+--==-+-∏

∏

插值多项式的公式为：

1

111

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

r

k

r

k i r j p x y l x l y f x y ++=-=-=

∑∑

注：本步进行插值运算的是(,)t u ，利用(,)i j x y 与(,)t u 的对应关系就可以得到z 与(,)i j x y 的对应关系。 1.2.3 曲面拟合

根据插值得到的数表,,(,)i j i j x y f x y 进行曲面拟合的过程： （1）根据拟合节点和基底函数写出矩阵B 和G ：

010*********

()()()()()()()()()k k k n

n n x x x x x x B x x x ⎛⎫

⎪ ⎪= ⎪ ⎪ ⎪⎝⎭ 01

0000

11110

1

()()()()()()()()()k k k m

m m y y y y y y G y y y ⎛⎫

⎪ ⎪

=

⎪ ⎪ ⎪⎝⎭

（2）计算 11()()T T T C B B B UG G G --=。

（3）对于每一个(,)i j x y ，

*

00(,)()()k

k

r s i j rs i j r s p x y C x y ===∑∑

拟合需要达到的精度条件为：

*2700

[(,)]10n m

i j ij i j p x y u σ-===-≤∑∑

其中ij u 对应着插值得到的数表,,(,)i j i j x y f x y 中(,)i j f x y 的值。

让k 逐步增加，每一次重复执行以上几步，直到710σ-≤ 成立。此时的k 值就是要求解最小的k 。

2 源程序代码

#include

#include

#include

#define n 4

static double x[n],x_[n],t[11][21],u[11][21],U[11][21],tt[8][5],uu[8][5],UU[8][5],c[9][9]; int k;

//求数组中的最大值

double MAX(double a[n])

{

int i;

double max;

max=fabs(a[0]);

for(i=0;i

if(fabs(a[i])>max)

max=fabs(a[i]);

return(max);

}

//选主元的DooLittle 分解法求线性方程组

void DooLittle(double a[n][n],double b[n])

{

int i,j,k,t,i_k,m[n];

double u[n][n],l[n][n],s[n],y[n];

double u_x,l_u,max,temp;

for(k=1;k<=n;k++)

{

for(i=k;i<=n;i++)

{

l_u=0;

for(t=1;t<=k-1;t++)

l_u=l_u+l[i-1][t-1]*u[t-1][k-1];

s[i-1]=a[i-1][k-1]-l_u;

}

max=fabs(s[k-1]);

i_k=k;

m[k-1]=k;

for(i=k;i<=n;i++)

if(fabs(s[i-1])>max)

{

max=fabs(s[i-1]);

i_k=i;

m[k-1]=i_k;

}