Thomas algorithm tri diagonal matrix

조회 수: 21 (최근 30일)
Asha
Asha 2023년 12월 21일
편집: John D'Errico 2023년 12월 21일
In this code i want a tridiagonal matrix with different diagonal elements. but after some point i got stuck with the elments.i have attached the form of the matrix in image. please help me with this code.
h = 0.01;
x0 = 0;x1 = 1;
b = 1;
t0 = 0; t1 = 0.1;
k = 0.01;
x = x0: h: x1;
sizex = x;
t = t0: k: t1;
M = length(x);
N = length(t);
function YYp = Tridiag( p0,h,k,M)
alf=0.5;
B1=@(x,t)2;
B1x=@(x,t)0;
B2=@(x,t)4;
B2x=@(x,t)0;
R1=-3*B1x/h+6*B1/h^2;
R2=-3*B2x/h+6*B2/h^2;
S1=-12*B1/h^2;
S2=-12*B2/h^2;
T1=6*B1/h^2+3/h*B1x;
T2=6*B2/h^2+3/h*B2x;
A = zeros(M,M);
%Mvalue = M;
R(1)=0;
P=zeros(1,N);
Q=zeros(1,N-1);
R=zeros(1,N);
Y=zeros(1,N-1);
L1 = 0;
for i = 1:N
X(i) = (i-1) * h;
if X(i) <= alf
L1 = i;
end
end
for j=1:M
T(j)=(j-1)*k;
end
for j=2:M
for i = 2:L1-1
P(i, i) = S1;
end
for i = 2:L1-2
Q(i, i+1) = R1;
end
for i = 2:L1-2
R(i+1, i) = T1;
end
for i = L1+2:N
P(i, i) = S2;
end
for i = L1+2:N-1
Q(i, i+1) = R2;
end
for i = L1+2:N-1
R(i+1, i) = T2;
end
i=L1
P(i, i) = pi*S1;
Q(i, i+1) =6* R1;
R(i+1, i) = pi*T1;
i=L1+1
P(i, i) = pi*S2;
Q(i, i+1) =6* R2;
R(i+1, i) = pi*T2;
end
A = diag([4*R1 P(2:L1-1) R2-T2, 0]) + diag([0 R(2:L1-2)], 1) + diag([Q(2:L1-2) 0], -1);
ft = p0;
alpha0 = 0;
alpha1 = 0;
yp0 = alpha0;
yp1 = alpha1;
dm = h^2*ft;
dm(1) = dm(1) - yp0*(P(1));
dm(M) = dm(M) - yp1*(P(M));
Amatrix = A;
dm = dm';
sizeA = size(A);
sizedm = size(dm);
cc = A\dm;
ccm1 = 0 - 4*cc(1) - cc(2);
ccmp1 = 0 - 4*cc(M) - cc(M-1);
ccc = [ccm1 cc' ccmp1];
YYp = ccc(1:M) + 4*ccc(2:M+1) + ccc(3:M+2);
end

답변 (1개)

John D'Errico
John D'Errico 2023년 12월 21일
편집: John D'Errico 2023년 12월 21일
DON'T WRITE YOUR OWN CODE FOR A PROBLEM LIKE THIS.
Instead, just create the tridiagonal matrix, as a sparse matrix. Then use backslash to do the solve. All of that is literally trivial, requiring not much more than two calls, one to spdiags, and a second to backslash. Instead, you spent probably hours writing messy, heavily looped, unreadable, uncommented, and un-debuggable code that you got wrong. and now cannot fix.
Use MATLAB as it is designed to be used.
Create an nx3 array of the elements on the three diagonals. Call spdiags. Use backslash, or perhaps linsolve or decomposition as you wish.

카테고리

Help CenterFile Exchange에서 Linear Algebra에 대해 자세히 알아보기

제품


릴리스

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by