计算方法上机题.docx
- 文档编号:29915355
- 上传时间:2023-08-03
- 格式:DOCX
- 页数:17
- 大小:161.69KB
计算方法上机题.docx
《计算方法上机题.docx》由会员分享,可在线阅读,更多相关《计算方法上机题.docx(17页珍藏版)》请在冰豆网上搜索。
计算方法上机题
1.计算
要求误差小于
给出实现算法。
算法组织:
利用绝对值的递增相加可以减小误差
具体算法:
%%公式求和%%
S=0;
k=100000;
%%利用for循环对公式进行求和
forn=1:
100000
S=S+1/(k*k);
k=k-1;
end
%%精确到小数点后6位,以保证误差范围
vpa(S,7)
运行程序后,得到的结果是1.644924,精确到了小数点后6位。
2.编写实现对N阶非奇矩阵A进行LU分解的程序。
算法组织:
(1)根据定理可知若N阶矩阵的各阶顺序主子式行列式不为零,则存在唯一的单位下三角矩阵L和上三角矩阵U,满足A=LU,因此首先需要判断A中的元素akk(k=0,1,2…,n)不为零。
(2)根据如下的式子计算L和U中的元素
当i=1时,得到u1j=a1j(j=1,2,3,…,n),而当j=1时,得li1=ai1/u11,由此可以计算出L的第一列元素和U的第一行元素。
然后可以依据下式计算L和U中的其他元素。
具体算法:
[n,n]=size(A);%检测A的阶数
L=eye(n);%初始化单位下三角阵L
U=zeros(n,n);%初始化上三角阵U
fori=1:
n
U(1,i)=A(1,i);%给上三角阵U的第一行赋值
L(i,1)=A(i,1)/U(1,1);%给单位下三角阵L的第一列赋值
end
fori=2:
n
forj=i:
n
fork=1:
i-1
M(k)=L(i,k)*U(k,j);
end
U(i,j)=A(i,j)-sum(M);%给上三角阵U的第i行赋值
end
forj=i+1:
n
fork=1:
i-1
M(k)=L(j,k)*U(k,i);
end
L(j,i)=(A(j,i)-sum(M))/U(i,i);%给单位下三角阵L的第i列赋值
end
end
%分别输出L,U的矩阵
L
U
3.编写程序实现大规模方程组的列主元高斯消去法程序,并对所附的方程组进行求解。
具体的要求参见所附的相关文档。
针对本专业中所碰到的实际问题,提炼一个使用方程组进行求解的例子,并对求解过程进行分析、求解。
具体算法:
fid=fopen('fun001.dat','r');%获取数据文件文件头
headfile=fread(fid,3,'uint32');%确认数据文件标识号是否正确
id=headfile
(1);
ifid~=hex2dec('F1E1D1A0')
fprintf('不是正确的数据文件.\n');
return;
end
%获取数据文件版本号
ver=headfile
(2);
ifver==hex2dec('101')
%获取方程阶数
n=headfile(3)
%获取稀疏矩阵
A=fread(fid,[n,n],'float32');
%获取右端系数
B=fread(fid,[n,1],'float32');
%对系数矩阵A和右端向量的方程组使用高斯方法求解
x=GAUSSPP(A,B,n);
end
ifver==hex2dec('102')
%系数矩阵为条状对角阵
%获取方程阶数
n=headfile(3)
%获取稀疏矩阵
A=fread(fid,[n,n],'float32');
%获取右端系数
B=fread(fid,[n,1],'float32');
%对系数矩阵A和右端向量的方程组使用高斯方法求解
x=GAUSSPP(A,B,n);
end
ifver==hex2dec('201')
%系数矩阵为压缩格式下的条状对角阵
%获取方程阶数
n=headfile(3)
%获取条状矩阵的上下带宽
headfile2=fread(fid,2,'uint32');
p=headfile2
(1);
q=headfile2
(2);
%压缩系数矩阵
A=fread(fid,[p+q+1,n],'float32');
%获取右端系数
B=fread(fid,[n,1],'float32');
%对系数矩阵A和右端向量的方程组使用高斯方法求解
x=GAUSSPP(A,B,n,p,q);
end
fclose(fid);
GAUSSPP(A,B,n)的程序如下:
functionx=GAUSSPP(A,B,n)
fork=1:
n-1
fori=k+1:
n
l(i,k)=A(i,k)/A(k,k);
forj=k:
n
A(i,j)=A(i,j)-l(i,k)*A(k,j);
end
B(i)=B(i)-l(i,k)*B(k);
end
end
x(n)=B(n)/A(n,n);
fork=n-1:
-1:
1
s=0;
forj=k+1:
n
s=s+A(k,j)*x(j);
end
x(k)=(B(k)-s)/A(k,k);
end
GAUSSPP(A,B,n,p,q)的程序如下:
functionx=GAUSSPP(A,B,n,p,q)
fork=1:
n-1
fori=k+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);
end
B(i)=B(i)-l(i,p-i+k+1)*B(k);
end
end
x(n)=B(n)/A(n,p+1);
fork=n-1:
-1:
1
s=0;
forj=1:
min(q,n-k)
s=s+A(k,1+p+j)*x(j+k);
end
x(k)=(B(k)-s)/A(k,p+1);
end
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:
12;
A(i,j)=YB^(j-1);
end
YB=YB+10;
end
%根据数据构建矩阵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);
forn=1:
12
a(1,n)=x(13-n,1);
end
e1=polyval(a,[2005,2015])
T=1900:
2010
plot(Y,B,'r')%红色表示原数据图像
holdon
plot(T,polyval(a,T),'g')%绿色表示多项式插值近似图像
holdon
最终得到的结果是
e1=
342456
即2005年的产量为342,2015年的产量为456。
多项式插值的结果如图
(1)所示:
图
(1)多项式插值的图形
图中红色曲线表示原始数据,绿色曲线表示多项式插值的数据。
2.三次样条函数插值
算法组织:
由于分段多项式插值不够光滑,因此需要较高次的分段多项式,常用三项多项式,即在每一个子区间上市一个三次多项式,这样的分段三次多项式就是三次样条函数。
具体算法:
n=12;
%构造数据矩阵Y
Y=1900:
10:
2010;
%构造数据矩阵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];
%构造步长矩阵h
h=zeros(1,n-1);
fork=1:
11
h(k)=Y(k+1)-Y(k);
end
%构造M的系数矩阵A
A=zeros(n);
%由第三类边界条件可得
A(1,2)=-2;
%由第三类边界条件可得
A(n,n-1)=-2;
fori=1:
n
A(i,i)=2;
end
fori=2:
n-1
A(i,i+1)=h(i)/(h(i-1)+h(i));
A(i,i-1)=h(i-1)/(h(i-1)+h(i));
end
%构造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:
n-1
d(j)=6*ChaShang_3(Y,B,j-1,j,j+1);
end
M=Gauss_pivot(A,d);
%构造在不同区间内三次样条拟合结果的矩阵S
T=1900:
2010;
S=zeros(1,(n-1)*10+1);
fori=1:
n-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;
end
end
%给出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');
holdon
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
41
29
27
25
24
22
20
18
17
16
考虑使用二次函数、三次函数、四次函数以及指数型的函数
计算相应的系数,估算误差,并作图比较各种函数之间的区别。
算法组织:
给定数据点{(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:
24;
%构建最小二乘拟合法的系数矩阵A
A=zeros(25,6);
fori=1:
25;
forj=1:
6;
A(i,j)=A2(i)^(j-1);
end
end
x1=(A'*A)\(A'*B');
a1=zeros(1,6);%初始化系数矩阵a1
%将系数矩阵转化为标准矩阵形式
forn=1:
6
a1(1,n)=x1(7-n,1);
end
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,'*')
holdon
plot(x,b1,'r')%二次多项式的图像
holdon
plot(x,b2,'g')%三次多项式的图像
holdon
plot(x,b3,'b')%四次多项式的图像
拟合图像如下图(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
Coefficients(with95%confidencebounds):
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
Coefficients(with95%confidencebounds):
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];
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];
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文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 上机
![提示](https://static.bdocx.com/images/bang_tan.gif)