数值分析matlab程序

数值分析（matlab程序）

曹德欣曹璎珞

第一章绪论

数值稳定性程序，计算P11 试验题一积分

function try_stable

global n

global a

a=input('a=');

N = 20;

I0 = log(1+a)-log(a);

I = zeros(N,1);

I(1) = -a*I0+1;

for k = 2:N

I(k) = - a*I(k-1)+1/k;

end

II = zeros(N,1);

if a>=N/(N+1)

II(N) = (1+2*a)/(2*a*(a+1)*(N+1));

else

II(N) =(1/((a+1)*(N+1))+1/N)/2;

end

for k = N:-1:2

II(k-1) = ( - II(k) +1/k) / a;

end

III = zeros(N,1);

for k = 1:N

n = k;

III(k) = quadl(@f,0,1);

end

clc

fprintf('\n 算法1结果算法2结果精确值') for k = 1:N,

fprintf('\nI(%2.0f) %17.7f %17.7f %17.7f',k,I(k),II(k),III(k)) end

function y = f(x)

global n

global a

y = x.^n./(a+x);

return

第二章非线性方程求解

下面均以方程y=x^4+2*x^2-x-3为例：

1、二分法

function y=erfen(a,b,esp)

format long

if nargin<3 esp=1.0e-4;

end

if fun(a)*fun(b)<0

n=1;

c=(a+b)/2;

while c>esp

if fun(a)*fun(c)<0

b=c;

c=(a+b)/2;

elseif fun(c)*fun(b)<0

a=c;

c=(a+b)/2;

else y=c; esp=10000;

end

n=n+1;

end

y=c;

elseif fun(a)==0

y=a;

elseif fun(b)==0

y=b;

else disp('these,nay not be a root in the intercal')

end

n

function y=fun(x)

y=x^4+2*x^2-x-3;

2、牛顿法

function y=newton(x0)

x1=x0-fun(x0)/dfun(x0);

n=1;

while (abs(x1-x0)>=1.0e-4) & (n<=100000000)

x0=x1;

x1=x0-fun(x0)/dfun(x0);

n=n+1;

end

y=x1

n

function y=fun(x)

y=x^4+2*x^2-x-3; 3、割线法

function y=gexian(x0,x1)

x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %根据初始XO 和X1求X2 n=1;

while (abs(x1-x0)>=1.0e-4) & (n<=100000000) %判断两个条件截止 x0=x1; %将x1赋给x0 x1=x2; %将x2赋给x1 x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %迭代运算 n=n+1; end y=x2 n

function y=fun(x)

y=x^4+2*x^2-x-3;

第四章

题目：推导外推样条公式：???

?

??

??

? ??=???????? ?????????? ?

?--------1232123211

22

3

32

2~~~~

~22~

n n n n n n n n d d d d M M M M

δμ

λμλμλδ，并编写程序与Matlab 的Spline 函数结果进行对比，最后调用追赶法解方程组。 解：设()s x 的二阶导数值''()(0,1,

,)i i s x M i n ==，

因为()s x 在[]1,i i x x +上是三次多项式，所以''()i s x 在[]1,i i x x +上是一次多项式，且可表示为

''11()i i

i i

i i i

x x x x s x M M h h ++--=+ （1） 对式（1）进行两次积分可得：

33111()()()()()66k k k k k k k k k k k

M M

s x x x x x P x x Q x x h h +++=

-+-+-+- （2） 由式（2）根据()k k y s x =和11()k k y s x ++=可得系数为：

6k k k k k y m h P h =

-

，116

k k k k k y m h

Q h ++=- （3） 将式（3）代入式（2）可得()s x 为：

3311111()()()()()()()6666

k k k k k k k k k k k k k k k k k M M y m h y m h

s x x x x x x x x x h h h h +++++=

-+-+--+-- 分别对()k k s x 和1()k k s x -进行求导，并根据两者相等整理得：

112k k k k k k M M M d μλ-+++= （4）

其中，11k k k k h h h μ--=

+，1k k k k h h h λ-=+，1111

6

()k k k k k k k k k y y y y d h h h h +-----=-+

但是要解给定的方程组，还需要两个另外的条件，而外推样条插值的条件可通过下面推理得出：令式（1）中的x 等于1x ，可得：

32

112122()M M M h h h h h =

