声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4059|回复: 6

[编程技巧] 求matlab或vb关于追赶法的源程序

[复制链接]
发表于 2006-10-14 15:26 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

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

[ 本帖最后由 eight 于 2007-2-6 23:09 编辑 ]
回复
分享到:

使用道具 举报

发表于 2006-10-14 15:43 | 显示全部楼层
function x = tridisolve(a,b,c,d)
%   TRIDISOLVE  Solve 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([[a; 0] b [0; c]],[-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为右端向量

评分

1

查看全部评分

 楼主| 发表于 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([[a; 0] b [0; c]],[-1 0 1],n,n)
这三句是不是表示在“abs(b) is much smaller than abs(a)+abs(c)”这种情况下,程序要想继续运行,则需要用x = T\d或 x = S\d 来代替x = d?
发表于 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 直接解估计是不行的,用高斯消元法或其他方法解吧
你先弄清楚你要解的是不是严格对角占优吧
个人意见,仅供参考
 楼主| 发表于 2006-10-14 21:09 | 显示全部楼层

致谢

jimin  (敏敏) 真是位人心的大侠啊!小弟对你的佩服之心犹如滔滔江水连绵不绝啊:)呵呵!谢过了啊
 楼主| 发表于 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)=[0,e(1:n)];
%————追赶法求解弯矩主程序————
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=[a(1:n-1)];
c=[c(1:n-1)];
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=[b(1:n)];
d=[d(1:n)]';
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=[4.4385    4.4032],但为何域运行出来是 b=[4.4385    4.0859]??????

有哪位好心的大侠原意帮我看看,调试一下?小弟不胜感激啊!!!
发表于 2007-4-19 07:26 | 显示全部楼层
  1.     function y=chase(a,b,c,d)   %追赶法
  2.     [n,n]=size(d);
  3.     d(1)=d(1)/b(1);
  4.     c(1)=c(1)/b(1);
  5.     for k=2:n
  6.         b(k)=b(k)-a(k)*c(k-1);
  7.         c(k)=c(k)/b(k);
  8.         d(k)=(d(k)-a(k)*d(k-1))/b(k);
  9.     end
  10.     for k=n-1:-1:1
  11.         d(k)=d(k)-c(k)*d(k+1);
  12.     end
  13.     for k=1:n
  14.         disp(d(k));
  15.     end
复制代码


