数值分析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 ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=621a 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 文件放在同一个文件夹下。

然后就可以运行。

第六章
题目：求积分方程21
02
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 221
1222221
11122111
2、实验结果及分析
采用五等分和七等分的高斯—柯特斯公式得到的结果如图1
图1 五、七等分的高斯—柯特斯公式积分结果图
五等分的高斯—柯特斯公式运行结果和精确值分别为：
0.99997600726849 1.22137345337173 1.49178890469185 1.82207508278336 2.22548753168652 2.71821660945298
精确值为：
1.00000000000000 1.22140275816017 1.49182469764127 1.82211880039051
2.22554092849247
2.71828182845905
七等分的高斯—柯特斯公式运行结果和精确值分别为：
0.99999987238056 1.15356484767778 1.33071202762260 1.53506281335132
1.77079472644729
2.04272680957445 2.35641814165885 2.71828148155343
精确值为：
1.00000000000000 1.15356499489511 1.33071219744735 1.53506300925521
1.77079495243515
2.04272707026614 2.35641844238366 2.71828182845905
八等分的高斯—柯特斯公式运行结果和精确值分别为：
0.99999999923687 1.13314845220208 1.28402541570786 1.45499141350785 1.64872126944194
1.86824595600650
2.11700001499712 2.39887529213644 2.71828182638464
精确值为：
1.00000000000000 1.13314845306683 1.28402541668774 1.45499141461820 1.64872127070013 1.86824595743222
2.11700001661267 2.39887529396710 2.71828182845905
从图1中可以看出，不管是五等分公式还是七等分公式都具有良好的精度，画出的曲线图和精确值
完全重合。

但是从得出的数据结果来看，七等分的精度要好于五等分的精度，八等分的精度要好于七等分。

虽然八等分的高斯—柯特斯公式稳定性不是很好，但是对于计算本题确是稳定的。

因此，在使用八等分公式之前要进行分析。

如果可用，则具有较好的精度；否则在一般情况下不要使用。

