若函数S（x）满足下列条件：
1）S（x）在每个区间[xi,xj]上是不高于3次的多项式。
2）S（x）及其2阶导数在[a,b]上连续。
则称S（x）使关于分划⊿的三次样条函数。
（2）matlab源程序：
clc,clear
x=0:
1:
20;
y=[9.018.967.967.978.029.0510.1311.1812.2613.2813.3212.6111.2910.229.157.97.958.869.8110.8010.93];
n=length（x）;
l
（1）=0;m（n）=0;
h=diff（x）;
df=diff（y）./diff（x）;
d
（1）=0;d（n）=0;
forj=2:
n-1
l（j）=h（j）/（h（j-1）+h（j））;
m（j）=h（j-1）/（h（j-1）+h（j））;
d（j）=6*（df（j）-df（j-1））/（h（j-1）+h（j））;
end
m=m（2:
end）;
u=diag（m,-1）;r=diag（l,1）;a=diag（2*ones（1,n））;
A=u+r+a;
M=inv（A）*d';
symsg
forj=1:
n-1
s（j）=M（j）*（x（j+1）-g）^3/（6*h（j））+M（j+1）*（（g-x（j））^3/（6*h（j）））+（y（j）-M（j）*h（j）^2/6）*（x（j+1）-g）/h（j）+（y（j+1）-M（j+1）*h（j）^2/6）*（g-x（j））/h（j）;
end
s
r=0;
forj=1:
n-1
df=diff（s（j）,g）;
warningoffall;
q=int（sqrt（1+df.^2）,g,j-1,j）;
r=r+q;
end
L=vpa（r,8）;
disp（'thelengthofthelabelisL='）;
disp（L）;
forj=1:
n-1
S（j,:
）=sym2poly（s（j））;
end
forj=1:
n-1
x1=x（j）:
0.1:
x（j+1）;
y1=polyval（S（j,:
）,x1）;
ifj==1
y2=y1;
else
fori=1:
11
k=（j-1）*10+i;
y2（k）=y1（i）;
end
end
end
x2=x
（1）:
0.1:
x（n）;
plot（x,y,'o'）
grid
holdon
plot（x2,y2,'r'）
（3）实验结果分析
图2.1铺设河底电缆长度
图2.2铺设河底光缆的曲线图
由三次样条插值得出的函数曲线的长度和即铺设河底电缆的长度为26.498514。
为了提高插值精度，用三次样条插值可以增加插值节点的办法来满足要求，而且在给定节点数的条件下，三次样条插值的精度要优于多项式插值以及线性分段插值，虽然舍弃了降低误差这个优点，但是其曲线的光滑性要好一些。
3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。
时刻
0
1
2
3
4
5
6
7
8
9
10
11
12
平均气温
15
14
14
14
14
15
16
18
20
20
23
25
28
时刻
13
14
15
16
17
18
19
20
21
22
23
24
平均气温
31
34
31
29
27
25
24
22
20
18
17
16
（1）解题思想和数学原理：
对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。
又因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。
最小二乘法，形成法方程是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。
对于给定的测量数据（xi,fi）（i=1,2,…，n），设函数分布为
特别的，取
为多项式
（j=0,1,…，m）
则根据最小二乘法原理，可以构造泛函
令
（k=0,1,…，m）
则可以得到法方程
求该解方程组，则可以得到解
，因此可得到数据的最小二乘解
（2）matlab源程序：
x=[0123456789101112131415161718192021222324];%给出题目数据（时间）
y=[15141414141516182020232528313231292725242220181716];%给出题目数据（温度）
plot（x,y,'m*'）%画出各个离散数据点
holdon
forn=2:
4;%2、3、4代表拟合函数最高项次数
alltemp=25;%alltemp代表数据点总共有25个
A=zeros（n+1,n+1）;%定义初始正规方程组的系数矩阵A
C=ones（n+1,1）;%定义初始正规方程组的系数向量C
D=zeros（n+1,1）;%定义初始正规方程组的向量D
fori=1:
n+1
forj=1:
n+1
fork1=1:
alltemp
A（i,j）=A（i,j）+（x（k1）.^（i-1+j-1））;
end
end
fork2=1:
alltemp
D（i,1）=D（i,1）+（x（k2）.^（i-1））.*（y（k2））;
end
end%以上为计算出正规方程组矩阵A、D的所有元素的程序
tol=1.0e-12;
maxit=1000;
C=bicg（A,D,tol,maxit）;%使用bicg迭代算出正规方程组的系数向量C
p=0;%误差分量
E=0;%误差总量
ifn==2
b=0:
24;
f=C
（1）+C
（2）.*b+C（3）.*（b.^2）;
p=y（b+1）-f;
forv=1:
25
E=E+（p（v））.^2;
end
plot（b,f,'r-'）
end%以上是对2阶拟合函数的图形处理与误差计算
ifn==3
b=0:
24;
f=C
（1）+C
（2）.*b+C（3）.*（b.^2）+C（4）.*（b.^3）;
p=y（b+1）-f;
forv=1:
25
E=E+（p（v））.^2;
end
plot（b,f,'g-'）
end%以上是对3阶拟合函数的图形处理与误差计算
ifn==4
b=0:
24;
f=C
（1）+C
（2）.*b+C（3）.*（b.^2）+C（4）.*（b.^3）+C（5）.*（b.^4）;
p=y（b+1）-f;
forv=1:
25
E=E+（p（v））.^2;
end
plot（b,f,'b-'）
end%以上是对4阶拟合函数的图形处理与误差计算
C,E
end
n=2;%重新对n赋值，进行指数函数拟合
A=zeros（n+1,n+1）;%重新对A矩阵赋初值
C=zeros（n+1,1）;%重新对C向量赋初值
D=zeros（n+1,1）;%重新对D向量赋初值
fori=1:
n+1
forj=1:
n+1
fork=1:
alltemp
A（i,j）=A（i,j）+（x（k）.^（i-1+j-1））;
end
end
forl=1:
alltemp
D（i,1）=D（i,1）+（（x（l）.^（i-1））.*（log（y（l））））;
end
end%计算出A矩阵、D向量各元素数值
C=bicg（A,D,tol,maxit）;%利用bicg迭代求解系数
b=0:
24;
p=0;
E=0;
f=exp（C
（1）+C
（2）.*b+C（3）.*（b.^2））;
p=y（b+1）-f;
forv=1:
25
E=E+（p（v））.^2;
end
plot（b,f,'c-'）%对指数拟合函数进行图形处理和误差计算
b=-C（3）;
c=C
（2）/（2*b）;
a=exp（b*（c^2）+C
（1））;%算出题设要求的指数拟合函数的各个系数
a,b,c,E
gridon
legend（'测量数据','二次函数','三次函数','四次函数','指数拟合','Location','NorthWest'）
holdoff%holdon与holdoff联合使用，表示将各个曲线画在同一个图中
图3.1二次曲线拟合系数与2范数误差
图3.2三次曲线拟合系数与2范数误差
图3.3四次曲线拟合系数与2范数误差
图3.4指数曲线拟合系数与2范数误差
图3.5数据原始点与拟合曲线对比图
（3）实验结果分析：
根据国家有关规定，用的是02、08、14、20时四个观测时次的数据做平均，最有代表性。
从图中可以看出并不是多项式次数越高越好，随着次数的增高，曲线所呈现出的给定点处和实际的吻合度越好，但对于其他地方的吻合度降低了。
4.设计算法，求出非线性方程
的所有根，并使误差不超过
。
解：
（1）解题思想和算法实现：
对于一个非线性方程的数值解法很多。
在此介绍两种最常见的方法：
二分法和Newton法。
首先要研究函数的形态，确定根的数量和大致区间的位置。
对于二分法，其数学实质就是说对于给定的待求解的方程f（x），其在[a,b]上连续，f（a）f（b）<0，且f（x）在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f（a）f（c）<0是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。
重复运行计算，直至满足精度为止。
这就是二分法的计算思想。
Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式
产生逼近解x*的迭代数列{xk}，这就是Newton法的思想。
当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。
另外，若将该迭代公式改进为
其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。
本题中用matlab画曲线
如下：
图4.1y=6*x.^5-45*x.^2+20的曲线
由曲线可以看出，该方程有三个实根，根所在区间依次为：
[-1,-0.5][0.5,1][1.5,2]
采用Newton迭代法，依据迭代式：
，k=0,1,2,…（8-1）
该方程的迭代式为：
（8-2）
分别选用迭代初值
、
、
进行迭代，分别求得满足精度的三个实根。
（2）matlab源程序：
clc;clear
x=-2:
0.1:
2;
y=6*x.^5-45*x.^2+20;
plot（x,y,'g'）%画相应的曲线
grid
title（'y=6*x.^5-45*x.^2+20曲线'）
formatlong
roots（[600-45020]）%根的真实解
clear
x0=input（'请输入迭代初值x0='）;
formatlong
i=1;
e=10^-4;%精度
x=x0-（6*x0.^5-45*x0.^2+20）/（30*x0.^4-90*x0）;
whileabs（x-x0）>e
i=i+1;
x0=x;
x=x0-（6*x0.^5-45*x0.^2+20）/（30*x0.^4-90*x0）;%迭代
end
display（'方程的根'）
x
display（'迭代的次数'）
i
（3）实验结果分析：
图4.2运行结果
对于Newton迭代法，三个初值x0都使得迭代收敛，这是非常重要的。
考虑Newton法迭代的收敛性条件：
（1）存在一个区间，满足。
由曲线和所选的三个区间可知这一条件满足。
（2）f（x）是[a,b]上的单调函数，即对一切不变号。
经验证所选的三个区间满足这一条件。
（3）f（x）的凹向在[a,b]上不变，即在[a,b]上不改变符号。
经验证所选的三个区间满足这一条件。
（4）
-1
1
1.5
>0
>0
>0
另外，可以看出Newton迭代法收敛速度也很快，且很快达到很高的精度，源于它一般是超线性收敛的。
5.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。
针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。
解:
（1）算法原理
由于一般线性方程在使用Gauss消去法求解时，从求解的过程中可以看到，若
=0，则必须进行行交换，才能使消去过程进行下去。
有的时候即使
0，但是其绝对值非常小，由于机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不正确。
因此有必要进行列主元技术，以最大可能的消除这种现象。
这一技术要寻找行r，使得
并将第r行和第k行的元素进行交换，以使得当前的
的数值比0要大的多。
这种列主元的消去法的主要步骤如下：
1.消元过程
对k=1,2,…,n-1,进行如下步骤。
1）选主元，记
很小，这说明方程的系数矩阵严重病态，给出警告，提示结果可能不对。
2）交换增广阵A的r，k两行的元素。
（j=k,…,n+1）
3）计算消元
（i=k+1,…,n;j=k+1,……,n+1）
2.回代过程
对k=n,n-1,…,1,进行如下计算
至此，完成了整个方程组的求解。
（2）matlab源程序：
%非压缩,dat51.dat、dat53.dat
clear;clc
fp=fopen（'dat53.dat','rb'）;
id=fread（fp,1,'int32'）;
ver=fread（fp,1,'int32'）;
N=fread（fp,1,'int32'）;
q=fread（fp,1,'int32'）;
p=fread（fp,1,'int32'）;
fori=1:
N
A（i,:
）=fread（fp,N,'float'）;
end
fori=1:
N
d（i）=fread（fp,1,'float'）;
end
%正向消去过程
fori=1:
N-q
fork=1:
p
ll=A（i+k,i）/A（i,i）;
forj=i:
i+q
A（i+k,j）=A（i+k,j）-ll*A（i,j）;
end
d（i+k）=d（i+k）-ll*d（i）;
end
end
fori=N-q+1:
N
fork=1:
N-i
ll=A（i+k,i）/A（i,i）;
forj=i:
N
A（i+k,j）=A（i+k,j）-ll*A（i,j）;
end
d（i+k）=d（i+k）-ll*d（i）;
end
end
%回代过程
x（N）=d（N）/A（N,N）;
fori=N-1:
-1:
1
S=0;
ifi+q>N
cv=N;%cv-criticalvalue
elsecv=i+q;
end
forj=i+1:
cv
S=S+A（i,j）*x（j）;
end
x（i）=（d（i）-S）/A（i,i）;
end
x
%压缩,dat52.dat、dat54.dat
clear;clc
fp=fopen（'dat54.dat','rb'）;
id=fread（fp,1,'int32'）;
ver=fread（fp,1,'int32'）;
N=fread（fp,1,'int32'）;
q=fread（fp,1,'int32'）;
p=fread（fp,1,'int32'）;
fori=1:
N
A（i,:
）=fread（fp,p+q+1,'float'）;
end
fori=1:
N
d（i）=fread（fp,1,'float'）;
end
%正向消去过程
fori=1:
N
ifi+pcv=p;%cv-criticalvalue
else
cv=N-i;
end
fork=1:
cv
ll=A（i+k,p+1-k）/A（i,p+1）;
forj=p+1:
p+q+1
A（i+k,j-k）=A（i+k,j-k）-ll*A（i,j）;
end
d（i+k）=d（i+k）-ll*d（i）;
end
end
%回代过程
x（N）=d（N）/A（N,p+1）;
fori=N-1:
-1:
1;
S=0;
ii=i+1;
jj=p+2;
while（ii<=N）&&（jj<=p+q+1）
S=S+A（i,jj）*x（ii）;
ii=ii+1;
jj=jj+1;
end
x（i）=（d（i）-S）/A（i,p+1）;
end
x
（3）实验结果：
非压缩矩阵求解结果（部分）
压缩矩阵求解结果（部分）
（4）分析心得：
采用Gauss消去法时，如果在消元时对角线上的元素始终较大，那么本方法不需要进行列主元计算，计算结果一般就可以达到要求，否则必须进行列主元这一步，以减少机器误差带来的影响，使方法得出的结果正确。
（5）实例
化学反应方程式严格遵守质量守恒定律，书写化学反应方程式写出反应物和生成物后，往往左右两边各原子数目不相等，不满足质量守恒定律，这就需要通过计算配平来解决。
对于化学反应方程式X1KMnO4+x2MnSO4+x3H2O=x4MnO2+x5K2SO4+x6H2SO4，可按以下方式配平：
上述化学反应式中包含5种不同的原子（钾、锰、氧、硫、氢），按照原子守恒有:
K:
x1=2x5
Mn:
x1+x2=x4
O:
4x1+4x2+x3=2x4+4x5+4x6
S:
x2=x5+x6
H:
3x2=2x6
上述方程组中有5个方程，6个未知量，因此有无数解。
可先令x6=1，于是形成方程组。
使用Gauss消去法求解此方程组得:
x=[1.0000,1.5000,1.0000,2.5000,0.5000,1.0000]；由于化学方程式通常取最简的正整数，因此将x乘2得配平的化学方程式：
2KMnO4+3MnSO4+2H2O=5MnO2+K2SO4+2H2SO4
程序如下：
clear;clc;
a=[1000-20;
110-100;
441-2-4-4;
0100-1-1;
00200-2;
00000-1];
b=[00000-1]';
n=length（b）;
fork=1:
n-1
ifa（k,k）==0
disp（'error'）;
return;
end
fori=k+1:
n
a（i,k）=a（i,k）/a（k,k）;
forj=k+1:
n
a（i,j）=a（i,j）-a（i,k）*a（k,j）;
end
b（i）=b（i）-a（i,k）*b（k）;
end
end
x（n）=b（n）/a（n,n）;
fork=n-1:
-1:
1
S=b（k）;
forj=k+1:
n
S=S-a（k,j）*x（j）;
end
x（k）=S/a（k,k）;
end
x
升级会员 