论坛 › Matlab讨论区 › 求matlab或vb关于追赶法的源程序

sv12345 发表于 2006-10-14 15:26

求matlab或vb关于追赶法的源程序

小弟最近在作论文时遇到了“追赶法求解三对角矩阵”的问题，在此小弟向各位大哥大姐求追赶法的matlab或vb源程序，哪位大侠快救救俺啊！小弟要是这块解决不了，下面就没有办法进行了啊 ！！！！！！！！！！！

===========eight============
请不要使用“跪求”等字眼
==========================

[ 本帖最后由 eight 于 2007-2-6 23:09 编辑 ]

jimin 发表于 2006-10-14 15:43

function x = tridisolve(a,b,c,d)
%   TRIDISOLVESolve tridiagonal system of equations.
%   x = TRIDISOLVE(a,b,c,d) solves the system of linear equations
%   b(1)*x(1) + c(1)*x(2) = d(1),
%   a(j-1)*x(j-1) + b(j)*x(j) + c(j)*x(j+1) = d(j), j = 2:n-1,
%   a(n-1)*x(n-1) + b(n)*x(n) = d(n).
%
%   The algorithm does not use pivoting, so the results might
%   be inaccurate if abs(b) is much smaller than abs(a)+abs(c).
%   More robust, but slower, alternatives with pivoting are:
%   x = T\d where T = diag(a,-1) + diag(b,0) + diag(c,1)
%   x = S\d where S = spdiags([ b ],[-1 0 1],n,n)

x = d;
n = length(x);

for j = 1:n-1
   mu = a(j)/b(j);
   b(j+1) = b(j+1) - mu*c(j);
   x(j+1) = x(j+1) - mu*x(j);
end

x(n) = x(n)/b(n);
for j = n-1:-1:1
   x(j) = (x(j)-c(j)*x(j+1))/b(j);
end
说明：a为主对角线下面的次对角线，b为主对角线，c为主对角线上面的次对角线，d为右端向量

sv12345 发表于 2006-10-14 19:42

请教高人

%   More robust, but slower, alternatives with pivoting are:
%   x = T\d where T = diag(a,-1) + diag(b,0) + diag(c,1)
%   x = S\d where S = spdiags([ b ],[-1 0 1],n,n)
这三句是不是表示在“abs(b) is much smaller than abs(a)+abs(c)”这种情况下，程序要想继续运行，则需要用x = T\d或 x = S\d 来代替x = d？

jimin 发表于 2006-10-14 20:06

实际问题提出的三对角方程组，系数矩阵往往严格对角占优，所以不选主元，就可保证顺利，稳定运行
至于那些不是对角占优的abs(b) is much smaller than abs(a)+abs(c)”，此方法求出的误差就会比较大，此方法可能会失效
T = diag(a,-1) + diag(b,0) + diag(c,1)这个是转换成正常的格式
x = T\d 直接解估计是不行的，用高斯消元法或其他方法解吧
你先弄清楚你要解的是不是严格对角占优吧
个人意见，仅供参考

sv12345 发表于 2006-10-14 21:09

致谢

jimin(敏敏) 真是位人心的大侠啊！小弟对你的佩服之心犹如滔滔江水连绵不绝啊：）呵呵！谢过了啊

sv12345 发表于 2006-11-4 17:56

求助：各位大侠帮我来调一下程序

我编了如下一个程序（主要是用追赶法求解），但就在程序算参数 b 的时候，每次 b(1) 的值正确，而 b(2) 以后的值都不对（我想实现 b(i)=2*(Y(i)+L(i+1)*I(i)*Y(i+1)/(L(i)*I(i+1)))，为啥还会出错？郁闷），求大侠快来救命！！！小弟在线等!!!!!!!!!

输入参数:
gaoding(2,147000,1000,10,0.244,0.244,0.184)
输入各跨的长度：
10
10
10
输入个跨惯性矩：
5e-5
5e-5
5e-5

运行程序：
%————————定义输入参数————————
function gaoding(n,PB,omega,bita,D0,DS,DC)
%————————定义已知条件————————
for j=1:n+1
    L(j)=input('请输入各跨的长度\n');
end
for j=1:n+1
    I(j)=input('请输入各跨的惯性矩\n');
end
for j=1:n
    pi=3.14159265;
    P(1)=PB;
    w=omega;
    bt=bita*pi/180;
    E=2e11;
    q=w*sin(bt);
    P(j+1)=P(j)-0.5*w*L(j)*cos(bt)-0.5*w*L(j+1)*cos(bt);
