若函数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
11
12
平均气温
5
16
8
时刻
8
19
20
21
22
23
24
平均气温
3
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
升级会员 