TCC Connection design matrix

I have a problem with forming a matrix in MATLAB according to the following formulas:

Could someone please write a complete MATLAB code that generates this matrix for an arbitrary number of elements

[Torsten]

What are C(0,1), a(0), C(n,n+1) and a(n+1) if V is assumed to be of size nxn ?

[dpb]

The three generating equations don't seem to be self-consistent without an additional conditon of for which i,j each is to be used...otherwise, why couldn't the middle one of V(i,j) be used for all i, j?

[Star Strider]

Something is clearly wrong --

C = sym('C',[4 4], 'real')

C = 

D = sym('D',[1 4], 'real')

D = 

a = sym('a',[1 4], 'real')

a = 

for i = 1:4
    if i > 1
        v(i, i-1) = -C(i-1,i) * a(i-1);
        V1 = v(i, i-1)
    end
    if i > 1 & i < 4
        v(i, i) = (C(i-1, i) + C(i, i+1) + D(i)) * a(i);
        V_diag = v(i, i)
    end
    if i < 4
        v(i, i+1) = -C(i, i+1) * a(i+1);
        V2 = v(i, i+1)
    end
end

V2 = 

V1 = 

V_diag = 

V2 = 

V1 = 

V_diag = 

V2 = 

V1 = 

v

v = 

So 's' is only defined for i=[2 3] as well, for the same reasons.

None of the 'i' subsceipts in the image look like they should be 'j' so there is only one actual index counter.

.

[Torsten]

@Star Strider

Usually these tridiagonal matrices result from the discretization of 1d- partial differential equations.

The values for C(0,1), a(0), C(n,n+1) and a(n+1) (which would make "s" computable also for i = 1 and i = 4) come from boundary conditions.

[Star Strider]

Noted.

In this instance, I would expect those values to be supplied as well.

[dpb]

"...there is only one actual index counter. "

My old eyes got fooled into reading the ",i" as a "j"

[Filip]

For boundary, i.e., outer elements of the matrix, the terms that would require the index i-1 do not exist and are therefore not included in the calculation of that element.

[Star Strider]

It would seem that they have to be elements of the matrix eventually. They should be supplied, or somehow noted.

[dpb]

By the defining equations, they have to have some definition in order to build V(1,1), for example. They may all be identically zero in value, but they must exist/be defined.

[Filip]

These elements must have a defined value, but elements labeled, for example, C01 or C45 should be ignored, because they definitely do not exist and do not participate in the calculation.”

The essence is that the matrix must be complete with its valid elements, while nonexistent or invalid elements should be formally skipped.

[dpb]

"... C01 or C45 should be ignored, because they definitely do not exist and do not participate in the calculation.”

That certainlly is NOT what the defining equations you posted imply.

It is unlear without any more context the source of the recursion relationships, but as @Torsten notes, if come from a set of PDEs, the BCs would set the values.

[Umar]

I’d like to offer two brief observations in light of the discussion above. In the formal derivation of tridiagonal systems—particularly those arising from PDE discretisations—every coefficient in the defining equations must be assigned a value, including those associated with the boundaries. The boundary conditions provide these values; if they are homogeneous, the corresponding couplings simply evaluate to zero.

From an implementation perspective, it is perfectly acceptable to omit terms for non-existent neighbours by treating their coefficients as zero. However, from a mathematical standpoint, defining them explicitly maintains clarity, preserves the structure of the recurrence relations, and makes later changes to the boundary conditions straightforward.

This supports the point raised that while certain couplings may appear “non-existent” in the physical mesh, they nonetheless have a defined mathematical role. Explicitly setting them—even to zero—keeps both the derivation and the computational implementation consistent.

Answer 1

[Umar]

Hi @Filip,

I’ve prepared the MATLAB script you asked for. It:

• Implements your provided formulas for V and s exactly.

• Works for any number of elements (just update n, a, D, and C).

• Builds the full V matrix and s vector step-by-step with printed intermediate results.

• Solves gamma = V \ s in a numerically stable way.

• Output format matches your example for easy verification.

