计算方法上机题Word文件下载.docx
- 文档编号:8572883
- 上传时间:2023-05-11
- 格式:DOCX
- 页数:15
- 大小:161.74KB
计算方法上机题Word文件下载.docx
《计算方法上机题Word文件下载.docx》由会员分享,可在线阅读,更多相关《计算方法上机题Word文件下载.docx(15页珍藏版)》请在冰点文库上搜索。
i-1
M(k)=L(i,k)*U(k,j);
end
U(i,j)=A(i,j)-sum(M);
%给上三角阵U的第i行赋值
forj=i+1:
M(k)=L(j,k)*U(k,i);
L(j,i)=(A(j,i)-sum(M))/U(i,i);
%给单位下三角阵L的第i列赋值
%分别输出L,U的矩阵
L
U
3.编写程序实现大规模方程组的列主元高斯消去法程序,并对所附的方程组进行求解。
具体的要求参见所附的相关文档。
针对本专业中所碰到的实际问题,提炼一个使用方程组进行求解的例子,并对求解过程进行分析、求解。
fid=fopen('
fun001.dat'
'
r'
);
%获取数据文件文件头
headfile=fread(fid,3,'
uint32'
%确认数据文件标识号是否正确
id=headfile
(1);
ifid~=hex2dec('
F1E1D1A0'
)
fprintf('
不是正确的数据文件.\n'
return;
%获取数据文件版本号
ver=headfile
(2);
ifver==hex2dec('
101'
%获取方程阶数
n=headfile(3)
%获取稀疏矩阵
A=fread(fid,[n,n],'
float32'
%获取右端系数
B=fread(fid,[n,1],'
%对系数矩阵A和右端向量的方程组使用高斯方法求解
x=GAUSSPP(A,B,n);
ifver==hex2dec('
102'
%系数矩阵为条状对角阵
%获取方程阶数
n=headfile(3)
201'
%系数矩阵为压缩格式下的条状对角阵
n=headfile(3)
%获取条状矩阵的上下带宽
headfile2=fread(fid,2,'
p=headfile2
(1);
q=headfile2
(2);
%压缩系数矩阵
A=fread(fid,[p+q+1,n],'
x=GAUSSPP(A,B,n,p,q);
fclose(fid);
GAUSSPP(A,B,n)的程序如下:
functionx=GAUSSPP(A,B,n)
fork=1:
n-1
fori=k+1:
l(i,k)=A(i,k)/A(k,k);
forj=k:
A(i,j)=A(i,j)-l(i,k)*A(k,j);
B(i)=B(i)-l(i,k)*B(k);
x(n)=B(n)/A(n,n);
fork=n-1:
-1:
1
s=0;
forj=k+1:
s=s+A(k,j)*x(j);
x(k)=(B(k)-s)/A(k,k);
GAUSSPP(A,B,n,p,q)的程序如下:
functionx=GAUSSPP(A,B,n,p,q)
n-1
min(k+p,n)
l(i,p-i+k+1)=A(i,p-i+k+1)/A(k,p+1);
forj=0:
q
A(i,p-i+k+1+j)=A(i,p-i+k+1+j)-l(i,p-i+k+1)*A(k,p+1+j);
B(i)=B(i)-l(i,p-i+k+1)*B(k);
x(n)=B(n)/A(n,p+1);
forj=1:
min(q,n-k)
s=s+A(k,1+p+j)*x(j+k);
x(k)=(B(k)-s)/A(k,p+1);
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
1.多项式插值
利用多项式插值进行数据近似时要选取一个函数gk(x),通常选取如下的函数:
因此,欲寻求的函数就是n次多项式
根据差值条件可以得到下式:
我们可以把这个方程抽象成一个线性方程组:
这样就把求解数据近似的问题转化成了求解线性方程组的问题。
同时观察上式可以看出方程的系数是一个Vandermonde矩阵。
%多项式插值近似%
Y=1900:
10:
2010;
A=zeros(12);
YB=1900;
%根据数据构建11阶Vandermonde矩阵A
fori=1:
12;
forj=1:
A(i,j)=YB^(j-1);
YB=YB+10;
%根据数据构建矩阵B
B=[75.995;
91.972;
105.711;
123.203;
131.699;
150.697;
179.323;
203.212;
226.505;
251.525;
291.854;
325.433];
%调用列主元Gauss消去法函数,求多项式系数矩阵
x=Gauss_pivot(A,B);
a=zeros(1,12);
12
a(1,n)=x(13-n,1);
e1=polyval(a,[2005,2015])
T=1900:
plot(Y,B,'
)%红色表示原数据图像
holdon
plot(T,polyval(a,T),'
g'
)%绿色表示多项式插值近似图像
最终得到的结果是
e1=
342456
即2005年的产量为342,2015年的产量为456。
多项式插值的结果如图
(1)所示:
图
(1)多项式插值的图形
图中红色曲线表示原始数据,绿色曲线表示多项式插值的数据。
2.三次样条函数插值
由于分段多项式插值不够光滑,因此需要较高次的分段多项式,常用三项多项式,即在每一个子区间上市一个三次多项式,这样的分段三次多项式就是三次样条函数。
具体算法:
n=12;
%构造数据矩阵Y
%构造数据矩阵B
%构造步长矩阵h
h=zeros(1,n-1);
fork=1:
11
h(k)=Y(k+1)-Y(k);
%构造M的系数矩阵A
A=zeros(n);
%由第三类边界条件可得
A(1,2)=-2;
A(n,n-1)=-2;
fori=1:
A(i,i)=2;
fori=2:
A(i,i+1)=h(i)/(h(i-1)+h(i));
A(i,i-1)=h(i-1)/(h(i-1)+h(i));
%构造M的结果矩阵d
d=zeros(n,1);
d
(1)=-12*h
(1)*ChaShang_4(Y,B,1,2,3,4);
d(n)=12*h(n-1)*ChaShang_4(Y,B,n-3,n-2,n-1,n);
forj=2:
d(j)=6*ChaShang_3(Y,B,j-1,j,j+1);
M=Gauss_pivot(A,d);
%构造在不同区间内三次样条拟合结果的矩阵S
S=zeros(1,(n-1)*10+1);
num=1;
forx=1:
10
S(1,(i-1)*10+num)=
M(i)*((Y(i+1)-T((i-1)*10+num))^3)/(6*h(i))+M(i+1)*((T((i-1)*10+num)-Y
(i))^3)/(6*h(i))+(B(i)-((M(i)*(h(i)^2))/6))*(Y(i+1)-T((i-1)*10+num))/
h(i)+(B(i+1)-((M(i+1)*(h(i)^2))/6))*(T((i-1)*10+num)-Y(i))/h(i);
num=num+1;
%给出S的边界条件
S(1,111)=M(11)*((Y(12)-T(111))^3)/(6*h(11))+M(12)*((T(111)-Y(11))^3)/(6*h(11))+(B(11)-((M(11)*(h(11)^2))/6))*(Y(12)-T(111))/h(11)+(B(12)-((M(12)*(h(11)^2))/6))*(T(111)-Y(11))/h(11);
plot(T,S,'
m'
5.假定某天的气温变化记录如下表所示,试用最小二乘法找出这一天的气温变化的规律,试计算这一天的平均气温。
时刻
2
3
4
5
6
7
8
9
平均气温
15
14
16
18
20
23
25
28
13
17
19
21
22
24
31
32
41
29
27
考虑使用二次函数、三次函数、四次函数以及指数型的函数
计算相应的系数,估算误差,并作图比较各种函数之间的区别。
给定数据点{(xi,yi)}(i=1,2,3…,m)和一组函数gk(x),求出一系列数a1,a2,a3,…,an,使函数
满足
可以将上式,改写为如下的形式:
也就是要求出使上式成立的ak的值。
因此得到了使上式成立的必要条件:
最终可以得到以下的一个方程组:
解这个方程就可以得到最小二乘法的系数。
%根据采样数据组构建矩阵
A2=0:
24;
B=[15,14,14,14,14,15,16,18,20,20,23,25,28,31,32,41,29,27,25,24,22,20,18
17,16];
x=0:
0.1:
%构建最小二乘拟合法的系数矩阵A
A=zeros(25,6);
25;
6;
A(i,j)=A2(i)^(j-1);
x1=(A'
*A)\(A'
*B'
a1=zeros(1,6);
%初始化系数矩阵a1
%将系数矩阵转化为标准矩阵形式
a1(1,n)=x1(7-n,1);
plot(x,polyval(a1,x),'
b'
)holdon
最终的拟合结果如图(3)所示:
图(3)最小二乘法拟合结果
从图中可以得出结论,当天气温变化规律0-5时气温略有小波动,但基本保持不变,并于4时达到当天气温最低值;
5-15时气温不断升高,并于15时达到当天气温最高值;
15-24时气温不断降低。
利用不同的函数进行拟合,拟合的结果如下:
二次函数拟合、三次函数拟合和四次函数拟合:
x=[0:
1:
24];
y=[15,14,14,14,14,15,16,18,20,20,23,25,28,31,32,41,29,27,25,24,22,20,18,17,16];
a1=polyfit(x,y,2)%二次多项式拟合
a2=polyfit(x,y,3)%三次多项式拟合
a3=polyfit(x,y,4)%四次多项式拟合
b1=polyval(a1,x)
b2=polyval(a2,x)
b3=polyval(a3,x)
r1=sum((y-b1).^2)%二次多项式的误差平方和
r2=sum((y-b2).^2)%三次多项式的误差平方和
r3=sum((y-b3).^2)%四次多项式的误差平方和
plot(x,y,'
*'
)
plot(x,b1,'
)%二次多项式的图像
plot(x,b2,'
)%三次多项式的图像
plot(x,b3,'
)%四次多项式的图像
拟合图像如下图(4)所示:
图(4)多项式拟合的结果
r1=442.8523,r2=254.9607,r3=170.5871。
从误差上可以看出,多项式拟合的次数越高误差越小。
最终,拟合而成的多项式为:
二次多项式:
LinearmodelPoly2:
f(x)=p1*x^2+p2*x+p3
Coefficients(with95%confidencebounds):
p1=-0.1(-0.1401,-0.05989)
p2=2.775(1.779,3.772)
p3=7.815(2.651,12.98)
三次多项式:
LinearmodelPoly3:
f(x)=p1*x^3+p2*x^2+p3*x+p4
p1=-0.009389(-0.01435,-0.004426)
p2=0.238(0.05662,0.4194)
p3=-0.4038(-2.255,1.447)
p4=13.52(8.491,18.54)
四次多项式:
LinearmodelPoly4:
f(x)=p1*x^4+p2*x^3+p3*x^2+p4*x+p5
p1=0.001012(0.0003408,0.001683)
p2=-0.05796(-0.09044,-0.02548)
p3=0.9777(0.464,1.491)
p4=-4.168(-7.11,-1.226)
p5=17.2(12.32,22.08)
指数函数拟合具体算法:
x=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24];
opts=fitoptions('
Method'
Nonlinear'
Normalize'
On'
ftype=fittype('
a*exp(-b*(x-c).^2)'
options'
opts);
[fresult,gof]=fit(x'
y'
ftype)
plot(x,fresult(x),x,y,'
*'
拟合的结果如图(5)所示:
图(5)指数拟合的结果
指数拟合的参数为:
Generalmodel:
fresult(x)=a*exp(-b*(x-c).^2)
wherexisnormalizedbymean12andstd7.36
Coefficients(with95%confidencebounds):
a=28.98(26.34,31.62)
b=0.3431(0.2337,0.4526)
c=0.3008(0.1552,0.4463)
所以最终的拟合而成的指数函数为:
从两幅图的比较可以看出四次多项式的拟合效果比较理想。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 上机