- function y=chase(a,b,c,d) %追赶法
- [n,n]=size(d);
- d(1)=d(1)/b(1);
- c(1)=c(1)/b(1);
- for k=2:n
- b(k)=b(k)-a(k)*c(k-1);
- c(k)=c(k)/b(k);
- d(k)=(d(k)-a(k)*d(k-1))/b(k);
- end
- for k=n-1:-1:1
- d(k)=d(k)-c(k)*d(k+1);
- end
- for k=1:n
- disp(d(k));
- end
复制代码
再加几个其他的数值分析方面的程序
- % exp3_1.m --- 解线性方程组左除命令‘\’的学习
- % ---------- example1 -----------
- % Ax = b (x,b是列向量),当A是可逆矩阵时 x = A\b 产生该方程组的解
- % 其算法基于LU分解相当于列主元Gauss消去法
- A = [ 1 5 -9
- 0 6 4
- 1 1 1];
- b = [-16 24 6]';
- x1 = A\b
- % [注] 该命令适用求解中小型稠密线性方程,而且性态是较好的(非病态),是最常用的命令
- % 对于大型矩阵或病态的还有其它一些命令pcg,gmres,qmr等
- % ---------- example2 -----------
- % Ax = b ,当A是列满秩矩阵时 x = A\b 产生该方程组唯一的最小二乘解
- % 其算法基于解法方程组 A'*A x = A'*b (见P125-126, 例11)
- A = [2 1 1
- 1 -2 -3
- 3 -4 1
- 1 3 -2];
- b = [-4 5 3 -2]';
- x2 = A\b
- % ----------- example3 ----------
- % 解矩阵方程AX = B ,当A可逆或列满秩
- A = [ 1 3
- 1 4];
- B = [ 1 2
- 3 4];
- X = A\B
- % exp3_2.m --- 矩阵分解命令的学习
- % ---------- LU分解 ----------
- % [简介] 设A是非奇异矩阵,则有如下列主元的LU分解
- % PA = LU
- % 其中P是置换矩阵(又称排列矩阵),L是单位下三角矩阵且|l(ij)|<=1,U是上三角矩阵
- % 矩阵A如果有了上述分解则求解线性方程组 Ax = b 就等价于分别求解两个三角方程组
- % Ly = Pb 和 Ux=y
- % 而求解三角方程组是非常容易的.用这种方法求解方程组与列主元的Gauss消去法是等价的
- % 但 LU 分解还有其它优点. 参见 P65 例9 和 P70 实验课题(一)
- A = [10 -7 0
- -3 2 6
- 5 -1 5];
- [L U P] = lu(A)
- % 如果再给b,则下面是解两个三角方程组的程序. 参见P49(3-15)式
- b = [7 4 6]';
- [n n]=size(A);
- x = zeros(n,1);
- % 前代求解:Ly = Pb(解用x储存)
- b = P*b;
- x(1) = b(1);
- for k = 2:n
- x(k) = b(k)-L(k,1:k-1)*x(1:k-1);
- end
- % 回代求解:Ux = y (这里先y=x,解仍用x储存)
- x(n) = x(n)/U(n,n);
- for k = n-1:-1:1
- x(k) = ( x(k)-U(k,k+1:n)*x(k+1:n) ) / U(k,k);
- end
- x
- % ---------- Cholesky分解 ----------
- % [简介] 设A是(实对称)正定矩阵,则有如下唯一的分解
- % A = R'*R
- % 其中 R 是上三角矩阵且主对角元全为正
- % 例
- A = [4 -1 1
- -1 4.25 2.75
- 1 2.75 3.5];
- R = chol(A)
- % exp3_3.m Gauss列主元消去法(示范程序)
- function try_gauss
- A = [ 1 5 -9
- 0 6 4
- 1 1 1];
- b=[-3;10;3];
- x=gauss(A,b) % 调用下面函数
- function x = gauss(A,b)
- % 输入: Ax = b 的系数矩阵 A (可逆)和右端项 b
- % 输出:方程组的解 x
- [n,n] = size(A);
- x = zeros(n,1);
- Aug = [A,b]; % 增广矩阵
- for k = 1:n-1
- [piv,r] = max(abs(Aug(k:n,k))); % 找列主元所在子矩阵的行 r
- r = r + k - 1; % 列主元所在大矩阵的行
-
- if r>k
- Aug([k,r],:) = Aug([r,k],:); % 对增广矩阵实施行交换
- end
-
- if Aug(k,k)==0,
- error('这是奇异矩阵!') % 程序遇到error会中断执行并显示其中的提示内容
- end
-
- % 把增广矩阵消元成为上三角
- for p = k+1:n
- mult = Aug(p,k)/Aug(k,k); % 消元乘子
- Aug(p,k:n+1) = Aug(p,k:n+1) - mult * Aug(k,k:n+1);
- end
- end
- % 解上三角方程组
- A = Aug(:,1:n); b = Aug(:,n+1);
- x(n) = b(n)/A(n,n);
- for k = n-1:-1:1
- x(k) = ( b(k) - A(k,k+1:n)*x(k+1:n) ) / A(k,k);
- end
- % exp3_4.m --- 病态矩阵试验
- % Hilbert矩阵是严重的病态的矩阵(见P56),求解病态方程组是相当困难
- % 一般的方法都会失效,例如列主元Gauss消去法,要采取特殊的方法和技术
- % 例如 (1) pcg (预处理共轭梯度法,只适用正定方程组),
- % (2) gmres (广义最残差法,适用一般的大型疏矩阵)
- n = 15;
- A = hilb(n); % n 阶Hilbert矩阵
- x = ones(n,1); % 指定精确解全是1
- b = A*x;
- c = cond(A,inf); % 求条件数,范数取最大值范数
- x1 = A\b; % 用Gauss列主元消去法求解或者调用你编的 x1 = gauss(A,b)
- x2 = pcg(A,b); % 用pcg法求解(A必须是正定矩阵)
- x3=gmres(A,b); % gmres法求解结果
- % 显示结果
- clc
- fprintf('A 的条件数 = %f\n',c)
- fprintf('\nGauss法求解结果 pcg法求解结果 gmres法求解结果')
- for i = 1:n
- fprintf('\n %6.2f %6.2f %6.2f',x1(i),x2(i),x3(i))
- end
- % 注:范德蒙(Vandermonde)矩阵也是严重病态矩阵
- % A = vander(V),其中V是一向量
- % 对于维数较大的V,看一下范德蒙矩阵的条件数
- % exp3_5.m --- SOR迭代法收敛速度受松驰因子的影响试验
- function try_sor_and_relaxation
- % 参见 P64例8 和 P65表3-4
- A = [ 2 -1 0 0
- -1 2 -1 0
- 0 -1 2 -1
- 0 0 -1 2];
- b = [1 0 1 0]';
- tol = 1e-6;
- maxiter = 1000;
- P = [1 0 1 0]';
- clc
- fprintf('** SOR迭代:松驰因子与迭代次数的关系 **\n'),
- fprintf('\n 松驰因子 迭代次数')
- fprintf('\n---------------------------------')
- for W = 1:0.1:1.9
- [X,iternum] = sor(A,b,tol,maxiter,P,W);
- fprintf('\n %3.2f %4.0f',W,iternum)
- end
- % 通过以上结果,你认为最佳松驰因子大致为多少?
- % ---------- 说明 ----------
- % 对于对称正定三对角矩阵(上面矩阵就是),其最佳松驰因子是可以求得的
- % Wopt = 2 / ( 1+sqrt(1-rho^2) )
- % 其中rho是Jacobi迭代矩阵 Mj 的谱半径
- % 下面求之,应该与我们实验的结果差不多吧,看看
- D = diag(diag(A));
- Mj = eye(4) - inv(D)*A;
- rho = max(abs(eig(Mj)));
- Wopt = 2 / ( 1+sqrt(1-rho^2) );
- fprintf('\n---------------------------------')
- fprintf('\n理论上最佳松驰因子是 Wopt = %3.2f',Wopt)
- % ---------- SOR 迭代法 -----------
- function [X,iternum] = sor(A,B,tol,maxiter,P,W)
- % [X,iternum] = sor(A,B,tol,maxiter,P,W) SOR迭代法解线性方程组 AX=B
- % 输入 ---- tol: 相对残向量的容差,当norm(B-AX)/norm(B) <= tol 终止迭代
- % maxiter: 最大迭代次数;
- % P: 初值
- % W: 松驰因子,必须 0 < W < 2(当 W = 1 时为G-S迭代法)
- % 输出 ---- X: 解向量
- % iternum: 收敛迭代次数
- N = length(B); X = P; tol = tol*norm(B);
- for k = 1:maxiter
- for i=1:N
- if i == 1
- X(1) = (B(1)-A(1,2:N)*P(2:N))/A(1,1);
- elseif i == N
- X(N) = (B(N)-A(N,1:N-1)*X(1:N-1))/A(N,N);
- else
- X(i) = (B(i)-A(i,1:i-1)*X(1:i-1)-A(i,i+1:N)*P(i+1:N))/A(i,i);
- end
- X(i) = (1-W)*P(i) + W*X(i);
- end
- if norm(B-A*X) <= tol, break,end;
- P=X;
- end
- iternum=k;
- return
复制代码
本贴转自:星火学术论坛 |