1、数值分析课程设计课 程 设 计 数值分析课程设计 201 年 月 日设计题目 数值分析课程设计成绩课程设计主要内容实验一.1 2 ，实验二.1 2 3 4，实验三.1 2 3，实验四.1 2 3 4，实验五.1 2 3 4，实验六.2 4，实验七.5 6 7，实验八.1 2 6，共计25题。题目通过matlib程序解答。指导教师评语建议：从学生的工作态度、工作量、设计（论文）的创造性、学术性、实用性及书面表达能力等方面给出评价。签名： 200 年 月 日1.1 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。
2、第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。解：算法分析：根据题意可以采用递归算法进行计算。设设椰子起初的数目为x0,第一至第五次猴子在夜里藏椰子后，椰子的数目分别为x0，x1，x2，x3，x4，再设最后每个人分得x个椰子，设计以下算法：n=input(n=);for x=1:n t=5*x+1;
3、 for k=1:5 t=5*t/4+1;endif t=fix(t) break endend disp(x,t) 结果：输入n=10000，得到没人分1024个椰子，椰子总数为15621。改变x的范围，得到 2047 31246；3071 46871；.等解。1.2 当时，选择稳定的算法计算积分.解：题目分析：观察积分中的函数，发现下列特性： 由此得到递归式： 并且由上式可知求时，的误差的影响被缩小了。当n=100时的近似值为0。则得到matli程序：fprintf(稳定算法:n)y0=0;n=100;plot(n,y0,r*);hold onfprintf(y100=%10.6f,y0)
4、;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 0 disp (因为RA=RB，所以此方程组无解.) returnendif RA = RB if RA = N disp (因为RA=RB=N，所以此方程组有唯一解.) X = zeros(N, 1); C = zeros(1, N+1); for i = 1 : N-1 Y k = max( abs(B(i: N, i) ); C=B(i,:); B(i, :) = B(k+i-1, :); B(k+i-1, :) =
5、 C; for j = 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); for i = 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=RBN，所以此方程组有无穷多解.) endendX=X;End因为RA=RB=N，所以此方程组有唯一解.RA = 5R
6、B = 5N = 5X = 1.0e-004 * 0.000006437215780 -0.000001957441606 0.000002799074569 -0.254776815358136 -0.001104956851916Jacobi 迭代法解线性方程组Ax=b%用Jacobi 迭代法解线性方程组Ax = b %Norm：范数的名称，Norm=1,2,inf;%error：近似解x的误差；%Max：迭代的最大次数;function X = Jacobi(A, b, X0, Norm, Error, Max)N N = size ( A );X = zeros ( N, 1 );fo
7、r i = 1 : Max for j = 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; if errX = Error disp ( 请注意：Jacobi迭代次数已经超过最大迭代次数Max. )endend在MATLAB命令窗口输入： A=2873496025 649907020 36747844 107210 12124;3417571600 1307048680 124970041 116920 2235;39
8、51253881 2131422972 287438116 125718 33908;4443822244 3132047408 551874064 133324 46984;4746383236 9492766472 4746383236 137788 137788;b=-1 -1 -1 -1 -1;X0=0 0 0 0 0;Jacobi(A, b, X0, inf, 0.001, 100);结果 X = 1.0e-005 * -0.000034800813758 -0.000076508244513 -0.000347900972187 -0.750052503675257 -0.7257
9、52605451854迭代次数i, 精确解 X1 和近似解 X 分别是:i = 1X1 = 1.0e-004 * 0.000006437215780 -0.000001957441606 0.000002799074569 -0.254776815358137 -0.001104956851916X = 1.0e-005 * -0.000034800813758 -0.000076508244513 -0.000347900972187 -0.750052503675257 -0.725752605451854用Gauss-Seidel迭代法解线性方程组Ax=b %用Gauss-Seidel迭
10、代法解线性方程组Ax=b%Norm:范数的名称，Norm=1,2,inf;%Error是近似解X的误差；Max是迭代的最大次数function X = G_S(A, b, X0, Norm, Error, Max)N N = size ( A );X = zeros ( N, 1 );for i = 1 : Max for j = 1 : N X(j) = 0; if j 1 X(j) = X(j) + A(j, 1:j-1) * X(1:j-1); end if j N X(j) = X(j) + A(j, j+1:N) * X0(j+1:N); end X(j) = ( b(j) - X(
11、j) ) / A(j, j); end X, errX = norm ( X-X0, Norm ); X0 = X; if errX = Error disp ( 请注意：Gauss-Seidel迭代次数已经超过最大迭代次数Max. )endendX = 1.0e-004 * -0.000003480081376 0.000001448627970 0.000002306743862 -0.002590235506932 -0.129368843047805迭代次数i, 精确解 X1 和近似解 X 分别是:i = 1X1 = 1.0e-004 * 0.000006437215780 -0.00
12、0001957441606 0.000002799074569 -0.254776815358137 -0.001104956851916X = 1.0e-004 * -0.000003480081376 0.000001448627970 0.000002306743862 -0.002590235506932 -0.1293688430478052.2 (1) 用Gauss列主元消去法、Gauss按比例列主元消去法、Cholesky分解求解下列线性方程组，并彼此互相验证。(2) 判断用Jacobi 迭代法、Gauss-Seidel迭代法、SOR法(分别取)解下列线性方程组的收敛性. 若收敛
13、，再用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 = A b;N = length(b);RA = rank(A);RB = rank(B);RDiff = RB
14、- RA;if RDiff 0 disp(因为RA=RB，所以此方程组无解.) returnendif RA = RB if RA = 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; for i = 1 : N-1 Y k = 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=
15、S(i); S(i) = S(K+i-1); S(K+i-1) = T; for j = 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); for i = 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=RBn，所以此方程组有无穷多解.)endX=XEnd2.用Choles
16、ky分解(平方根法)解线性方程组%用Cholesky分解解线性方程组% A是方程组的系数矩阵，b是方程组的右边向量function Cholesky(A, b)N N = size(A);RA = rank(A); if RA = N disp(因为A的n阶行列式D等于零，所以A不能进行Cholesky分解. A的秩RA如下:) RA, returnendif A = A disp(因为A不是对称矩阵，所以A不能进行Cholesky分解.) returnendif RA = N for i= 1 : N h(i) = det( A(1:i, 1:i) ); end D1=h(1:N); for
17、 i=1:N if h(1,i) 0 disp(因为A是对称正定矩阵，所以A能进行Cholesky分解. A的下三角矩阵G和方程组的解X如下：) for j = 1 : N U(1, j) = A(1, j); end for i = 2 : N L(1,1) = 1; L(i,i) = 1; L(i, 1) = A(i,1) / U(1,1); end for k = 2 : N for i = 2 : N for j = 2 : N if i 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) =
18、 ( 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); for i = 1 : N D1(i, i) = sqrt(U(i, i); end G=L*D1; Y(1) = b(1) / G(1,1); for i = 2 : N S1 = 0; for j = 1 : i-1 S1 = S1 + G(i,j) * Y(j); end Y(i) = (b(i)-S1) / G(i, i); end GT=G; X
19、(N)=Y(N) / GT(N, N); for i = N-1 : -1 : 1 S2 = 0; for j = 1 : N-i S2 = S2 + GT(i,i+j) * X(i+j); end X(i) = (Y(i)-S2) / GT(i,i); end G, X, endendendA=1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19;b=1 3 5 7; RA,RB,n,x=liezy(A,b)；RA,RB,n,x=bilizy(A,b)；cholesky(A,b) 结果：列主元因为RA=RB=n,所以此方程组有唯一解.RA = 4RB = 4n = 4
20、x = -8.000000000000005 0.333333333333332 3.666666666666668 2.000000000000000比例主元因为RA=RB=n,所以此方程组有唯一解.RA = 4RB = 4n = 4x = -8.000000000000000 0.333333333333333 3.666666666666667 2.000000000000000cholesky分解S1 = 0S1 = 0S1 = 0x = 0 0 3.666666666666673 2.000000000000003x = 0 0.333333333333329 3.666666666
21、666673 2.000000000000003x = -8.000000000000020 0.333333333333329 3.666666666666673 2.000000000000003 -8.000000000000020 0.333333333333329 3.666666666666673 2.000000000000003 对比可知，三种方法可相互验证。（2）先用谱半径判别Jacobi 迭代法的收敛性%用谱半径测试Jacobi迭代法的收敛性% A是方程组的系数矩阵function JacobiTest(A)m n = size(A);if m = n disp ( 系数矩
22、阵必须为方阵 ) returnendD=diag(diag(A);I=eye(n,n);B=I-inv(D)*A;E=eig(B);SRH=norm(E,inf);if SRH = 1 disp ( 因为谱半径不小于1，所以Jacobi迭代序列发散, 谱半径SRH和迭代矩阵的所有特征值如下： ) SRH, Eig=E,else disp ( 因为谱半径小于1，所以Jacobi迭代序列收敛,谱半径SRH和迭代矩阵的所有特征值如下： ) SRH, Eig=E,endend用谱半径测试Gauss-Seidel迭代法的收敛性% A是方程组的系数矩阵function G_STest(A)m n = si
23、ze(A);if m = n disp ( 系数矩阵必须为方阵 ) returnendD=diag(diag(A);U=-triu(A,1);L=-tril(A,-1);B=inv(D-L)*U;E=eig(B);SRH = norm(E, inf )if SRH = 1 disp ( 因为谱半径不小于1，所以Jacobi迭代序列发散,谱半径SRH和迭代矩阵的所有特征值如下： ) SRH, Eig=E,else disp ( 因为谱半径小于1，所以Jacobi迭代序列收敛,谱半径SRH和迭代矩阵的所有特征值如下： ) SRH, Eig=E,endend用谱半径测试SOR方法的收敛性% A是方程组的系数矩阵，w松弛因子function SORTest(A, w)m n = size(A);if m = n d