This should fully address your requirement of generating the matrix for arbitrary size without utilizing matlab toolbox.

Script

% ---------------------------------------------------------
% Build V, s, and gamma for any number of elements n
% from given vectors a, D and full coupling matrix C.
%
% Implements:
%   v(i,i-1) = -C(i,i-1) * a(i-1)
%   v(i,i)   = (C(i,i-1) + C(i,i+1) + D(i)) * a(i)
%   v(i,i+1) = -C(i,i+1) * a(i+1)
%
%   s(i) = -C(i,i+1)*(a(i+1)-a(i)) + C(i,i-1)*(a(i)-a(i-1))
%
% gamma = V \ s
% ---------------------------------------------------------

clear; clc;

% --- Example input (change as needed) ---
n = 5;

a = [1.0; 0.95; 1.05; 0.9; 1.1];
D = [0.2; 0.15; 0.1; 0.25; 0.3];

C = zeros(n); % full coupling matrix
C(1,2) = 1.2;
C(2,1) = 1.2; C(2,3) = 0.8;
C(3,2) = 0.8; C(3,4) = 1.0;
C(4,3) = 1.0; C(4,5) = 0.7;
C(5,4) = 0.7;

% --- Allocate ---
V = zeros(n);
s = zeros(n,1);

fprintf('--- Building V and s ---\n\n');

for i = 1:n
  % Get neighbor couplings (zero if none)
  C_left  = 0; if i > 1, C_left  = C(i, i-1); end
  C_right = 0; if i < n, C_right = C(i, i+1); end

    % Fill V
    if i > 1
        V(i, i-1) = -C_left * a(i-1);
    end
    V(i, i) = (C_left + C_right + D(i)) * a(i);
    if i < n
        V(i, i+1) = -C_right * a(i+1);
    end

    % Fill s
    term1 = 0; if i < n, term1 = -C_right * (a(i+1) - a(i)); end
    term2 = 0; if i > 1, term2 =  C_left  * (a(i) - a(i-1)); end
    s(i) = term1 + term2;

    % Debug print
    fprintf('Row %d:\n', i);
    fprintf('  C_left  = %g\n', C_left);
    fprintf('  C_right = %g\n', C_right);
    if i > 1, fprintf('  V(%d,%d) = %g\n', i, i-1, V(i,i-1)); end
    fprintf('  V(%d,%d) = %g\n', i, i, V(i,i));
    if i < n, fprintf('  V(%d,%d) = %g\n', i, i+1, V(i,i+1)); end
    fprintf('  s(%d) = %g\n\n', i, s(i));
  end

% --- Final results ---
disp('Matrix V:'); disp(V);
disp('Vector s:'); disp(s);

gamma = V \ s;
disp('Gamma:'); disp(gamma);

Please see attached.

Hope this is what you’re looking for.

[Walter Roberson]

@Umar

Unfortunately, the two files that you attached are empty.

This is separate from the fact that you used <> file inclusion on a mobile browser; even when I correct for that, the retrieved files are empty.

[Umar]

Hi @Walter,

That is really weird, thanks for letting me know. I do appreciate your help and feedback. Please see attached.

Answer 2

[Torsten]

You can code it with V and C being sparse. For simplicity, I used full matrices.

n = 50;
V = zeros(n);
s = zeros(n,1);
V(1,1) = (C(1,2) + D(1)) * a(1);
V(1,2) = -C(1,2) * a(2);
s(1) = -C(1,2) * (a(2) - a(1));
for i = 2:n-1
  V(i,i-1) = -C(i-1,i) * a(i-1);
  V(i,i) = (C(i-1,i) + C(i,i+1) + D(i)) * a(i);
  V(i,i+1) = -C(i,i+1) * a(i+1);
  s(i) = -C(i,i+1) * (a(i+1) - a(i)) + C(i-1,i) * (a(i) - a(i-1));
end
V(n,n-1) = -C(n-1,n) * a(n-1);
V(n,n) = (C(n-1,n) + D(n)) * a(n);
s(n) = C(n-1,n) * (a(n) - a(n-1));
gamma = V\s