end
%———————求u(j),X(j),Y(j),Z(j)———————
for j=1:n+1
    u(j)=L(j)*sqrt(P(j)/(E*I(j)))/2;
    X(j)=3*(tan(u(j))-u(j))/u(j)^3;
    Y(j)=1.5*(0.5/u(j)-1/tan(2*u(j)))/u(j);
    Z(j)=3*(1/(sin(2*u(j)))-0.5/u(j))/u(j);
end
%————在求d时要用到个e(0),故在此转化一下————
e(1:n+1)=0.5*(D0-DS);
f(1:n+1)=;
%————追赶法求解弯矩主程序————
for j=1:n-1
    a(j)=Z(j+1);
    c(j)=L(j+1)*I(j)*Z(j+1)/(L(j)*I(j+1));
end
a=;
c=;
for j=1:n

    b(j)=2*(Y(j)+L(j+1)*I(j)*Y(j+1)/(L(j)*I(j+1)));             ％这就是b的一般表达式
    d(j)=-q*L(j)^2*X(j)/4-q*L(j+1)^2*L(j+1)*I(j)*X(j+1)/(4*L(j)*I(j+1))+6*E*I(j)*(e(j)-f(j))/L(j)^2-6*E*I(j)*(e(j+1)-e(j))/(L(j)*L(j+1));
end
b=;
d=';
m = d;
k= length(m);

for j = 1:k-1
   mu = a(j)/b(j);
   b(j+1) = b(j+1) - mu*c(j);
   m(j+1) = m(j+1) - mu*m(j);
end
m(k) = m(k)/b(k);
for j = k-1:-1:1
   m(j) = (m(j)-c(j)*m(j+1))/b(j);
end
％————输出参数————
Pa=P(1)*e(1)/L(1)-m(1)/L(1)-q*L(1)/2;
P,u,X,Y,Z,a,b,c,d,m,Pa

得到结果：
P =

1.0e+005 *

    1.4700    1.3701    1.2703


u =

    0.6062    0.5853    0.5635


X =

    1.1727    1.1591    1.1458


Y =

    1.1141    1.1052    1.0964


Z =

    1.2028    1.1867    1.1710


a =

    1.1867


b =

    4.4385    4.0859


c =

    1.1867


d =

1.0e+003 *

   -3.0509
   -3.0156


m =

-543.4104
-538.4268


Pa =

-207.3387

按我编程的原意，算出来 b=，但为何域运行出来是 b=？？？？？？

有哪位好心的大侠原意帮我看看，调试一下？小弟不胜感激啊！！！

gghhjj 发表于 2007-4-19 07:26

    function y=chase(a,b,c,d)   %追赶法
    =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 = [ 15 -9
      064
      111];
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 = [211
   1 -2 -3
   3 -41
   13 -2];
b = [-4 5 3 -2]';
x2 = A\b

% -----------example3----------
% 解矩阵方程AX = B ,当A可逆或列满秩
A = [ 13
      14];
B = [ 12
      34];
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-70
   -3   26
      5-15];
= lu(A)

% 如果再给b,则下面是解两个三角方程组的程序. 参见P49(3-15)式
b = ';
=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
    -14.252.75
   12.75   3.5];
R = chol(A)

% exp3_3.mGauss列主元消去法(示范程序)
function try_gauss

A = [ 15 -9
      064
      111];

b=[-3;10;3];

x=gauss(A,b)                         % 调用下面函数


function x = gauss(A,b)
% 输入: Ax = b 的系数矩阵 A (可逆)和右端项 b
% 输出:方程组的解 x

= size(A);
x = zeros(n,1);

Aug = ;                         % 增广矩阵

for k = 1:n-1
    = max(abs(Aug(k:n,k)));% 找列主元所在子矩阵的行 r
    r = r + k - 1;                   % 列主元所在大矩阵的行
   
    if r>k
      Aug(,:) = Aug(,:); % 对增广矩阵实施行交换
    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 = ';
tol = 1e-6;
maxiter = 1000;
P = ';

clc
fprintf('** SOR迭代:松驰因子与迭代次数的关系 **\n'),
fprintf('\n    松驰因子         迭代次数')
fprintf('\n---------------------------------')

for W = 1:0.1:1.9
    = 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 = sor(A,B,tol,maxiter,P,W)
% = 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

本贴转自：星火学术论坛

查看完整版本: 求matlab或vb关于追赶法的源程序