uk=i;
end
end
ifabs（a（k,k））<=10^-6%判断a（k,k）是否接近0
error（'MATLAB:
Gausspp:
a（k,k）=0.SeeDoolittle.'）;
end
forj=1:
1:
n%交换列
ex=a（uk,j）;
a（uk,j）=a（k,j）;
a（k,j）=ex;
end
ex=b（uk）;
b（uk）=b（k）;
b（k）=ex;
fori=k+1:
1:
n%高斯消去过程
a（i,k）=a（i,k）/a（k,k）;
forj=k+1:
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
x（n）=b（n）/a（n,n）;%回代过程
fork=n-1:
-1:
1
S=b（k）;
forj=k+1:
1:
n
S=S-a（k,j）*x（j）;
end
x（k）=S/a（k,k）;
end
end
（2）分别读取文件fun001.dat至fun005.dat并采用列主元Gauss消去法求解的源代码,说明:
求解fun001.dat至fun005.dat，只需将第三行代码中的fun001.dat代换为fun002.dat至fun005.dat即可。
clc;
clear;
fid=fopen（'fun001.dat','rb'）;%打开fun00X.dat文件
[H,count]=fread（fid,3,'int'）;%读取前三个数据
ver=dec2hex（H
（2））;%判断版本
n=H（3）;%读取阶数
ifver=='201'%如果是压缩格式下的条状对角阵
[H,count]=fread（fid,2,'int'）;%读取上下带宽
p=H
（1）;%上带宽
q=H
（2）;%下带宽
A=[];%创建A矩阵
[H,count]=fread（fid,inf,'float'）;%读取剩余数据
t=1;%t代表读出的值的序号
fori=1:
n
forj=i:
i+p+q
A（i,j）=H（t）;
t=t+1;
end
end
A=A（1:
10,1+p:
n+p）;%求系数矩阵A
fori=1:
n
B（i）=H（t）;
t=t+1;
end%求B
else%如果不是压缩格式下的条状对角阵
[H,count]=fread（fid,inf,'float'）;%读取剩余数据
t=1;%t代表读出的值的序号
fori=1:
n
forj=1:
n
A（i,j）=H（t）;
t=t+1;
end
end%求系数矩阵A
fori=1:
n
B（i）=H（t）;
t=t+1;
end%求B
end
x=disanti_3_gauss（A,B）%调用函数求解
输出结果：
（1）fun001:
解为元素全为2.5的1024行列向量；
（2）fun002:
解为元素全为1.5的1560行列向量；
（3）fun003:
解为元素全为1.0的10行列向量；
（4）fun004:
解为元素全为1.895的10行列向量
（5）fun005:
解为元素全为1.895的1024行列向量；
4.已知某产品从1900年到2010年每隔10年的产量,用多项式插值和三次样条插值的方法，画出每隔一年的插值曲线的图形,试计算并比较在不同方法下的2005年以及2015年的产量。
年份
1900
1910
1920
1930
1940
1950
产量
75.995
91.972
105.711
123.203
131.699
150.697
年份
1960
1970
1980
1990
2000
2010
产量
179.323
203.212
226.505
251.525
291.854
325.433
算法思想：
多项式插值就是对于给定的平面上的n个数据点{（
）}（i=0,1,2,…,n），构造出一个函数,使得每个数据点都通过这条多项式曲线。
用这样的一条多项式曲线可以计算一些未测的中间值或者预测这些数据所在范围以外的值。
多项式插值函数的构造有拉格朗日形式和牛顿形式。
若所给定的n+1个数据点都不相同，则存在满足插值条件的唯一形式的次数不超过n的多项式P（x），因此通过拉格朗日形式和牛顿形式获得的多项式是相同的。
因为牛顿形式多项式插值具有拉格朗日形式所不具有的优点，即它可以用新获得的数据点来获得新的数据点，而不必重新计算，因此本题采用牛顿形式来求得插值多项式
三次样条插值也是构造多项式插值函数，不同之处在于三次样条插值是在每一个子区间[
，
]上构造一个三次多项式，并且这些多项式曲线在节点
处光滑地连接起来，这样的分段的三次多项式称之为三次样条函数。
三次样条插值先根据不同的边界条件确定μ、λ、d的值，再用《计算方法教程》P49页追赶法算法求解M的值，最后用P111页EVASPLINE算法和P102页FINDK算法求输入年份的产量并画图。
源程序：
Newton多项式插值
functionyk=disiti_4_newton（x,y,xk）
%本函数用Newton多项式插值计算yk
%x为插值点行向量，y为插值点函数值行向量
n=max（size（x））;
B=zeros（n,n+1）;%B是包含了节点的差商表
B（:
1）=x;
B（:
2）=y;
fori=3:
n+1%计算差商表B
forj=i-1:
n
B（j,i）=（B（j,i-1）-B（j-1,i-1））/（B（j,1）-B（j-i+2,1））;
end
end
yi=zeros（1,n）;
fori=1:
n%提取Newton插值系数
yi（i）=B（i,i+1）;
end
yk=yi（n）;
fori=n-1:
-1:
1%计算yk
yk=yk*（xk-x（i））+yi（i）;
end
%每隔1年计算函数值并画图
xm=1900:
1:
2010;
m=max（size（xm））;
ym=zeros（1,m）;
forj=1:
m
ym（j）=yi（n）;
fori=n-1:
-1:
1%计算yk
ym（j）=ym（j）*（xm（j）-x（i））+yi（i）;
end
end
plot（xm,ym）;%画图
end
求解输入为：
>>x=1900:
10:
2010;
y=[75.99591.972105.711123.203131.699150.697179.323203.212226.505251.525291.854325.433];
>>yk1=disiti_4_newton（x,y,2005）
>>yk2=disiti_4_newton（x,y,2015）
求解输出为：
Newton插值曲线图形：
由计算结果可知，2005年的产量为332.2477，高于2010年的产量，这个结果不太可靠，且2015年的产量为-17.8236，为负值，是错误的，所以此处采用Newton插值多项式拟合得到的曲线对于未来年份产量的预测不准确。
三次样条曲线插值
functionyk=disiti_4_sanciyangtiao（x,y,xk,flag,v1,vn）
%三次样条曲线插值求yk
%x为插值节点行向量，y为插值节点函数值行向量
%flag=1:
自然样条（端点二阶导数为0）；
%flag=2:
第一类边界条件（两端点一阶导数给定）；
%flag=3:
第二类边界条件；
%若flag=3,v1，vn为左，右端点导数值，若flag=1或2，v1=vn=0;
%初始化
n=max（size（x））;a=zeros（1,n）;b=zeros（1,n）;c=zeros（1,n）;h=zeros（1,n）;d=zeros（1,n）;
fori=2:
n-1%计算三阶差商求d
d（i）=（（y（i+1）-y（i））/（x（i+1）-x（i））-（y（i）-y（i-1））/（x（i）-x（i-1）））/（x（i+1）-x（i-1））;
end
d=6*d;
h
（2）=x
（2）-x
（1）;
fori=2:
1:
n-1
h（i+1）=x（i+1）-x（i）;
c（i）=h（i+1）/（h（i）+h（i+1））;%计算c
a（i）=1-c（i）;%计算a
b（i）=2;%对b赋值
end
b
（1）=2;b（n）=2;
ifflag==1%自然样条（两端点二阶导数为0）；
d
（1）=0;d（n）=0;c
（1）=0;a（n）=0;
end
ifflag==2%第一类边界条件（两端点一阶导数给定）；
c
（1）=1;a（n）=1;
d
（1）=6*（（y
（2）-y
（1））/（x
（2）-x
（1））-v1）/h
（2）;
d（n）=6*（vn-（y（n）-y（n-1））/（x（n）-x（n-1）））/h（n）;
end
ifflag==3%第二类边界条件；
c
（1）=-2;a（n）=-2;
d
（1）=-12*h
（2）*（（d（3）-d
（2））/（6*（x（4）-x
（1））））;
d（n）=-12*h（n）*（（d（n-1）-d（n-2））/（6*（x（n）-x（n-3））））;
end
%用追赶法计算M
u=zeros（1,n）;ym=zeros（1,n）;l=zeros（1,n）;M=zeros（1,n）;u
（1）=b
（1）;ym
（1）=d
（1）;
fork=2:
1:
n
l（k）=a（k）/u（k-1）;
u（k）=b（k）-l（k）*c（k-1）;
ym（k）=d（k）-l（k）*ym（k-1）;
end
M（n）=ym（n）/u（n）;
fork=n-1:
-1:
1
M（k）=（ym（k）-c（k）*M（k+1））/u（k）;
end
%查找xk所在区间
k=1;
fori=2:
n
ifxk<=x（i）
k=i;
break;
end
end
ifxk>=x（n）
k=n;
end
%计算yk
h0=x（k）-x（k-1）;
x1=x（k）-xk;
x2=xk-x（k-1）;
yk=（M（k-1）*x1^3/6+M（k）*x2^2/6+（y（k-1）-M（k-1）*h0^2/6）*x1+（y（k）-M（k）*h0^2/6）*x2）/h0;
end
求解输入为：
>>x=1900:
10:
2010;
y=[75.99591.972105.711123.203131.699150.697179.323203.212226.505251.525291.854325.433];
>>yk1=disiti_4_sanciyangtiao（x,y,2005,1,0,0）
>>yk2=disiti_4_sanciyangtiao（x,y,2015,1,0,0）
求解输出为：
由计算结果可知，2005年的产量为309.7236，低于2010年的产量，结果较为准确，且2015年的产量为641.1424，高于2010年的产量，符合预测，所以采用三次样条插值多项式拟合得到的曲线对于未来年份产量的预测较为准确。
由每隔1年计算函数值并画图源程序
xm=1900:
1:
2010;
m=max（size（xm））;
ym=zeros（1,m）;
fori=1:
m
ym（i）=disiti_4_sanciyangtiao（x,y,xm（i）,1,0,0）;
end
plot（xm,ym）;
三次样条曲线插值曲线图形
5.假定某天的气温变化记录如下表所示，试用最小二乘法找出这一天的气温变化的规律,试计算这一天的平均气温。
时刻
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
32
31
29
27
25
24
22
20
18
17
16
考虑使用二次函数、三次函数、四次函数以及指数型的函数
计算相应的系数，估算误差，并作图比较各种函数之间的区别。
算法思想：
在有的实际问题当中用插值方法并不合适，当数据点的数目很大时，要求p（x）通过所有的数据点，可能失去原数据所表示的规律。
如果数据点是由测量而来的，必然带来误差。
插值法要求准确通过这些不准确的数据点是不合适的。
在这种情况下，不用插值标准而用其他近似标准更合理，这就是最小二乘近似问题。
指数型函数拟合：
将指数型函数公式
进行变形，改写为
，将lnC当做因变量，t为自变量可以转化为利用多项式函数进行最小二乘拟合，再由二次函数拟合求出b,c,a的值。
源程序：
多项式拟合源程序：
function[xm,ave,p]=diwuti_5_duoxiangshinihe（x,y,n）
%多项式最小二乘法函数
%n为多项式系数个数，xm为计算出的多项式拟合系数，ave为平均温度，p为误差
m=max（size（x））;
g=zeros（m,n+1）;
omiga=zeros（m,1）;
xm=zeros（n,1）;
fori=1:
n%生成矩阵G
g（:
i）=x.^（i-1）;
end
g（:
n+1）=y;
g（1,1）=1;
fork=1:
n
sum0=0;%形成矩阵Qk
fori=k:
m
sum0=sum0+g（i,k）^2;
end
sigema=-sign（g（k,k））*sqrt（sum0）;
omiga（k）=g（k,k）-sigema;
forj=k+1:
m
omiga（j）=g（j,k）;
end
beta=sigema*omiga（k）;
g（k,k）=sigema;%变换G
forj=k+1:
n+1
sum1=0;
fori=k:
m
sum1=sum1+omiga（i）*g（i,j）;
end
t=sum1/beta;
fori=k:
m
g（i,j）=g（i,j）+t*omiga（i）;
end
end
end
xm（n）=g（n,n+1）/g（n,n）;%解三角方程
fori=n-1:
-1:
1
sum2=0;
forj=i+1:
n
sum2=sum2+g（i,j）*xm（j）;
end
xm（i）=（g（i,n+1）-sum2）/g（i,i）;
end
fori=1:
m%计算各时刻拟合温度
yk（i）=0;
forj=1:
n
yk（i）=yk（i）+xm（j）*x（i）^（j-1）;
end
end
ave=sum（yk）/m;%计算平均温度
plot（x,yk）%输出图像
p=0;
fori=n+1:
m%计算误差
p=p+g（i,n+1）^2;
end
end
输出结果：
二次多项式拟合结果
>>x=0:
24;
y=[15141414141516182020232528313231292725242220181716];
n=3;
>>[xm,ave,p]=diwuti_5_duoxiangshinihe（x,y,n）
输出为：
其中，xm=8.42742.5606-0.0920为拟合系数，由常数项到高阶；
ave=21.1200为平均温度；
p=253.6476为误差；
二次多项式拟合拟合曲线为：
三次多项式拟合
>>x=0:
24;
y=[15141414141516182020232528313231292725242220181716];
n=4;
>>[xm,ave,p]=diwuti_5_duoxiangshinihe（x,y,n）
输出为：
其中，xm=13.4072；-0.2164；0.2032；-0.0082为拟合系数，由常数项到高阶；
ave=21.1200为平均温度；
p=110.2954为误差；
三次多项式拟合拟合曲线为：
四次多项式拟合
>>x=0:
24;
y=[15141414141516182020232528313231292725242220181716];
n=5;
>>[xm,ave,p]=diwuti_5_duoxiangshinihe（x,y,n）
输出为：
其中，xm=16.6766；-3.5547；0.8592；-0.0513；0.0009为拟合系数，由常数项到高阶；
ave=21.1200为平均温度；
p=43.9298为误差；
四次多项式拟合拟合曲线为：
指数型函数拟合源代码
x=0:
24;
y=[15141414141516182020232528313231292725242220181716];
y=log（y）;%计算新的y值
n=3;%二次函数拟合
m=max（size（x））;
g=zeros（m,n+1）;
omiga=zeros（m,1）;
xm=zeros（n,1）;
fori=1:
n%生成矩阵G
g（:
i）=x.^（i-1）;
end
g（:
n+1）=y;
g（1,1）=1;
fork=1:
n
sum0=0;%形成矩阵Qk
fori=k:
m
sum0=sum0+g（i,k）^2;
end
sigema=-sign（g（k,k））*sqrt（sum0）;
omiga（k）=g（k,k）-sigema;
forj=k+1:
m
omiga（j）=g（j,k）;
end
beta=sigema*omiga（k）;
g（k,k）=sigema;%变换G
forj=k+1:
n+1
sum1=0;
fori=k:
m
sum1=sum1+omiga（i）*g（i,j）;
end
t=sum1/beta;
fori=k:
m
g（i,j）=g（i,j）+t*omiga（i）;
end
end
end
xm（n
升级会员 