+- （5） 1212122

()()n n n n n n n n M M

M h h h h h -------=

+- （6） 令式（4）中2k =，然后将式（5）代入得：

1123222

(2)(1)h h

M M d h h +

+-= （7） 令式（4）中1k n =-，然后将式（6）代入得：

1112122

(2)(1)n n n n n n n h h

M M d h h -------+

+-= （8） 将方程组（4）和式（7）、式（8）联立展开即得题目所求。

按照推导出的外推样条插值公式编程可得如下M 文件（spline_wt.m ） function spline_wt X=[0 1 2 3 4 5];

Y=[0 0.5 2 1.5 3.5 1.9]; % 调用自编程序

S = spline_w(X,Y); %调用matlab 提供程序 S1=spline(X,Y);

fnplt(S,'r',2); % 作图 hold on

fnplt(S1,'b',1);

hold on

plot(X,Y,'or'); % 画上节点

title('红线为自编程序曲线，蓝线为自带程序曲线')

% pp的第j行表示第j个三次多项式的4系数并写出分段多项式

pp = S.coefs

P1 = poly2str(pp(1,:),'(x-0)')

P2 = poly2str(pp(2,:),'(x-1)')

P3 = poly2str(pp(3,:),'(x-2)')

ppp=S1.coefs

PP1 = poly2str(ppp(1,:),'(x-0)')

PP2 = poly2str(ppp(2,:),'(x-1)')

PP3 = poly2str(ppp(3,:),'(x-2)')

function sp = spline_w(X,Y)

n = length(X);

h = diff(X);

d = diff(Y)./h;

d1(2:n-1)=6*diff(d)./(h(1:n-2)+h(2:n-1));

mu(2:n-1)=h(1:n-2)./(h(1:n-2)+h(2:n-1));

mu(n-1)=1-h(n-2)/h(n-1);

la(2:n-1)=1-mu(2:n-1);

la(2)=1-h(1)/h(2);

% 计算三对角方程组

a = mu(3:n-1);

b = 2*ones(n-2,1);

b(1)=2+h(1)/h(2);

b(n-2)=2+h(n-2)/h(n-1);

c = la(2:n-2);

u(1:n-2) = d1(2:n-1);

% 调用追赶法解方程组

M(2:n-1) = tridiag(a,b,c,u); M(1)=(1+h(1)/h(2))*M(2)-M(3)*h(1)/h(2); M(n)=(1+h(n-1)/h(n-2))*M(n-1)-M(n-2)*h(n-1)/h(n-2);

% 下面计算分段多项式的四个系数

S=zeros(n-1,4);

for k=0:n-2

S(k+1,1)=(M(k+2)-M(k+1))/(6*h(k+1));S(k+1,2)=M(k+1)/2;

S(k+1,3)=d(k+1)-h(k+1)*(2*M(k+1)+M(k+2))/6;

S(k+1,4)=Y(k+1);

end

sp = mkpp(X,S);%转换成Matlab 格式

% -------- 求解三对角线性方程组的追赶法--------

function x=tridiag(a,b,c,d)

n = length(d);

x = zeros(n,1);

for k = 2:n

b(k) = b(k) - a(k-1)*c(k-1)/b(k-1);

d(k) = d(k) - a(k-1)*d(k-1)/b(k-1);

end

x(n) = d(n)/b(n);

for k = n-1:-1:1

x(k) = ( d(k)-c(k)*x(k+1) ) / b(k);

end

根据所编写的程序可得三次多项式及系数为（包括自编和提供）：pp =

-0.9511 3.3533 -1.9022 0

-0.9511 0.5000 1.9511 0.5000

1.7556 -

2.3533 0.0978 2.0000

-1.5711 2.9133 0.6578 1.5000

-1.5711 -1.8000 1.7711 3.5000

P1 =-0.95111 (x-0)^3 + 3.3533 (x-0)^2 - 1.9022 (x-0)

P2 =-0.95111 (x-1)^3 + 0.5 (x-1)^2 + 1.9511 (x-1) + 0.5

P3 =1.7556 (x-2)^3 - 2.3533 (x-2)^2 + 0.097778 (x-2) + 2

ppp =

-0.9511 3.3533 -1.9022 0

-0.9511 0.5000 1.9511 0.5000

1.7556 -

2.3533 0.0978 2.0000

