Loading [MathJax]/jax/output/HTML-CSS/fonts/TeX/fontdata.js

Generating nice complex systems of linear equations


The easiest way to come up with a system of linear equations that requires partial pivoting to solve them is to:

  1. Start with a matrix that is already in row-echelon form. Note that I do not require the leading digit to be a '1'. I simply describe the requisite entry not as the "leading one" but rather the "leading non-zero entry."
  2. Starting with jn1 and going down to 1, in Column j, add appropriate multiples of Row j onto each subsequent row.
  3. Having a dense matrix, multiply it by a nice solution to get an appropriate target vector.

Nice multiples are those z such that |z|<1 where each component has at most one digit after the decimal point, thus including:

    0        0.1j      0.2j      0.3j      0.4j      0.5j      0.6j      0.7j      0.8j      0.9j
    0.1  0.1+0.1j  0.1+0.2j  0.1+0.3j  0.1+0.4j  0.1+0.5j  0.1+0.6j  0.1+0.7j  0.1+0.8j  0.1+0.9j
    0.2  0.2+0.1j  0.2+0.2j  0.2+0.3j  0.2+0.4j  0.2+0.5j  0.2+0.6j  0.2+0.7j  0.2+0.8j  0.2+0.9j
    0.3  0.3+0.1j  0.3+0.2j  0.3+0.3j  0.3+0.4j  0.3+0.5j  0.3+0.6j  0.3+0.7j  0.3+0.8j  0.3+0.9j
    0.4  0.4+0.1j  0.4+0.2j  0.4+0.3j  0.4+0.4j  0.4+0.5j  0.4+0.6j  0.4+0.7j  0.4+0.8j  0.4+0.9j
    0.5  0.5+0.1j  0.5+0.2j  0.5+0.3j  0.5+0.4j  0.5+0.5j  0.5+0.6j  0.5+0.7j  0.5+0.8j
    0.6  0.6+0.1j  0.6+0.2j  0.6+0.3j  0.6+0.4j  0.6+0.5j  0.6+0.6j  0.6+0.7j
    0.7  0.7+0.1j  0.7+0.2j  0.7+0.3j  0.7+0.4j  0.7+0.5j  0.7+0.6j  0.7+0.7j
    0.8  0.8+0.1j  0.8+0.2j  0.8+0.3j  0.8+0.4j  0.8+0.5j
    0.9  0.8+0.1j  0.8+0.2j  0.8+0.3j  0.8+0.4j

I would suggest the following code:

Get a multiplier that is of the form z=±0.m±0.nj such that |z|0.75. This second condition ensures that it is reasonably easy to determine the maximum entry in absolute value.

function z = multiplier()
  do
    z = (round( rand()*19 ) - 10)*0.1 + (round( rand()*19 ) - 10)*0.1*j;
  until abs( z ) < 0.75
end

Get an entry of the matrix that is non-zero on the diagaonal entries and possibly zero on the strictly upper-triangular entries.

function z = entry( nonzero )
  do
    z = round( rand()*19 - 9.5) + round( rand()*19 - 9.5 )*1j;
  until ~nonzero || abs( z ) > 0;
end

Code to find the matrices is as follows, but we check to make sure that, for example, at most one entry of the solution vector is zero, and that no component in the target vector has more than two digits in the significand.

N = 3;    % This can be any positive integer
A = zeros( N, N );
for i = 1:N
  for j = i:N
    A(i, j) = entry( i == j );
  end
end

for j = (N - 1):-1:1
  for i = (j + 1):N
    A(i,:) += multiplier()*A(j,:);
  end
end

do
  do 
    u = 0.1*round( 19*rand( N, 1 ) - 9.5 ) + 0.1j*round( 19*rand( N, 1 ) - 9.5 );
  until sum( u == 0 ) <= 1;
  b = A*u;
  
  found = true;
  arb = abs(real(round( 100*b )));
  aib = abs(imag(round( 100*b )));
  
  % Make sure no entry in the target vector has more than
  % two digits in the significand, so mn.00, m.n0 or 0.mn.
  if    any( arb >=  100 &  10*round( arb/10 ) ~= arb ) ...
     || any( aib >=  100 &  10*round( aib/10 ) ~= aib ) ...
     || any( arb >= 1000 & 100*round( arb/100 ) ~= arb ) ...
     || any( aib >= 1000 & 100*round( aib/100 ) ~= aib )
    found = false;
  end;
until found;
A
u
b

In order to introduce required row swaps, start with Row n1 going down to Row 1, and at each step, optionally swap Row i with any Row k where k>i:

for i = (N - 1):-1:1
  k = i + floor( rand()*(N - i + 1) );
  A([i,k],:) = A([k,i],:);
  b([i,k]) = b([k,i]);
end

Without row swaps, some examples for N=2 are:

(43j5+4j3.41.3j11.3+10.8j)(0.6+0.5j0.7+0.4j)=(1+j2.2+9.6j)

(12j69j1.30.9j2.9+7.3j)(0.2j0.6+0.6j)=(1.4+9.2j6.32.9j)

(4+8j2+2j2.8+0.4j1.20.6j)(0.9+0.9j0.40.4j)=(3.6+9.2j3.62.4j)

(6+5j64j4.1+3.7j1.2+13.6j)(0.2+0.6j0.10.2j)=(5.61.8j1.2+4.8j)

Note that because the multipliers are of the form 0.m+0.nj, it follows that all the ratios 3.41.3j43j=(3.41.3j)(4+3j)25=175j25=0.70.2j, 1.30.9j12j=(1.30.9j)(1+2j)5=0.5+3.5j5=0.1+0.7j, 2.8+0.4j4+8j=(2.8+0.4j)(48j)80=8+24j80=0.1+0.3j and 4.1+3.7j6+5j=(4.1+3.7j)(65j)61=6.142.7j61=0.10.7j must simplify to such a value.

Again, without row swaps, some examples for N=3 are:

(28j91j9+5j0.83.6j0.40.2j2.6+10.8j3.4j1.20.2j5.2+11.8j)(0.60.5j0.6+0.4j0.5j)=(1.92.3j6.41.3j6.85j)

(98j2+4j81j2.28.4j11.83.4j7.22.4j1.28.6j0.4+2.4j10+2.4j)(0.9+0.1j0.50.7j0.5+0.3j)=(1.22.8j5.10.26j8.12.5j)

(1+3j2+9j12j0.31.9j3.3+2.4j3.6+2.3j0.41.2j1.71.2j8.81.9j)(0.70.5j0.3+0.4j0.2+0.1j)=(1.8+1.2j2.2+0.06j1.60.5j)

(44j5+2j76j3.6+0.4j8.3+5j1.1+1.2j4+0.8j4.82.5j5.55.9j)(0.40.2j0.30.6j0.90.2j)=(44j1.6+6.9j0.615.4j)

Finally, again without row swaps, two examples from N=4 are:

(5+5j4j9+1j7+9j22j9.67j2.4+10.6j11.61.8j3.50.5j0.6+0.3j8.1+3.8j14.11.8j4.52.5j1+1.3j1.6+7.3j1.315.5j)(0.4+0.8j0.7+0.8j0.70.5j0.30.4j)=(3.1+1.5j0.7+0.66j4.61.8j0.18+0.15j)

(7+7j1+9j48+9j2.8+2.8j5.6+8.6j4.6+9j6.25.4j3.5+2.1j4.3+1j8.98.9j6.80.1j4.2+1.4j7.53.4j9.4+8.2j4.617j)(0.30.4+0.7j0.1+0.3j0.40.3j)=(2.55j61.9j1+0.5j3.70.89j)