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

sv12345

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

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

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

jimin

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为右端向量

sv12345

%   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？

jimin

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

sv12345

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

sv12345

我编了如下一个程序（主要是用追赶法求解），但就在程序算参数 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]？？？？？？

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

gghhjj

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. 
   
3. % ----------  example1  -----------
   
4. % Ax = b (x,b是列向量),当A是可逆矩阵时 x = A\b 产生该方程组的解
   
5. % 其算法基于LU分解相当于列主元Gauss消去法
   
6. A = [ 1  5 -9
   
7.       0  6  4
   
8.       1  1  1];
   
9. b = [-16 24 6]';
   
10. x1 = A\b
    
11. % [注] 该命令适用求解中小型稠密线性方程,而且性态是较好的(非病态),是最常用的命令
    
12. %      对于大型矩阵或病态的还有其它一些命令pcg,gmres,qmr等
    
13. 
    
14. % ----------  example2  -----------
    
15. % Ax = b ,当A是列满秩矩阵时 x = A\b 产生该方程组唯一的最小二乘解
    
16. % 其算法基于解法方程组 A'*A x = A'*b (见P125-126, 例11)
    
17. A = [2  1  1
    
18.      1 -2 -3
    
19.      3 -4  1
    
20.      1  3 -2];
    
21. b = [-4 5 3 -2]';
    
22. x2 = A\b
    
23. 
    
24. % -----------  example3  ----------
    
25. % 解矩阵方程AX = B ,当A可逆或列满秩
    
26. A = [ 1  3
    
27.       1  4];
    
28. B = [ 1  2
    
29.       3  4];
    
30. X = A\B
    
31. 
    
32. % exp3_2.m --- 矩阵分解命令的学习
    
33. 
    
34. % ---------- LU分解 ----------
    
35. % [简介] 设A是非奇异矩阵,则有如下列主元的LU分解
    
36. %             PA = LU
    
37. %        其中P是置换矩阵(又称排列矩阵),L是单位下三角矩阵且|l(ij)|<=1,U是上三角矩阵
    
38. %        矩阵A如果有了上述分解则求解线性方程组 Ax = b 就等价于分别求解两个三角方程组
    
39. %            Ly = Pb  和  Ux=y
    
40. %        而求解三角方程组是非常容易的.用这种方法求解方程组与列主元的Gauss消去法是等价的
    
41. %        但 LU 分解还有其它优点. 参见 P65 例9 和 P70 实验课题(一)
    
42. 
    
43. A = [10  -7  0
    
44.      -3   2  6
    
45.       5  -1  5];
    
46. [L U P] = lu(A)
    
47. 
    
48. % 如果再给b,则下面是解两个三角方程组的程序. 参见P49(3-15)式
    
49. b = [7 4 6]';
    
50. [n n]=size(A);
    
51. x = zeros(n,1);
    
52. 
    
53. % 前代求解:Ly = Pb(解用x储存)
    
54. b = P*b;
    
55. x(1) = b(1);
    
56. for k = 2:n
    
57.     x(k) = b(k)-L(k,1:k-1)*x(1:k-1);
    
58. end
    
59. 
    
60. % 回代求解:Ux = y (这里先y=x,解仍用x储存)
    
61. x(n) = x(n)/U(n,n);
    
62. for k = n-1:-1:1
    
63.     x(k) = ( x(k)-U(k,k+1:n)*x(k+1:n) ) / U(k,k);
    
64. end
    
65. x
    
66. 
    
67. % ---------- Cholesky分解 ----------
    
68. % [简介] 设A是(实对称)正定矩阵,则有如下唯一的分解
    
69. %             A = R'*R
    
70. %        其中 R 是上三角矩阵且主对角元全为正
    
71. 
    
72. % 例
    
73. A = [4    -1     1
    
74.     -1  4.25  2.75
    
75.      1  2.75   3.5];
    
76. R = chol(A)
    
77. 
    
78. % exp3_3.m  Gauss列主元消去法(示范程序)
    
79. function try_gauss
    
80. 
    
81. A = [ 1  5 -9
    
82.       0  6  4
    
83.       1  1  1];
    
84. 
    
85. b=[-3;10;3];
    
86. 
    
87. x=gauss(A,b)                         % 调用下面函数
    
88. 
    
89. 
    
90. function x = gauss(A,b)
    
91. % 输入: Ax = b 的系数矩阵 A (可逆)和右端项 b
    
92. % 输出:方程组的解 x
    
93. 
    
94. [n,n] = size(A);
    
95. x = zeros(n,1);
    
96. 
    
97. Aug = [A,b];                         % 增广矩阵
    
98. 
    
99. for k = 1:n-1
    
100.     [piv,r] = max(abs(Aug(k:n,k)));  % 找列主元所在子矩阵的行 r
     
101.     r = r + k - 1;                   % 列主元所在大矩阵的行
     
102.    
     
103.     if r>k
     
104.         Aug([k,r],:) = Aug([r,k],:); % 对增广矩阵实施行交换
     
105.     end
     
106.    
     
107.     if Aug(k,k)==0,
     
108.         error('这是奇异矩阵!')        % 程序遇到error会中断执行并显示其中的提示内容   
     
109.     end
     
110.    
     
111.     % 把增广矩阵消元成为上三角
     
112.     for p = k+1:n
     
113.         mult = Aug(p,k)/Aug(k,k);    % 消元乘子
     
114.         Aug(p,k:n+1) = Aug(p,k:n+1) - mult * Aug(k,k:n+1);
     