-1.5711 2.9133 0.6578 1.5000

-1.5711 -1.8000 1.7711 3.5000

PP1 =-0.95111 (x-0)^3 + 3.3533 (x-0)^2 - 1.9022 (x-0)

PP2 =-0.95111 (x-1)^3 + 0.5 (x-1)^2 + 1.9511 (x-1) + 0.5

PP3 =1.7556 (x-2)^3 - 2.3533 (x-2)^2 + 0.097778 (x-2) + 2

根据所编写的程序可得生成的图形如图1所示

第五章

摘要：针对多项式检验的困难性，文中先将多项式线性化为多元线性问题。在此基础上，简要推导出多元线性拟合的数学模型及系数的最小二乘估计（以二元为例），同时为了检验建立模型的准确性，推导出检验模型的三个指标（误差差方和、相关检验和F 检验）计算公式，最后用某矿的煤样实例进行试验研究，得出各次多项式拟合的检验结果。

1 多元线性拟合的数学模型

设一个随机变量y 与m 个一般变量12,,

,m x x x 的内在联系是线性的，而且统计相关。

y 与这些(1,2,

,)i x i m =之间可用如下线性关系表示：

01122m m y b b x b x b x ε=+++

+ （1）

上式称为多元线性拟合的数学模型。式中01,,,m b b b 是(1)m +个待估参数，这些参数

被称为回归系数；ε是服从(0,)N σ的随机变量。

多项式拟合可以通过线性化来化为多元线性拟合，如多项式模型为：

2012m i i i m i i y b b x b x b x ε=++++ (1,2,

,)i n =

对上式线性化，令1i x x =，22i x x =，…，m m i x x =

即可得如式（1）所示的多元线性模型，因此我们只需讨论多元线性模型即可。为了方便我们采用二次多项式来推导，多次多项式可根据二次近似得到。

2 系数的最小二乘估计

设给定n 组实测数据12;,(1,2,

,)i i i y x x i n =，将其带入式（1）可得n 个观测方程：

01122i i i y b b x b x =++ (1,2,

,)i n = （2）

式中随机误差12,,

n εεε是n 个相互独立且服从(0,)N σ的随机变量。

令12n y y Y y ??????=??????，01m b b B b ??????=??????

，01

n εεεε??????=??????，111221221

211

1n n x x x x X x x ?????

?=?

??

???

则式（2）可写为：Y XB ε=+

设012,,b b b ∧

∧

∧

是采用最小二乘估计方法求得的012,,b b b 估值，则多元线性方程为：

01212y b b x b x ∧∧∧∧

=++，

则i y ∧

与实测值i y 之间的差值为：

01212()i i i i V y y y b b x b x ∧∧∧∧

=-=-++

用最小二乘估计，所有实验点上偏差平方和

2

1

1

()

n n

i i

i

i

i i VV y y ∧

===-∑∑最小的要求下，设

012[,,]B b b b ∧

∧

∧

∧

T =，则可得法方程用矩阵表示为：

N B X Y ∧

T =

于是，1

1

()B N X Y X X X Y ∧

-T

T

-T

== （3） 由（3）式解出012,,b b b ∧

∧

∧

，即可得多元线性方程模型。

3 方程模型的检验

在实际问题中，事先并不能断定随机变量y 与一般变量12,,

,m x x x 之间是否有线性关

系。用n 组实测数据按最小二乘的方法总能找出回归方程模型。但所得的方程模型是否具有实际意义，需要对所得到的方程模型进行统计检验。本文采用了误差平方和和两种统计方法进行检验。 3.1 误差平方和2Q

按下式可计算出总差方和：

2

1

()n

i Q y y -

==-∑，自由度为1f n =-

因为：

2

2

2

2

11111

()[()()]()()2()()n

n

n

n

n

i i i i i Q y y y y y y y y y y y y y y -

∧∧-

∧

∧-

∧∧-

======-=-+-=-+-+--∑∑∑∑∑

()()y y y y ∧

∧-

--=0

则Q 可分解为2

11

()

n

i Q y y ∧-

==

-∑和2

21

()n

i Q y y ∧

==

-∑，1

Q 称为回归差方和，反映自变量x 的

重要程度自由度为12f =；2Q 为误差差方和，反映试验误差及其他因素对试验结果的影响，自由度为2213f n n =--=-。 3.2 相关系数检验法

