数值分析课程设计.docx
- 文档编号:15834457
- 上传时间:2023-07-08
- 格式:DOCX
- 页数:120
- 大小:471.98KB
数值分析课程设计.docx
《数值分析课程设计.docx》由会员分享,可在线阅读,更多相关《数值分析课程设计.docx(120页珍藏版)》请在冰点文库上搜索。
数值分析课程设计
课程设计
数值分析课程设计
201年月日
设计题目
数值分析课程设计
成绩
课
程
设
计
主
要
内
容
实验一.12,实验二.1234,实验三.123,实验四.1234,实验五.1234,实验六.24,实验七.567,实验八.126,共计25题。
题目通过matlib程序解答。
指
导
教
师
评
语
建议:
从学生的工作态度、工作量、设计(论文)的创造性、学术性、实用性及书面表达能力等方面给出评价。
签名:
200年月日
1.1水手、猴子和椰子问题:
五个水手带了一只猴子来到南太平洋的一个荒岛上,发现那里有一大堆椰子。
由于旅途的颠簸,大家都很疲惫,很快就入睡了。
第一个水手醒来后,把椰子平分成五堆,将多余的一只给了猴子,他私藏了一堆后便又去睡了。
第二、第三、第四、第五个水手也陆续起来,和第一个水手一样,把椰子分成五堆,恰多一只猴子,私藏一堆,再去入睡,天亮以后,大家把余下的椰子重新等分成五堆,每人分一堆,正好余一只再给猴子,试问原先共有几只椰子?
试分析椰子数目的变化规律,利用逆向递推的方法求解这一问题。
解:
算法分析:
根据题意可以采用递归算法进行计算。
设设椰子起初的数目为x0,第一至第五次猴子在夜里藏椰子后,椰子的数目分别为x0,x1,x2,x3,x4,再设最后每个人分得x个椰子,设计以下算法:
n=input('n=');
forx=1:
n
t=5*x+1;
fork=1:
5
t=5*t/4+1;
end
ift==fix(t)
break
end
end
disp([x,t])
结果:
输入n=10000,得到没人分1024个椰子,椰子总数为15621。
改变x的范围,得到204731246;307146871;.......
等解。
1.2当
时,选择稳定的算法计算积分
.
解:
题目分析:
观察积分中的函数,发现下列特性:
由此得到递归式:
并且由上式可知求
时,
的误差的影响被缩小了。
当n=100时
的近似值为0。
则得到matli程序:
fprintf('稳定算法:
\n')
y0=0;
n=100;
plot(n,y0,'r*');
holdon
fprintf('y[100]=%10.6f',y0);
while
(1)
y1=1/10*[(1-exp(-n))/n-y0];
fprintf('y[%10.0f]=%10.6f',n-1,y1);plot(n-1,y1,'r*')
if(n<=1)break;end
y0=y1;n=n-1;
ifmod(n,3)==0,fprintf('\n'),end,end
得到图形:
由此看出这是稳定的算法。
2.1小行星轨道问题:
一天文学家要确定一颗小行星绕太阳运行的轨道,他在轨道平面内建立以太阳为原点的直角坐标系,在五个不同的对小行星作了五次观察,测得轨道上五个点的坐标数据(单位:
万公里)如下表所示:
P1
P2
P3
P4
P5
X坐标
53605
58460
62859
66662
68894
Y坐标
6026
11179
16954
23492
68894
由开普勒第一定律知,小行星轨道为一椭圆,椭圆的一般方程可表示为:
现需要建立椭圆的方程以供研究。
(1)分别将五个点的数据代入椭圆一般方程中,写出五个待定系数满足的等式,整理后写出线性方程组
以及方程组的系数矩阵和右端项b;
(2)用MARLAB求低阶方程的指令A\b求出待定系数
;
(3)分别用直接法、Jacobi迭代法、Gauss-Seidel迭代法求出待定系数
.
解:
(1)分别将五个点的数据代入椭圆一般方程中,写出五个待定系数满足的等式,整理后写出线性方程组
以及方程组的系数矩阵和右端项b;
线性方程组为:
=
先算出系数矩阵。
使用MATLIB计算。
X0=[5360558460628596666268894];
Y0=[606211179169542349268894];
A=zeros(5);
fori=1:
5
A(i,1)=X0(i)^2;A(i,2)=2*X0(i)*Y0(i);A(i,3)=Y0(i)^2;A(i,4)=2*X0(i);A(i,5)=2*Y0(i);
end;
formatlongg;A
得到结果
A=
28734960256499070203674784410721012124
341757160013070486801249700411169202235
395125388121314229722874381112571833908
4443822244313204740855187406413332446984
474638323694927664724746383236137788137788
b=[-1-1-1-1-1]
(2)用MARLAB求低阶方程的指令A\b求出待定系数;
B=[-1-1-1-1-1]';formatlongg;x=A\B
x=
-8.06820280371841e-011
-7.63620099622306e-011
-3.0801152978055e-010
-8.89025159419867e-006
2.02829368401655e-005
(3)分别用直接法、Jacobi迭代法、Gauss-Seidel迭代法求出待定系数.
直接法可使用用Gauss列主元消元法解线性方程组Ax=b
%用Gauss列主元消元法解线性方程组Ax=b
%RA,RB分别表示系数矩阵A和增广矩阵B的秩
%N表示向量b的维数,X是解向量
function[RA,RB,N,X’]=Gauss(A,b)
B=[Ab];
N=length(b);
RA=rank(A);
RB=rank(B);
Diff=RB-RA;
ifDiff>0
disp('因为RA~=RB,所以此方程组无解.')
return
end
ifRA==RB
ifRA==N
disp('因为RA=RB=N,所以此方程组有唯一解.')
X=zeros(N,1);
C=zeros(1,N+1);
fori=1:
N-1
[Yk]=max(abs(B(i:
N,i)));
C=B(i,:
);
B(i,:
)=B(k+i-1,:
);
B(k+i-1,:
)=C;
forj=i+1:
N
m=B(j,i)/B(i,i);
B(j,i:
N+1)=B(j,i:
N+1)-m*B(i,i:
N+1);
end
end
b=B(1:
N,N+1);
A=B(1:
N,1:
N);
X(N)=b(N)/A(N,N);
fori=N-1:
-1:
1
X(i)=(b(i)-sum(A(i,i+1:
N)*X(i+1:
N)))/A(i,i);
end
else
disp('因为RA=RB end end X=X’; End 因为RA=RB=N,所以此方程组有唯一解. RA= 5 RB= 5 N= 5 X= 1.0e-004* 0.000006437215780 -0.000001957441606 0.000002799074569 -0.254776815358136 -0.001104956851916 Jacobi迭代法解线性方程组Ax=b %用Jacobi迭代法解线性方程组Ax=b %Norm: 范数的名称,Norm=1,2,inf; %error: 近似解x的误差; %Max: 迭代的最大次数; functionX=Jacobi(A,b,X0,Norm,Error,Max) [NN]=size(A); X=zeros(N,1); fori=1: Max forj=1: N X(j)=(b(j)-A(j,[1: j-1,j+1: N])*X0([1: j-1,j+1: N]))/A(j,j); end X, errX=norm(X-X0,Norm); X0=X; iferrX X1=A\b; disp('迭代次数i,精确解X1和近似解X分别是: ') formatlong i,X1,X, return end end iferrX>=Error disp('请注意: Jacobi迭代次数已经超过最大迭代次数Max.') end end 在MATLAB命令窗口输入: A=[28734960256499070203674784410721012124;341757160013070486801249700411169202235; 3951253881213142297228743811612571833908;4443822244313204740855187406413332446984; 474638323694927664724746383236137788137788]; b=[-1-1-1-1-1]'; X0=[00000]'; Jacobi(A,b,X0,inf,0.001,100); 结果 X= 1.0e-005* -0.000034800813758 -0.000076508244513 -0.000347900972187 -0.750052503675257 -0.725752605451854 迭代次数i,精确解X1和近似解X分别是: i= 1 X1= 1.0e-004* 0.000006437215780 -0.000001957441606 0.000002799074569 -0.254776815358137 -0.001104956851916 X= 1.0e-005* -0.000034800813758 -0.000076508244513 -0.000347900972187 -0.750052503675257 -0.725752605451854 用Gauss-Seidel迭代法解线性方程组Ax=b %用Gauss-Seidel迭代法解线性方程组Ax=b %Norm: 范数的名称,Norm=1,2,inf; %Error是近似解X的误差;Max是迭代的最大次数 functionX=G_S(A,b,X0,Norm,Error,Max) [NN]=size(A); X=zeros(N,1); fori=1: Max forj=1: N X(j)=0; ifj>1 X(j)=X(j)+A(j,1: j-1)*X(1: j-1); end ifj X(j)=X(j)+A(j,j+1: N)*X0(j+1: N); end X(j)=(b(j)-X(j))/A(j,j); end X, errX=norm(X-X0,Norm); X0=X; iferrX X1=A\b; disp('迭代次数i,精确解X1和近似解X分别是: ') formatlong i,X1,X, return end end iferrX>=Error disp('请注意: Gauss-Seidel迭代次数已经超过最大迭代次数Max.') end end X= 1.0e-004* -0.000003480081376 0.000001448627970 0.000002306743862 -0.002590235506932 -0.129368843047805 迭代次数i,精确解X1和近似解X分别是: i= 1 X1= 1.0e-004* 0.000006437215780 -0.000001957441606 0.000002799074569 -0.254776815358137 -0.001104956851916 X= 1.0e-004* -0.000003480081376 0.000001448627970 0.000002306743862 -0.002590235506932 -0.129368843047805 2.2 (1)用Gauss列主元消去法、Gauss按比例列主元消去法、Cholesky分解求解下列线性方程组,并彼此互相验证。 (2)判断用Jacobi迭代法、Gauss-Seidel迭代法、SOR法(分别取 )解下列线性方程组的收敛性.若收敛,再用Jacobi迭代法、Gauss-Seidel迭代法、SOR法(分别取 )分别解线性方程组,并比较各种方法的收敛速度. (3)用Cholesky分解求解下列线性方程组 (方程组的精确解是 ) 解: (1)Gauss列主元消去法已在上题中给出。 1.用按比例主元消元法解线性方程组Ax=b %用按比例主元消元法解线性方程组Ax=b %RA,RB分别表示系数矩阵A和增广矩阵B的秩 %N表示向量b的维数,X是解向量 function[RA,RB,N,X]=Ratio(A,b) B=[Ab]; N=length(b); RA=rank(A); RB=rank(B); RDiff=RB-RA; ifRDiff>0 disp('因为RA~=RB,所以此方程组无解.') return end ifRA==RB ifRA==N disp('因为RA=RB=n,所以此方程组有唯一解.') X=zeros(N,1); C=zeros(1,N+1); T=0; D=B'; s=max(abs(D(1: N,: ))); S=s'; fori=1: N-1 [Yk]=max((abs(B(i: N,i)))./S(i: N)); K=min(k); C=B(i,: ); B(i,: )=B(K+i-1,: ); B(K+i-1,: )=C; T=S(i); S(i)=S(K+i-1); S(K+i-1)=T; forj=i+1: N m=B(j,i)/B(i,i); B(j,i: N+1)=B(k,i: N+1)-m*B(i,i: N+1); end end b=B(1: N,N+1); A=B(1: N,1: N); X(N)=b(N)/A(N,N); fori=N-1: -1: 1 X(i)=(b(i)-sum(A(i,i+1: N)*X(i+1: N)))/A(i,i); end else disp('请注意: 因为RA=RB end X=X’ End 2.用Cholesky分解(平方根法)解线性方程组 %用Cholesky分解解线性方程组 %A是方程组的系数矩阵,b是方程组的右边向量 functionCholesky(A,b) [NN]=size(A); RA=rank(A); ifRA~=N disp('因为A的n阶行列式D等于零,所以A不能进行Cholesky分解.A的秩RA如下: ') RA, return end ifA~=A' disp('因为A不是对称矩阵,所以A不能进行Cholesky分解.') return end ifRA==N fori=1: N h(i)=det(A(1: i,1: i)); end D1=h(1: N); fori=1: N ifh(1,i)<=0 disp('因为A的i阶主子式不大于零,所以A不能进行Cholesky分解.A的各阶顺序主子式值Dl依次如下: ') D1, return end end ifh(1,i)>0 disp('因为A是对称正定矩阵,所以A能进行Cholesky分解.A的下三角矩阵G和方程组的解X如下: ') forj=1: N U(1,j)=A(1,j); end fori=2: N L(1,1)=1; L(i,i)=1; L(i,1)=A(i,1)/U(1,1); end fork=2: N fori=2: N forj=2: N ifi>j L(1,1)=1; L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1); L(i,k)=(A(i,k)-L(i,1: k-1)*U(1: k-1,k))/U(k,k); else U(k,j)=A(k,j)-L(k,1: k-1)*U(1: k-1,j); end end end end D1=eye(N); fori=1: N D1(i,i)=sqrt(U(i,i)); end G=L*D1; Y (1)=b (1)/G(1,1); fori=2: N S1=0; forj=1: i-1 S1=S1+G(i,j)*Y(j); end Y(i)=(b(i)-S1)/G(i,i); end GT=G'; X(N)=Y(N)/GT(N,N); fori=N-1: -1: 1 S2=0; forj=1: N-i S2=S2+GT(i,i+j)*X(i+j); end X(i)=(Y(i)-S2)/GT(i,i); end G,X, end end end A=[1-121;-130-3;209-6;1-3-619];b=[1357]'; [RA,RB,n,x]=liezy(A,b); [RA,RB,n,x]=bilizy(A,b); cholesky(A,b) 结果: 列主元 因为RA=RB=n,所以此方程组有唯一解. RA= 4 RB= 4 n= 4 x= -8.000000000000005 0.333333333333332 3.666666666666668 2.000000000000000 比例主元 因为RA=RB=n,所以此方程组有唯一解. RA= 4 RB= 4 n= 4 x= -8.000000000000000 0.333333333333333 3.666666666666667 2.000000000000000 cholesky分解 S1= 0 S1= 0 S1= 0 x= 003.6666666666666732.000000000000003 x= 00.3333333333333293.6666666666666732.000000000000003 x= -8.0000000000000200.3333333333333293.6666666666666732.000000000000003 -8.0000000000000200.3333333333333293.6666666666666732.000000000000003 对比可知,三种方法可相互验证。 (2) 先用谱半径判别Jacobi迭代法的收敛性 %用谱半径测试Jacobi迭代法的收敛性 %A是方程组的系数矩阵 functionJacobiTest(A) [mn]=size(A); ifm~=n disp('系数矩阵必须为方阵') return end D=diag(diag(A)); I=eye(n,n); B=I-inv(D)*A; E=eig(B); SRH=norm(E,inf); ifSRH>=1 disp('因为谱半径不小于1,所以Jacobi迭代序列发散,谱半径SRH和迭代矩阵的所有特征值如下: ') SRH,Eig=E’, else disp('因为谱半径小于1,所以Jacobi迭代序列收敛,谱半径SRH和迭代矩阵的所有特征值如下: ') SRH,Eig=E’, end end 用谱半径测试Gauss-Seidel迭代法的收敛性 %A是方程组的系数矩阵 functionG_STest(A) [mn]=size(A); ifm~=n disp('系数矩阵必须为方阵') return end D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); B=inv(D-L)*U; E=eig(B); SRH=norm(E,inf) ifSRH>=1 disp('因为谱半径不小于1,所以Jacobi迭代序列发散,谱半径SRH和迭代矩阵的所有特征值如下: ') SRH,Eig=E’, else disp('因为谱半径小于1,所以Jacobi迭代序列收敛,谱半径SRH和迭代矩阵的所有特征值如下: ') SRH,Eig=E’, end end 用谱半径测试SOR方法的收敛性 %A是方程组的系数矩阵,w松弛因子 functionSORTest(A,w) [mn]=size(A); ifm~=n d
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 课程设计