北航研究生数值分析A作业三

《数值分析（A）》计算实习题目三

一、题目

关于x，y，t，u，v，w的下列方程组

0.5cost+u+v+w-x=2.67

t+0.5sinu+v+w-y=1.07

0.5t+u+cosv+w-x=3.74

t+0.5u+v+sinw-y=0.79

以及关于z，t，u的下列二维数表

确定了一个二元函数z=f(x，y)。

1．试用数值方法求出f(x，y)在区域D={(x,y)︱0≤x≤0.8,0.5≤y≤1.5}上的一个

近似表达式

,0

(,)k

r s

rs r s p x y c x y

==

∑

要求p(x,y)一最小的k 值达到以下的精度

10

20

27

((,)(,))10i

j

i j i j f x y

p x y σ-===

-≤∑∑

其中x i =0.08i ，y j =0.5+0.05j 。

2．计算f(x i *,y j *),p(x i *,y j *)(i=1,2,…,8;j=1,2,…,5)的值，以观察p(x,y)逼近f(x,y)

的效果，其中x i *=0.1i,y j *

=0.5+0.2j 。

二、算法设计方案

1.将0.08(0,1,,10)

i x i i *

== 和0.50.05(0,1,,20)j y j j *=+= 代入非线性方程

组中，用牛顿法解出i t 和j u ；

2.以采取分片二次插值，选择（m ，n ）满足

,232

2

m i m h h t t t m -<≤+

≤≤ ,23

2

2

n j n u u u n τ

τ

-

<≤+

≤≤

如果12

i

h t t ≤+

或42

i

h t t >-

，则m=1或4；如果12

j

u u τ

≤+

或42

j

u u τ

>-

，则n=1

或n=4。选择(,)(1,,1;1,,1)

k r t u k

m m m r n n n =-+=-+为插值节点，相应的Lagrange

形式的插值多项式为

),()(~

)(),(1

11

1

22r k r m m k n n r k u t f u l t l u t p ∑∑

+-=+-==

其中

1

1()m w k w m k w

w k

t t l t t t +=-≠-=

-∏

(k=m-1, m, m+1)

∏

+≠-=--=

1

1)(~

n r

w n w w

r w r y y y y u l (r=n-1, n, n+1)

并将i t 和j u 代入22(,)p t u ，便得到了数表,,

(,)

i j i j x y f x y 。

3.进行曲面拟和系数矩阵[]rs c *

=C ，1

1

()

()

T

T

T

--=C B B B U G G G

其中

001110

1011[()]1k

k r i k x x x x x x x ϕ⎡⎤

⎢⎥

⎢⎥==⎢⎥

⎢⎥⎢⎥⎣

⎦B

，001110

1011[()]1

k

k s j k y y y y G y y y ψ⎡⎤

⎢⎥⎢⎥==

⎢⎥⎢⎥⎢⎥

⎣⎦

[(,)]i j f x y =U

k 从0逐渐增大，直到7

10

σ

-≤，便得到了要求精度的系数rs c 。

4.由前面得到的函数关系，根据重新取值的x ，y 可以分别得到新的数表，比

较两组数据观察逼近效果

三、全部源程

#include "stdio.h"

#include "stdafx.h"

#include "math.h"

//高斯列主元法解线性方程组

void gauss(double a[4][4],double b[4],double x[4]) {

intk,i,j,ik;

double mik,temp,maxa,sum,tempb;

for (k=0;k<3;k++)

{

temp=fabs(a[k][k]);

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

if (fabs(a[i][k])>temp)

{

ik=i;

maxa=fabs(a[i][k]);

}

if(ik!=k)

{

for (j=k;j<4;j++)

{

temp=a[k][j];

a[k][j]=a[ik][j];

a[ik][j]=temp;

}

tempb=b[k];

b[k]=b[ik];

b[ik]=tempb;

}

for (i=k+1;i<4;i++)

{

mik=a[i][k]/a[k][k];

for (j=k+1;j<4;j++)

a[i][j]=a[i][j]-mik*a[k][j];

b[i]=b[i]-mik*b[k];

}

}

x[3]=b[3]/a[3][3];

for (k=2;k>=0;k--)

{