若函数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