3、程序代码：
五等分高斯—柯特斯公式：
function temp
clc
format long
I=eye(6);
A=[19/288 75/288 50/288 50/288 75/288 19/288];
x = linspace(0,1,6);
for i=1 : 6
for j= 1:6
k(i,j)=exp(x(i)+x(j));
end
end
for j=1 :6
for i=1:6
K(i,j)=k(i,j)*A(j);
end
end
F=f(x)';
Y=(I-K)\F
yyy=exp(x)
x1 =0:0.01:1;
yy=exp(x1);
plot(x,Y,'o',x1,yy,'b')
function s=f(x)
for i=1:6
s(i)=1.5*exp(x(i))-0.5*exp(x(i)+2);
end
七等分高斯—柯特斯公式：
function temp
clc
format long
I=eye(8);
A=[751/17280 3577/17280 1323/17280 2989/17280 2989/17280 1323/17280 3577/17280 751/17280];
x = linspace(0,1,8);
for i=1 : 8
for j= 1:8
k(i,j)=exp(x(i)+x(j));
end
end
for j=1 :8
for i=1:8
K(i,j)=k(i,j)*A(j);
end
end
F=f(x)';
Y=(I-K)\F
yyy=exp(x)
x1 =0:0.01:1;
yy=exp(x1);
plot(x,Y,'o',x1,yy,'b')
function s=f(x)
for i=1:8
s(i)=1.5*exp(x(i))-0.5*exp(x(i)+2);
end
八等分高斯—柯特斯公式：
function temp
clc
format long
I=eye(9);
A=[989/28350 5888/28350 -928/28350 10496/28350 -4540/28350 10496/28350 -928/28350 5888/28350 989/28350];
x = linspace(0,1,9);
for i=1 : 9
for j= 1:9
k(i,j)=exp(x(i)+x(j));
end
end
for j=1 :9
for i=1:9
K(i,j)=k(i,j)*A(j);
end
end
F=f(x)';
Y=(I-K)\F
yyy=exp(x)
x1 =0:0.01:1;
yy=exp(x1);
plot(x,Y,'o',x1,yy,'b')
function s=f(x)
for i=1:9
s(i)=1.5*exp(x(i))-0.5*exp(x(i)+2);
end
%----------------------------------------5点高斯-勒让德公式
A1=[0.2369268851 0.4786286705 0.5688888889 0.4786286705 0.2369268851];
I1=eye(5);
t=[-0.9061798459 -0.5384693101 0.0000000000 0.5384693101 0.9061798459];
x=0.5.*t+0.5;
for i=1 : 5
for j= 1:5
k1(i,j)=exp(x(i)+0.5.*t(j)+0.5);
end
end
for j=1 :5
for i=1:5
K1(i,j)=k1(i,j)*A1(j);
end
end
F1=f1(x);
F1=F1';
Y1=(I1-0.5.*K1)\F1
yyyy=exp(x)
x1 =0:0.01:1;
yy=exp(x1);
plot(x,Y1,'mo',x1,yy,'g')
function s1=f1(x)
for i=1:5
s1(i)=1.5*exp(x(i))-0.5*exp(x(i)+2);
end
%----------------------------------------复化的辛甫生公式
I=eye(10);
x=linspace(0,1,11);
for i=1 :10
for j=1:10
K(i,j)=exp(x(i)+x(j))+4*exp((x(i)+x(i+1))/2+(x(j)+x(j+1)) /2)+exp(x(i+1)+x(j+1));
end
end
F2=f2(x);
Y2=(I-(1/60)*K)\F2'
x2=x(1:10);
yyy=exp(x2)
plot(x2,Y2,'*',x2,exp(x2),'g')
function s2=f2(x)
for i=1:10
s2(i)=1.5*exp(x(i))-0.5*exp(x(i)+2);
end
老师提供例子
1.1
% exp1_1.m --- 算法的数值稳定性实验
% 见P1 例1 求积分
% I(n) = exp(-1)*int( x^n*exp(x),0,1 )
% 参见P2表1-1 和P5表1-2
function try_stable
global n % 定义全局变量n ,见后自定义函数f 中的参量
N = 20; % 计算N 个值
%--------------------------------------------
% 【算法1】用递推公式（P2 (1-2)式）
% I(k) = 1 - k*I(k-1)
% 取初值I0=1-exp(-1) (约有15位有效数字)
% [注] 如果取I0 = 0.6321,结果同书P2 表1-1,这里初值更精确
I0 = 1-exp(-1); % 初值
I = zeros(N,1); % 创建N x 1 零矩阵(列向量),用于储存递推的值
I(1) = 1-I0;
for k = 2:N
I(k) = 1 - k*I(k-1);
end
%--------------------------------------------
% 【算法1】用递推公式( P5 (1-6)式)
% II(k-1) = ( 1 - II(k) ) / k
% 取初值II(N) = 0
II = zeros(N,1);
II(N) = 0;
for k = N:-1:2
II(k-1) = ( 1 - II(k) ) / k;
end
%--------------------------------------------
% 调用matlab 高精度数值积分命令quadl 计算以便比较
% [注] 该命令以后学习,你现在要模仿使用
III = zeros(N,1);
e_1 = exp(-1);
for k = 1:N
n = k; % 给函数f 中的参量n 赋值
III(k) = e_1 * quadl(@f,0,1); % 求函数f 在[0,1]上的定积分( P1 (1-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 % 参量n 为全局变量
y = x.^n.*exp(x); % ★注意:这里一定要'点' 运算
return
%--------------------------------------------
1.2
% exp1_2.m --- 多项式求值的Horner 算法
% [简介] matlab中用多项式的系数(按降幂顺序)组成的向量来表示多项式,
% 以5 次多项式为例
% y = c1*x^5 + c2*x^4 + c3*x^3 + c4*x^2 + c5*x + c6
% 用向量
% P = [c1,c2,...,c6]
% 来表示,注意它是6 维向量(可以是行向量,也可以是列向量).
% [方法] Horner 算法(参见P7)
% y = c6+(c5+(c4+(c3+(c2+c1*x)*x)*x)*x)*x
% 计算量为5 次乘法,5 次加法,比按自然顺序计算大大减少了计算量
% [调用] y = polyval(P,x) 这里x 是矩阵(向量,标量),y 与x 同维数
function try_horner_method
P = [2 1 0 4 -5 6]; % 表示多项式
% P(x) = 2 x^5 + x^4 + 4 x^2 - 5 x + 6
x = [1,2,3]; % 自变量取值
clc, disp('调用自编程序计算结果:')
y1 = mypolyval(P,x) % 调用自编的Horner 算法程序
% 求y1 = [P(x1),P(x2),P(x3)] 的值
disp('调用matlab命令计算结果:')
y2 = polyval(P,x) % 调用matlab 中多项式求值命令
% 该命令也是用Horner 算法编的
% ------- 多项式求值的Horner 嵌套算法-------
function y = mypolyval(P,x)
% y = mypolyval(P,x) --- 多项式求值的Horner 嵌套算法
% P --- 向量(表示多项式)
% x --- 矩阵或向量或标量
% y = P(x)(维数同x)
np = length(P); % 向量P 的维数
[m,n] = size(x); % 矩阵x 的维数( m 是行数,n 是列数)
y(1:m,1:n) = P(1); % 产生矩阵y 它与x 同维数, 且每个元素都用P(1) 赋值
% 以上两句可合并写为y(size(x))=P(1);
for k = 2:np
y = y.*x + P(k); % ★注意: 这里是点乘,另外矩阵加一个数等于矩阵每个元素加这个数
end
% ----------------------------------------
2.1
% exp2_1.m --- 学习非线性方程f(x) = 0 数值求根命令fzero
% [调用] X = fzero(FUN,X0), FUN 是函数,求X0 附近的根
% 或X = fzero(FUN,[a b]), 求[a,b] 之间的根,这时FUN(a) 与FUN(b) 符号必须相反
% [注] 该命令集二分法, 割线法, 逆二次内插法于一体,
% 该命令只适用于奇重根,即在根的附近左右函数值反号
% example1:求sin(x)在3 附近的根
X1 = fzero(@sin,3)
% example2: 求f1 在[0,1]之间的根, 其中f1 是用内联方式定义的函数
f1 = inline('cos(x)-x');
X2 = fzero(f1,[0 1])
% example3: 求f2 = x-exp(-x) = 0 的根
%第一步: 通过作图估计根的大致位置(近似值)
f2 = inline('x-exp(-x)'); % 定义函数
x = linspace(0,1,101); % 采样: 把[0,1] 100 等分,x 为101 维的行向量
y = f2(x); % 计算f2 的函数值,其维数同x
plot(x,y,'r') % 作图(就是把各个点连起来),颜色红色
grid on % 显示坐标刻度网格
title('f(x) = x - e^{-x} = 0 的根') % 加标题
text(0.3,0,'根的大致位置\rightarrow') % 加文字说明
% 通过作图知道在[0.5,0.6]之间有一个根
%第二步: 调用fzero 命令求更精确的近似值,参见P23例6
X3 = fzero(f2,[0.5 0.6])
2.2
% exp2_2.m --- 学习非线性方程组F(X) = 0 数值求根命令fsolve
% [调用方法] X = fsolve(F,X0)
% 这里F 是向量函数(就是方程组),X0 是迭代初值(向量)
% [注] 该命令采用优化方法求解,默认的是一维搜索法
function study_fsolve
% 下面仅以一个例子说明其使用方法.求下面方程组的九个根
% 7x^3-10x-1-y = 0
% 8y^3-11y-1+x = 0
% 首先学习隐函数作图，通过作图观察根的大致位置(两个曲线的交点)
clf % 清图像
ezplot('7*x^3 - 10*x - 1-y',[-4,6,-6,4]); %第一个方程所确定的隐函数
% 后面的区间是x,y的范围，一开始可能取不到合适的，偿试几次就可以了。

hold on % 锁定前面的图，让后面的图叠加上去，否则后面的图就会把前面的冲掉。

ezplot('8*y^3 - 11*y - 1 + x',[-4,6,-6,4]);
grid on
% 通过作图发现在(0,-1) 附近有一个根(交点),调用fsolve 求更精确的解
clc
X0 = [0,-1]; % 初值
X = fsolve(@myfun,X0)
% ------------ 定义函数组myfun ----------------
function F = myfun(X)
F = [7*X(1)^3-10*X(1)-1-X(2);
8*X(2)^3-11*X(2)-1+X(1)];
% [注] 这里F 必须是列向量;
% ----------------------------------------------
2.3
% exp2_3.m --- (1)学习多项式求根命令roots
% (2)了解多项式求根的病态性
% [简介] roots 命令是把多项式求根转化为矩阵的特征值问题
% 参见P209 实验课题(六) 有详细说明
% 由于多项式求根大部分都是病态问题,此方法对20 次以下多项式效果较好% 例1:
% 求f(x) = 2 x^3 + x^2 + 1 的三个根
clc
C = [2,1,0,1]; % 表示上面多项式
F = roots(C)
%例2：代数方程求根的病态性实验( 见P33 例13 )
C = [1 2 3 4 5 6 7];
P = poly(C); % 以1,2,...,7为根的首一多项式,它和roots 是互为逆命令
% 作图( 即P33 图2-15(a) )
figure(1)
x = 0 : 0.1 : 8; y = polyval(P,x);
plot(x,y,'b',[0,8],[0,0],'r')
axis([0 8 -300 300])
%再求根,当然应该是1,2,...,7
F1=roots(P)
% 把系数P(2) 给一个小扰动再作图( 即P33 图2-15(b) )
figure(2)
P(2) = P(2) - 0.002;
x = 0 : 0.1 : 8; y = polyval(P,x);
plot(x,y,'b',[0,8],[0,0],'r')
axis([0 8 -300 300])
%再求根,看看根的变化是否剧烈
F2=roots(P)
3.1
% exp3_1.m --- 解线性方程组左除命令‘\’的学习
% ---------- example1 -----------
% Ax = b (x,b是列向量),当A是可逆矩阵时x = A\b 产生该方程组的解
% 其算法基于LU分解相当于列主元Gauss消去法
A = [ 1 5 -9
0 6 4
1 1 1];
b = [-16 24 6]';
x1 = A\b
% [注] 该命令适用求解中小型稠密线性方程,而且性态是较好的(非病态),是最常用的命令% 对于大型矩阵或病态的还有其它一些命令pcg,gmres,qmr等
% ---------- example2 -----------
% Ax = b ,当A是列满秩矩阵时x = A\b 产生该方程组唯一的最小二乘解
% 其算法基于解法方程组A'*A x = A'*b (见P125-126, 例11)
A = [2 1 1
1 -
2 -3
3 -
4 1
1 3 -2];
b = [-4 5 3 -2]';
x2 = A\b
% ----------- example3 ----------
% 解矩阵方程AX = B ,当A可逆或列满秩
A = [ 1 3
1 4];
B = [ 1 2
3 4];
X = A\B
% 自学与矩阵计算有关的一些常用函数
% det(A) --- 求方阵A的行列式
% inv(A) --- 求方阵的逆
% rank(A) --- 求矩阵A的(数值)秩
% rref(A) --- 化矩阵A为最简阶梯形
% norm(A,p) --- 求矩阵A的p-范数
% cond(A,p) --- 求矩阵A的p-范数的条件数
% eig(A) --- 求方阵A的特征值
3.2
% exp3_2.m --- 矩阵分解命令的学习
% ---------- LU分解----------
% [简介] 设A是非奇异矩阵,则有如下列主元的LU分解
% PA = LU
% 其中P是置换矩阵(又称排列矩阵),L是单位下三角矩阵且|l(ij)|<=1,U是上三角矩阵% 矩阵A如果有了上述分解则求解线性方程组Ax = b 就等价于分别求解两个三角方程组
% Ly = Pb 和Ux=y
% 而求解三角方程组是非常容易的.用这种方法求解方程组与列主元的Gauss消去法是等价的
% 但LU 分解还有其它优点. 参见P65 例9 和P70 实验课题(一)
A = [10 -7 0
-3 2 6
5 -1 5];
[L U P] = lu(A)
% 如果再给b,则下面是解两个三角方程组的程序. 参见P49(3-15)式
b = [7 4 6]';
[n n]=size(A);
x = zeros(n,1);
% 前代求解:Ly = Pb(解用x储存)
b = P*b;
x(1) = b(1);
for k = 2:n
x(k) = b(k)-L(k,1:k-1)*x(1:k-1);
end
% 回代求解:Ux = y (这里先y=x,解仍用x储存)
x(n) = x(n)/U(n,n);
for k = n-1:-1:1
x(k) = ( x(k)-U(k,k+1:n)*x(k+1:n) ) / U(k,k);
end
x
% ---------- Cholesky分解----------
% [简介] 设A是(实对称)正定矩阵,则有如下唯一的分解
% A = R'*R
% 其中R 是上三角矩阵且主对角元全为正
% 例
A = [4 -1 1
-1 4.25 2.75
1 2.75 3.5];
R = chol(A)
% [注1] 有了Chol分解,求解Ax = b 又可转化为分别求解两个三角方程组% 这是容易的,请同学们自己完成(常称为平方根法)
% [注2] P51 的chol 分解A = LDL'
% 避免了前面分解中的开方运算(开方运算比四则运算速度是慢的) 3.3
% exp3_3.m Gauss列主元消去法(示范程序)
function try_gauss
A = [ 1 5 -9
0 6 4
1 1 1];
b=[-3;10;3];
x=gauss(A,b) % 调用下面函数
function x = gauss(A,b)
% 输入: Ax = b 的系数矩阵A (可逆)和右端项b
% 输出:方程组的解x
[n,n] = size(A);
x = zeros(n,1);
Aug = [A,b]; % 增广矩阵
for k = 1:n-1
[piv,r] = max(abs(Aug(k:n,k))); % 找列主元所在子矩阵的行r
r = r + k - 1; % 列主元所在大矩阵的行
if r>k
Aug([k,r],:) = Aug([r,k],:); % 对增广矩阵实施行交换
end。