数值分析课程设计.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=RBend
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;
iferrXX1=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
ifjX(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;
iferrXX1=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=RBend
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