115.     end
     
116. end
     
117. 
     
118. % 解上三角方程组
     
119. A = Aug(:,1:n); b = Aug(:,n+1);
     
120. x(n) = b(n)/A(n,n);
     
121. for k = n-1:-1:1
     
122.     x(k) = ( b(k) - A(k,k+1:n)*x(k+1:n) ) / A(k,k);
     
123. end
     
124. 
     
125. % exp3_4.m --- 病态矩阵试验
     
126. 
     
127. % Hilbert矩阵是严重的病态的矩阵(见P56),求解病态方程组是相当困难
     
128. % 一般的方法都会失效,例如列主元Gauss消去法,要采取特殊的方法和技术
     
129. % 例如 (1) pcg (预处理共轭梯度法,只适用正定方程组),
     
130. %      (2) gmres (广义最残差法,适用一般的大型疏矩阵)
     
131. 
     
132. n = 15;
     
133. A = hilb(n);       % n 阶Hilbert矩阵
     
134. x = ones(n,1);     % 指定精确解全是1
     
135. b = A*x;
     
136. 
     
137. c = cond(A,inf);   % 求条件数,范数取最大值范数
     
138. 
     
139. x1 = A\b;          % 用Gauss列主元消去法求解或者调用你编的 x1 = gauss(A,b)
     
140. 
     
141. x2 = pcg(A,b);     % 用pcg法求解(A必须是正定矩阵)
     
142. 
     
143. x3=gmres(A,b);     % gmres法求解结果
     
144. 
     
145. % 显示结果
     
146. clc
     
147. fprintf('A 的条件数 =  %f\n',c)
     
148. fprintf('\nGauss法求解结果    pcg法求解结果       gmres法求解结果')
     
149. for i = 1:n
     
150.     fprintf('\n  %6.2f            %6.2f              %6.2f',x1(i),x2(i),x3(i))
     
151. end
     
152. 
     
153. % 注：范德蒙(Vandermonde)矩阵也是严重病态矩阵
     
154. %       A = vander(V),其中V是一向量
     
155. % 对于维数较大的V,看一下范德蒙矩阵的条件数
     
156. 
     
157. % exp3_5.m --- SOR迭代法收敛速度受松驰因子的影响试验
     
158. 
     
159. function try_sor_and_relaxation
     
160. % 参见 P64例8 和 P65表3-4
     
161. A = [ 2  -1   0   0
     
162.      -1   2  -1   0
     
163.       0  -1   2  -1
     
164.       0   0  -1   2];
     
165. b = [1 0 1 0]';
     
166. tol = 1e-6;
     
167. maxiter = 1000;
     
168. P = [1 0 1 0]';
     
169. 
     
170. clc
     
171. fprintf('** SOR迭代:松驰因子与迭代次数的关系 **\n'),
     
172. fprintf('\n    松驰因子         迭代次数')
     
173. fprintf('\n---------------------------------')
     
174. 
     
175. for W = 1:0.1:1.9
     
176.     [X,iternum] = sor(A,b,tol,maxiter,P,W);
     
177.     fprintf('\n      %3.2f            %4.0f',W,iternum)
     
178. end
     
179. % 通过以上结果,你认为最佳松驰因子大致为多少?
     
180. 
     
181. % ----------  说明  ----------
     
182. % 对于对称正定三对角矩阵(上面矩阵就是),其最佳松驰因子是可以求得的
     
183. %                  Wopt = 2 / ( 1+sqrt(1-rho^2) )
     
184. %                  其中rho是Jacobi迭代矩阵 Mj 的谱半径
     
185. % 下面求之,应该与我们实验的结果差不多吧,看看
     
186. D = diag(diag(A));
     
187. Mj = eye(4) - inv(D)*A;
     
188. rho = max(abs(eig(Mj)));
     
189. Wopt = 2 / ( 1+sqrt(1-rho^2) );
     
190. fprintf('\n---------------------------------')
     
191. fprintf('\n理论上最佳松驰因子是 Wopt = %3.2f',Wopt)
     
192. 
     
193. % ---------- SOR 迭代法 -----------
     
194. function [X,iternum] = sor(A,B,tol,maxiter,P,W)
     
195. % [X,iternum] = sor(A,B,tol,maxiter,P,W) SOR迭代法解线性方程组 AX=B
     
196. % 输入 ---- tol:     相对残向量的容差,当norm(B-AX)/norm(B) <= tol 终止迭代
     
197. %           maxiter: 最大迭代次数;
     
198. %           P:       初值
     
199. %           W:       松驰因子,必须 0 < W < 2(当 W = 1 时为G-S迭代法)
     
200. % 输出 ---- X:       解向量
     
201. %           iternum: 收敛迭代次数
     
202. 
     
203. N = length(B);  X = P;  tol = tol*norm(B);
     
204. for k = 1:maxiter
     
205.     for i=1:N
     
206.         if i == 1
     
207.             X(1) = (B(1)-A(1,2:N)*P(2:N))/A(1,1);
     
208.         elseif i == N
     
209.             X(N) = (B(N)-A(N,1:N-1)*X(1:N-1))/A(N,N);
     
210.         else
     
211.             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);
     
212.         end
     
213.         X(i) = (1-W)*P(i) + W*X(i);
     
214.     end
     
215.     if norm(B-A*X) <= tol, break,end;
     
216.     P=X;
     
217. end
     
218. iternum=k;
     
219. return

复制代码

本贴转自：星火学术论坛