一元线性方程中相关系数定义为：xy

x y

σρσσ=，其中x σ，y σ，xy σ分别表示,x y 的均

方差和协方差。因为它们是未知的，故用子样均方差和协方差来代替。最后得相关系数ρ的

估值：()()

n

i

i

xy x y

x x y y S R S S -

-

--=

=

∑ 其中x S ，y S ，xy S 分别为子样的方差和协

方差，x S =

，y S =

，1

()()

n

i

i

i xy x x y y S n

--

=--=

∑

在多元线性方程中，定义复相关系数：R =

=1R →时，线性关系密切，而当0R →时，线性关系不密切，甚至不存在线性关系。可通过下式计算来得到R

值：

2

1

22

(1)

(1)

(1)

Q R m

m

F R Q R n m n m =

=?=

-----

3.3 F 检验

如果变量y 与12,,

,m x x x 之间无线性关系，则模型的一次项系数012,,b b b 均为零。所

以，要检验变量y 与变量12,,

,m x x x 之间是否有线性关系，就是要检验假设0H 是否成立 012:0H b b ==

从1Q ，2Q 得定义可看出两者是独立的，且222122

2

2

(1),

(),

(1)Q Q Q n m n m χχχσσσ---，则

在矩阵满秩和原假设0H 成立的条件下

12

(,1)(1)

Q m F F m n m Q n m =

----

在α水平下，若计算的F F α>计（查F 分布表），认为0H 不可信，方程模型效果显著；反之，则认为0H 可信，线性效果不显著，所求方程没有实际意义。

4 MATLAB 实现及分析

实验中取用某煤矿的18个煤样，所得的18组数据列于表1，分析容重和灰分之间的关系。

表1取样数据表

样品号

1

2

3

4

5

6

7

8

9

10

11

12

容重x 1.5 1.2 1.7 1.4 1.8 1.3 1.3 1.5 1.7 1.3 1.5 1.5 灰分y% 25 4 30 20 36 7 5

24

33

4

17

24

样品号 13

14

15

16

17

18

容重x 1.6 1.4 1.6 1.5 1.4 1.5 灰分y%

25

6

26

24

20

9

2，

表2 多项式拟合检验结果表

多项式次数 1

2

3 4 临界值 误差差方和 378.2118 373.5294 369.4278 358.566 相关检验 0.893 0.8944 0.8956 0.8988 0.47 F 检验

62.9611

29.977

18.9112

13.6677

4.49

拟合结果图形见图1

图1 多项式拟合效果图

从表2中可以看出，随着次数的增加，误差差方和逐渐减少，说明拟合越接近于实际曲线。这也可以通过相关检验的变化得以证明，次数越高，相关系数越接近于1且大于临界值，说明多项式线性化后的模型显著且逐渐接近于线性关系。而F检验的结果（方程显著且大于临界值）却说明当次数增加到一定程度时，拟合模型可能出现不显著，也即拟合模型不可用。

从上述分析可以看出，对于一个特定的问题，并不是建立模型的次数越多越接近于真实曲线。当多项式的次数超过F检验结果可用的情况下，高次数的多项式拟合反而成为一种错误的结果。因此，当我们在建立模型时，对建立的模型结果要进行检验，不能只根据某些实验数据拟合结果可用，就用它来进行下一步的预测研究。

附录：matlab源码

function poly_nihe

%多元线性拟合及检验

x=[1.5 1.2 1.7 1.4 1.8 1.3 1.3 1.5 1.7 1.3 1.5 1.5 1.6 1.4 1.6 1.5 1.4 1.5];

y=[25 4 30 20 36 7 5 24 33 4 17 24 25 6 26 24 20 9];

nn=length(x);

n=1;

p=polyfit(x,y,n);

xi=linspace(1.2,1.8,100);

z=polyval(p,xi);

subplot(2,2,1),plot(x,y,'o',xi,z)

title('一次多项式拟合')

ay=mean(y);

dy=y-ay;

q=0;

for i=1:nn

q=q+dy(i)^2;

end

ax=mean(x);

dx=x-ax;

qx=0;

for i=1:nn

qx=qx+dx(i)^2;

end

cxy=0;

for i=1:nn

cxy=cxy+dx(i)*dy(i);