再加几个其他的数值分析方面的程序

  1. % exp3_1.m --- 解线性方程组左除命令‘\’的学习

  2. % ----------  example1  -----------
  3. % Ax = b (x,b是列向量),当A是可逆矩阵时 x = A\b 产生该方程组的解
  4. % 其算法基于LU分解相当于列主元Gauss消去法
  5. A = [ 1  5 -9
  6.       0  6  4
  7.       1  1  1];
  8. b = [-16 24 6]';
  9. x1 = A\b
  10. % [注] 该命令适用求解中小型稠密线性方程,而且性态是较好的(非病态),是最常用的命令
  11. %      对于大型矩阵或病态的还有其它一些命令pcg,gmres,qmr等

  12. % ----------  example2  -----------
  13. % Ax = b ,当A是列满秩矩阵时 x = A\b 产生该方程组唯一的最小二乘解
  14. % 其算法基于解法方程组 A'*A x = A'*b (见P125-126, 例11)
  15. A = [2  1  1
  16.      1 -2 -3
  17.      3 -4  1
  18.      1  3 -2];
  19. b = [-4 5 3 -2]';
  20. x2 = A\b

  21. % -----------  example3  ----------
  22. % 解矩阵方程AX = B ,当A可逆或列满秩
  23. A = [ 1  3
  24.       1  4];
  25. B = [ 1  2
  26.       3  4];
  27. X = A\B

  28. % exp3_2.m --- 矩阵分解命令的学习

  29. % ---------- LU分解 ----------
  30. % [简介] 设A是非奇异矩阵,则有如下列主元的LU分解
  31. %             PA = LU
  32. %        其中P是置换矩阵(又称排列矩阵),L是单位下三角矩阵且|l(ij)|<=1,U是上三角矩阵
  33. %        矩阵A如果有了上述分解则求解线性方程组 Ax = b 就等价于分别求解两个三角方程组
  34. %            Ly = Pb  和  Ux=y
  35. %        而求解三角方程组是非常容易的.用这种方法求解方程组与列主元的Gauss消去法是等价的
  36. %        但 LU 分解还有其它优点. 参见 P65 例9 和 P70 实验课题(一)

  37. A = [10  -7  0
  38.      -3   2  6
  39.       5  -1  5];
  40. [L U P] = lu(A)

  41. % 如果再给b,则下面是解两个三角方程组的程序. 参见P49(3-15)式
  42. b = [7 4 6]';
  43. [n n]=size(A);
  44. x = zeros(n,1);

  45. % 前代求解:Ly = Pb(解用x储存)
  46. b = P*b;
  47. x(1) = b(1);
  48. for k = 2:n
  49.     x(k) = b(k)-L(k,1:k-1)*x(1:k-1);
  50. end

  51. % 回代求解:Ux = y (这里先y=x,解仍用x储存)
  52. x(n) = x(n)/U(n,n);
  53. for k = n-1:-1:1
  54.     x(k) = ( x(k)-U(k,k+1:n)*x(k+1:n) ) / U(k,k);
  55. end
  56. x

  57. % ---------- Cholesky分解 ----------
  58. % [简介] 设A是(实对称)正定矩阵,则有如下唯一的分解
  59. %             A = R'*R
  60. %        其中 R 是上三角矩阵且主对角元全为正

  61. % 例
  62. A = [4    -1     1
  63.     -1  4.25  2.75
  64.      1  2.75   3.5];
  65. R = chol(A)

  66. % exp3_3.m  Gauss列主元消去法(示范程序)
  67. function try_gauss

  68. A = [ 1  5 -9
  69.       0  6  4
  70.       1  1  1];

  71. b=[-3;10;3];

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


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

  76. [n,n] = size(A);
  77. x = zeros(n,1);

  78. Aug = [A,b];                         % 增广矩阵

  79. for k = 1:n-1
  80.     [piv,r] = max(abs(Aug(k:n,k)));  % 找列主元所在子矩阵的行 r
  81.     r = r + k - 1;                   % 列主元所在大矩阵的行
  82.    
  83.     if r>k
  84.         Aug([k,r],:) = Aug([r,k],:); % 对增广矩阵实施行交换
  85.     end
  86.    
  87.     if Aug(k,k)==0,
  88.         error('这是奇异矩阵!')        % 程序遇到error会中断执行并显示其中的提示内容   
  89.     end
  90.    
  91.     % 把增广矩阵消元成为上三角
  92.     for p = k+1:n
  93.         mult = Aug(p,k)/Aug(k,k);    % 消元乘子
  94.         Aug(p,k:n+1) = Aug(p,k:n+1) - mult * Aug(k,k:n+1);
  95.     end
  96. end

  97. % 解上三角方程组
  98. A = Aug(:,1:n); b = Aug(:,n+1);
  99. x(n) = b(n)/A(n,n);
  100. for k = n-1:-1:1
  101.     x(k) = ( b(k) - A(k,k+1:n)*x(k+1:n) ) / A(k,k);
  102. end

  103. % exp3_4.m --- 病态矩阵试验

  104. % Hilbert矩阵是严重的病态的矩阵(见P56),求解病态方程组是相当困难
  105. % 一般的方法都会失效,例如列主元Gauss消去法,要采取特殊的方法和技术
  106. % 例如 (1) pcg (预处理共轭梯度法,只适用正定方程组),
  107. %      (2) gmres (广义最残差法,适用一般的大型疏矩阵)

  108. n = 15;
  109. A = hilb(n);       % n 阶Hilbert矩阵
  110. x = ones(n,1);     % 指定精确解全是1
  111. b = A*x;

  112. c = cond(A,inf);   % 求条件数,范数取最大值范数

  113. x1 = A\b;          % 用Gauss列主元消去法求解或者调用你编的 x1 = gauss(A,b)

  114. x2 = pcg(A,b);     % 用pcg法求解(A必须是正定矩阵)

  115. x3=gmres(A,b);     % gmres法求解结果

  116. % 显示结果
  117. clc
  118. fprintf('A 的条件数 =  %f\n',c)
  119. fprintf('\nGauss法求解结果    pcg法求解结果       gmres法求解结果')
  120. for i = 1:n
  121.     fprintf('\n  %6.2f            %6.2f              %6.2f',x1(i),x2(i),x3(i))
  122. end

  123. % 注:范德蒙(Vandermonde)矩阵也是严重病态矩阵
  124. %       A = vander(V),其中V是一向量
  125. % 对于维数较大的V,看一下范德蒙矩阵的条件数

  126. % exp3_5.m --- SOR迭代法收敛速度受松驰因子的影响试验

  127. function try_sor_and_relaxation
  128. % 参见 P64例8 和 P65表3-4
  129. A = [ 2  -1   0   0
  130.      -1   2  -1   0
  131.       0  -1   2  -1
  132.       0   0  -1   2];
  133. b = [1 0 1 0]';
  134. tol = 1e-6;
  135. maxiter = 1000;
  136. P = [1 0 1 0]';

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

  141. for W = 1:0.1:1.9
  142.     [X,iternum] = sor(A,b,tol,maxiter,P,W);
  143.     fprintf('\n      %3.2f            %4.0f',W,iternum)
  144. end
  145. % 通过以上结果,你认为最佳松驰因子大致为多少?

  146. % ----------  说明  ----------
  147. % 对于对称正定三对角矩阵(上面矩阵就是),其最佳松驰因子是可以求得的
  148. %                  Wopt = 2 / ( 1+sqrt(1-rho^2) )
  149. %                  其中rho是Jacobi迭代矩阵 Mj 的谱半径
  150. % 下面求之,应该与我们实验的结果差不多吧,看看
  151. D = diag(diag(A));
  152. Mj = eye(4) - inv(D)*A;
  153. rho = max(abs(eig(Mj)));
  154. Wopt = 2 / ( 1+sqrt(1-rho^2) );
  155. fprintf('\n---------------------------------')
  156. fprintf('\n理论上最佳松驰因子是 Wopt = %3.2f',Wopt)

  157. % ---------- SOR 迭代法 -----------
  158. function [X,iternum] = sor(A,B,tol,maxiter,P,W)
  159. % [X,iternum] = sor(A,B,tol,maxiter,P,W) SOR迭代法解线性方程组 AX=B
  160. % 输入 ---- tol:     相对残向量的容差,当norm(B-AX)/norm(B) <= tol 终止迭代
  161. %           maxiter: 最大迭代次数;
  162. %           P:       初值
  163. %           W:       松驰因子,必须 0 < W < 2(当 W = 1 时为G-S迭代法)
  164. % 输出 ---- X:       解向量
  165. %           iternum: 收敛迭代次数

  166. N = length(B);  X = P;  tol = tol*norm(B);
  167. for k = 1:maxiter
  168.     for i=1:N
  169.         if i == 1
  170.             X(1) = (B(1)-A(1,2:N)*P(2:N))/A(1,1);
  171.         elseif i == N
  172.             X(N) = (B(N)-A(N,1:N-1)*X(1:N-1))/A(N,N);
  173.         else
  174.             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);
  175.         end
  176.         X(i) = (1-W)*P(i) + W*X(i);
  177.     end
  178.     if norm(B-A*X) <= tol, break,end;
  179.     P=X;
  180. end
  181. iternum=k;
  182. return
复制代码


本贴转自:星火学术论坛

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 01:20 , Processed in 0.068941 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表