1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
| function [A,C] = cqj_InverseMatrixOfTridiagonalMatrices(a,b,c)
C = []; flag = 0;
A = diag(a,-1)+diag(b)+diag(c,1); a = [0 a];
n = size(A,1); alpha(1) = b(1); if alpha(1) == 0 fprintf('算法2.1失败!正在尝试算法2.2...\n'); syms x alpha = x; flag = 1; end for i = 2:n t(i-1) = c(i-1)/alpha(i-1); alpha(i) = b(i)-a(i)*t(i-1); gamma(i) = a(i)/alpha(i-1); if alpha(i) == 0 fprintf('算法2.1失败!正在尝试算法2.2...\n'); syms x alpha = [alpha(1:i-1) x]; t = [t x]; gamma = [gamma x]; flag = 1; end end if alpha(n) == 0 fprintf('矩阵是奇异的!\n'); return; end
if flag P = 1; for i = 1:n P = P*alpha(i); end P = expand(P); P_poly = sym2poly(P); if ~P_poly(end) fprintf('矩阵是奇异的\n'); return; end end
C = zeros(n); if flag C = [C(1:end-1) x]; C = reshape(C,n,n); end C(n,n) = 1/alpha(n);
for i = n-1:-1:1 C(i,i) = 1/alpha(i)+t(i)*gamma(i+1)*C(i+1,i+1); end
for i = n:-1:2 for j = i-1:-1:1 C(i,j) = -gamma(j+1)*C(i,j+1); end end for i = n:-1:2 for j = i-1:-1:1 C(j,i) = -t(j)*C(j+1,i); end end
if flag C = expand(C); str = '@(x) ['; for i = 1:n for j = 1:n temp = factor(C(i,j)); str = [str ' ' char(prod(temp))]; end str = [str ';']; end str = [str ']']; Cfun = str2func(str); C = Cfun(0); return; end
|