end

r=cxy/sqrt(qx*q) %自由度为16及显著性水平为0.05查相关系数临界值表，得0.468 xp=[1.5 1.2 1.7 1.4 1.8 1.3 1.3 1.5 1.7 1.3 1.5 1.5 1.6 1.4 1.6 1.5 1.4 1.5];

yp1=polyval(p,xp);

v=y-yp1;

q2=0;

for i=1:nn

q2=q2+v(i)^2;

end

q2

q1=q-q2;

f=q1*(nn-2)/q2 %自由度为16及显著性水平为0.05查F分布表，得4.49

n=2;

p=polyfit(x,y,n);

xi=linspace(1.2,1.8,100);

z=polyval(p,xi);

subplot(2,2,2),plot(x,y,'o',xi,z)

title('二次多项式拟合')

%F检验

yp2=polyval(p,xp);

v=y-yp2;

q2=0;

for i=1:nn

q2=q2+v(i)^2;

end

q2

q1=q-q2;

f=q1*(nn-n-1)/(q2*n)

%R相关系数检验

r_2=q1/q;

ff=r_2*(nn-n-1)/((1-r_2)*n);

r=sqrt(n*f/(nn-n-1+n*f))n=3;

p=polyfit(x,y,n);

xi=linspace(1.2,1.8,100);

z=polyval(p,xi);

subplot(2,2,3),plot(x,y,'o',xi,z)

title('三次多项式拟合')

yp3=polyval(p,xp);

v=y-yp3;

q2=0;

for i=1:nn

q2=q2+v(i)^2;

end

q2

q1=q-q2;

f=q1*(nn-n-1)/(q2*n)

%R相关系数检验

r_2=q1/q;

ff=r_2*(nn-n-1)/((1-r_2)*n);

r=sqrt(n*f/(nn-n-1+n*f))n=4;

p=polyfit(x,y,n);

xi=linspace(1.2,1.8,100);

z=polyval(p,xi);

subplot(2,2,4),plot(x,y,'o',xi,z)

title('四次多项式拟合')

yp4=polyval(p,xp);

v=y-yp4;

q2=0;

for i=1:nn

q2=q2+v(i)^2;

end

q2

q1=q-q2;

f=q1*(nn-n-1)/(q2*n)

%R相关系数检验

r_2=q1/q;

ff=r_2*(nn-n-1)/((1-r_2)*n);

r=sqrt(n*f/(nn-n-1+n*f))

第一题：

一：运行图形和结果：

1969、1995和2000年的人口为：

左图的计算为三次多项式拟合的函数计算而得的结果：单位（亿）

1969年人口为：g1969_3 =8.10404781077523

1969年人口为：g1995_3 =11.0391780295176

1969年人口为：g2000_3 =10.7846676285844

偏差平方和为：sumerror3 =0.472697009288482

右图的计算为二次多项式拟合的函数计算而得的结果：单位（亿）

1969年人口为：g1969_2 =8.05176837343242

1969年人口为：g1995_2 =12.74409380623

1969年人口为：g2000_2 =13.8025210025012

偏差平方和为：sumerror2 = 0.709032276353744

二：源代码

（注：把t的值存放在‘t.txt‘中，p的值存放在‘p.txt’中，并且把它们和M文件放到同一个文件夹里）

clear

x=load('x.txt');

y=load('y.txt');

n=length(x);

c(1:n)=1;

A=[c',x',x.^2',x.^3'];

[Q,R]=qr(A,0);

a=R\(Q'*y');

