西安交通大学计算方法上机实习.docx
- 文档编号:24115418
- 上传时间:2023-05-24
- 格式:DOCX
- 页数:38
- 大小:259.03KB
西安交通大学计算方法上机实习.docx
《西安交通大学计算方法上机实习.docx》由会员分享,可在线阅读,更多相关《西安交通大学计算方法上机实习.docx(38页珍藏版)》请在冰豆网上搜索。
西安交通大学计算方法上机实习
计算方法上机报告
学院:
机械工程学院
班级:
姓名:
学号:
目录
上机实习题目
(一)1
第1题1
第2题1
第3题3
第4题7
第5题13
上机实习题目
(二)22
第6题22
第7题24
第8题26
小结29
上机实习题目
(一)
1.计算
要求误差小于
给出实现算法。
算法思想:
采用循环语句累加求和,程序中k从100000循环到1,因为在浮点数系中,浮点数运算存在舍入误差,大数和小数相加时会出现“大数吃小数”的现象,而此采用从100000开始求和直至1结束。
源程序:
>>formatlong
s=0;%初始化结果值
fork=100000:
-1:
1%设置循环结构,为避免大数吃小数,采用倒着运算
s=s+1/(k^2);%累加表达式
end
s%先输出一次精度较高的值
s=sprintf('%.6f',s)%在小数点后某七位四舍五入,误差小于10^-6
输出结果:
2.编写实现对N阶非奇矩阵A进行LU分解的程序。
算法思想:
在线性代数中,LU分解是矩阵分解的一种,可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积。
Gauss消去法的消去过程将矩阵逐步变换成一个新的矩阵。
从矩阵运算的观点来看,消去的每一步就等价于用一个初等三角阵去左乘方程的两端。
对于非奇矩阵A可以表示成一个单位下三角矩阵与一个上三角矩阵的乘积,即A=LU。
参照《计算方法教程》P38页,LU的计算公式:
源程序:
%第2题对N阶非奇矩阵A进行LU分解
function[L,U]=dierti_2_LU(A)
nj=size(A);%求矩阵A的行数和列数
ifnj
(1)~=nj
(2)%检查矩阵A是否为n阶方阵
disp('不是方阵,错误');
elseifnj
(1)~=rank(A)%检查矩阵A是否为n阶非奇异方阵
disp('A为奇异阵,错误');
else
n=nj
(1);
U=zeros(n);%初始化U矩阵,建立与A矩阵相同阶数的零矩阵
L=eye(n);%初始化L矩阵,建立与A矩阵相同阶数的单位矩阵
U(1,:
)=A(1,:
);%求U矩阵的第一行
L(2:
n,1)=A(2:
n,1)/U(1,1);%从第二行开始,求L矩阵的第一列
fork=2:
n%设置循环结构
U(k,k:
n)=A(k,k:
n)-L(k,1:
k-1)*U(1:
k-1,k:
n);
L(k+1:
n,k)=(A(k+1:
n,k)-L(k+1:
n,1:
k-1)*U(1:
k-1,k))/U(k,k);
%%以上两个式子,是利用公式和来求L和%U矩阵的剩余项
end
end
输出结果验证
>>A=[114;2-15;212];
[L,U]=dierti_2_LU(A)
输出为:
3.编写程序实现大规模方程组的列主元高斯消去法程序,并对所附的方程组进行求解。
具体的要求参见所附的相关文。
针对本专业中所碰到的实际问题,提炼一个使用方程组进行求解的例子,并对求解过程进行分析、求解。
算法思想:
本算法主要依据为列主元高斯消去法。
消去法的中心就是“降维”,即:
将求解n元方程组的问题转化为先解n-1元方程组,一旦这个n-1元方程组的解取得,则剩余的一个未知量自然可以求得。
这样逐步减少未知量个数的方法,便是求解多元方程组的一个重要思想。
一般地,消去法包含两个主要步骤:
消去和回带。
另外,稀疏矩阵具有大量零元的矩阵,同样可用列主元高斯消去法。
源程序:
(1)列主元Gauss消去法的源代码
functionx=disanti_3_gauss(a,b)
%本函数用disanti_3_gauss高斯列主元消去法求解线性方程组
n=size(a,2);%求矩阵a的列数
x=linspace(0,0,n)';%列向量x全部赋值为0
fork=1:
1:
n-1
uk=k;
fori=k:
1:
n
ifabs(a(uk,k)) 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 西安交通大学 计算方法 上机 实习