a=fliplr(a');

A1=[c',x',x.^2'];

[Q,R]=qr(A1,0);

a1=R\(Q'*y');

a1=fliplr(a1');

t = 1949:2000; % 三次多项式拟合

g = polyval(a,t);

plot(x,y,'o',t,g,'m')

hold on

y1=a(1)*x.^3+a(2)*x.^2+a(3)*x+a(4);

g1969_3 = polyval(a,1969)

g1995_3= polyval(a,1995)

g2000_3 = polyval(a,2000)

error3=y-y1;

sumerror3=error3*error3'

g = polyval(a1,t); %二次多项式拟合

plot(x,y,'o',t,g,'r')

hold on

y2=a1(1)*x.^2+a1(2)*x+a1(3);

g1969_2 = polyval(a1,1969)

g1995_2= polyval(a1,1995)

g2000_2 = polyval(a1,2000)

error2=y-y2;

sumerror2=error2*error2'

第二题：

一：图形及结果：

运算结果如下：

以上五种拟合方法它们的均方误差如下：

sum1 为指数函数拟合 y=a*exp(x)+b 的均方误差。

sum2为三次多项式拟合的均方误差。

sum3为二次多项式拟合的均方误差。

sum4为y^2=a*x^2+b*x+c 拟合的均方误差。

sum 为指数 y=a*exp(b*x) 拟合的均方误差。

从均方误差可以看出来，三次多项式的拟合效果较好！

二：源代码：

（注：把t的值存放在‘t.txt‘中，p的值存放在‘p.txt’中，并且把它们和M文件放到同一个文件夹里）

clear

t=load('t.txt');

p=load('p.txt');

%------------------指数函数拟合 y=a*exp(x)+b

m=length(t);

c(1:m)=1;

A1=[c',exp(t)'];

[Q,R]=qr(A1,0);

a1=R\(Q'*p');

x=-69.7:100;

y=a1(1)+a1(2)*exp(x);

axis([-80,100,-500,4000]);

plot(x,y,'m');

hold on

y1=a1(1)+a1(2)*exp(t);

error1=p-y1;

sum1=sqrt((error1*error1')/m)

%---------------------三次多项式拟合

A2=[c',t',t.^2',t.^3'];

[Q,R]=qr(A2,0);

a2=R\(Q'*p');

a2=fliplr(a2');

x=-69.7:100;

y = polyval(a2,x);

axis([-80,100,-500,4000]);

plot(x,y,'r')

y1=a2(1)*t.^3+a2(2)*t.^2+a2(3)*t+a2(4);

error2=p-y1;

sum2=sqrt((error2*error2')/m)

%---------------------二次多项式拟合

A3=[c',t',t.^2'];

[Q,R]=qr(A3,0);

a3=R\(Q'*p');

a3=fliplr(a3');

x=-69.7:100;

y = polyval(a3,x);

axis([-80,100,-500,4000]);

plot(x,y,'b')

hold on

y1=a3(1)*t.^2+a3(2)*t+a3(3);

error3=p-y1;

sum3=sqrt((error3*error3')/m)

%-------------------------- y^2=a*x^2+b*x+c 拟合y=p.^2;

x=t;

c(1:m)=1;

A4=[c',x',x.^2'];

[Q,R]=qr(A4,0);

a4=R\(Q'*y');

x=-69.7:100;

y=sqrt(a4(3)*x.^2+a4(2)*x+a4(1));

axis([-80,100,-500,4000]);

plot(t,p,'o',x,y,'black');

y1=sqrt(a4(3)*t.^2+a4(2)*t+a4(1));

error4=y1-p;

sum4=sqrt((error4*error4')/m)

%------------------指数 y=a*exp(b*x) 拟合

y=log(p);

x=t;

A=[c',x'];

[Q,R]=qr(A,0);

a=R\(Q'*y'); x=-69.7:100; a(1)=exp(a(1)); y=a(1)*exp(a(2)*x);

axis([-80,100,-500,4000]); plot(t,p,'o',x,y,'g'); hold on

y1=a(1)*exp(a(2)*t); error=y1-p;

sum=sqrt((error*error')/m)

附加题第一题

二次曲面拟合地表（结合自己专业知识）

背景：

一般地表是连续的曲面组成的，我们可以根据一些地面已知点（知道三维坐标）用曲面来拟合地表。

步骤：

首先：我们可以到实地去取样点，就是所说的去实地测量地表点的（X ，Y ，Z ）坐标。在地表均匀地

选取一些点进行三维坐标（X ，Y ，Z ），然后我们可以根据所测量的点来拟合地表曲面。 其次：根据地表模型，一般用二次曲面来拟合地表，选择拟合公式。

Z=f （X ，Y ）=a 1 + a 2*X + a 3*Y + a 4*XY + a 5^X 2 + a 6*Y 2 再次：运用数值分析中的曲线拟合知识。

采样点：（x 1,y 1,z 1）(x 2,y 2,z 2)……(x m,y m ,z m )

基函数：)(,),(),(21x x x n ??? （n m ≥，一般n m >>）

函数族：)},(),()()(|)({2211y x a y x a y x a x s x s n n ???+++==Φ ，

记号：n m m n m m n n x x x x x x x x x ?????????????=)()()()()()()()()(212222111211????????? A ，????????????=n a a a 21a ，?????

???????=m y y y 21y 在这里：

采样点：我们已经实地测量出（X ，Y ，Z ）值。

基函数：由 Z=f （X ，Y ）=a 1 + a 2*X + a 3*Y + a 4*XY + a 5^X 2 + a 6*Y 2 我们6 个基函数：

)(1y x ，?=1，)(2y x ，?=x )(3y x ，? =y

)(4y x ，?=x*y )(5y x ，?=x 2 )(6y x ，?=y 2

我们可以求出 A

n m m m m m m

m y x y x y x y x y x y x y x y x y x ????

?????????=)()()()()()()()()(621226222121116112111，，，，，，，，，????????? A

????????????=62

1a a a a ， ?????

???????=m z z z z 21

我们用Cholesky 分解法求解，但直接计算A A T

是极其数值不稳定的，不宜采用。最常用的方法是

正交三角分解法（简称QR 分解），它间接地对A A T

进行了Cholesky 分解，其数值稳定性相当好。这

里我们借用Matlab 分解命令来求解。方法如下：

)0,(],[A qr R Q =

得到n n n m n m R Q A ???=，其中Q 的列是规范正交的，即n T I Q Q =，R 是非奇异上三角矩阵。 此时R R QR Q R A A T T T T ==（间接做了Cholesky 分解），法方程为

z Q Ra z Q R Ra R z A Aa A T T T T T T =?=?=

这样只需求解上面上三角方程组，可用

)(\z Q R a T =

这样就求出了二次曲面函数的系数。

最后：二次曲面函数绘图。调用MA TLAB 里的mesh （x,y,z ）绘出曲面。

绘出的图形如下：

注：已知的（x,y,z ）坐标是编的！没有进行实地测量！

程序代码如下：

clear

x=load('x1.txt'); y=load('y1.txt'); z=load('z1.txt'); m=length(x); c(1:m)=1;

%---------------------------用二次曲面拟合地表

%-------------z=a(1)+a(2)*x+a(3)*y + a(4)*x.*y+a(5)*x.^2+a(6)*y.^2 A=[c',x',y',(x.*y)',x.^2',y.^2'];

[Q,R]=qr(A,0); a=R\(Q'*z');

T = 1:40; S= 1:40;

[t,s]=meshgrid(T,S);

g=a(1)+a(2)*t+a(3)*s + a(4)*t.*s+a(5)*t.^2+a(6)*s.^2 ; mesh(t,s,g)

注：以上程序中

X1=[-6 1 4 2 6 9 12 15 16 10 7 3 8 12 15 11 17 20 14 21 25 24

29 32 28 35 38 30]

Y1=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

23 24 25 26 27 28]

Z1=[25 27 26 29 35 28 27 15 29 38 46 50 26 20 32 16 29 35 40

56 50 48 43 40 35 31 28 20]

运行代码之前，把编好的x 坐标放在‘x1.txt ‘，y 坐标放在‘y1.txt ‘，z 坐标放在‘z1.txt ‘并把它们和M 文件放在同一个文件夹下。然后就可以运行。

第六章

题目：求积分方程2

1

2

123)()(++-+

=

?

x x s x e e ds s y e x y 的数值解，并与精确解比较（精确解为x e x y =)(）

解：

1、理论分析

设],[b a 离散化为n x x x ,,,21 ，下面求)(i x y 的近似值i y 。

对方程中的积分用某个求积公式：

∑?

=≈n

j j j b

a

x F A dx x F 1

)()(近似代替得：

)()(),()(1

x f x y x x k A x y n

j j j j +≈∑=λ

记)(),,(),(i i j i ij i i x f f x x k k x y y ==≈得：

n i f y k A y i

n

j j ij j i ,,2,11

=+=∑=λ

这样转化为关于n y y y ,,,21 的线性方程组：

F Y K I =-)(λ

其中I 是单位矩阵，T

n y y y Y ],,,[21 =，T

n f f f F ],,,[21 =

?

??????

?????=nn n n n n n n n k A k A k A k A k A k A k A k A k A K

2

21

1222

221

11122111