Jtdcftul

I have the following code:

max = 1;



PotentialTilde[V0_] := 

  SparseArray[{Band[{1, 1}] -> V0/1, Band[{2, 1}] -> V0/2, 

    Band[{1, 2}] -> V0/2}, {2*max + 1, 2*max + 1 }];



a = KroneckerProduct[PotentialTilde[1], PotentialTilde[1]] // 

  MatrixForm

That generates this matrix:

$$ $$">

How do I set every element that is not equal to 1 or 1/2 to 0?

Be careful with assigning a "form" (e.g. MatrixForm) to a variable. Another observation is that a is a SparseArray uses rules to store the values and since ReplaceAll works with the FullForm of an expression care has to be taken (using Normal will make the array a regular matrix again).

This should work:

max = 1;

PotentialTilde[V0_] := SparseArray[{Band[{1, 1}] -> V0/1, Band[{2, 1}] -> V0/2, 

    Band[{1, 2}] -> V0/2}, {2*max + 1, 2*max + 1}];



a = KroneckerProduct[PotentialTilde[1], PotentialTilde[1]];

(* a // MatrixForm *)



aMod = (a // Normal) /. x_?NumericQ /; Not@MatchQ[x, 1 | 1/2] -> 0;

aMod // MatrixForm

Lukas Lang's suggestion is even nicer:

aMod = (a // Normal) /. Except[1 | 1/2, _?NumericQ] -> 0;

UPDATE

@HenrikSchumacher is of course right in warning against the use of Normal for a SparseArray -- after all there may have been a good reason for using it in the first place. So for a pattern replacement solution the use of ArrayRules may be safer:

aMod = ( a // ArrayRules ) // RightComposition[

    ReplaceAll @ Rule[ 

        Rule[{row_?NumericQ, col_?NumericQ}, Except[ 1 | 1/2, _?NumericQ ],

        Rule[{row,col}, 0]

    ]

    , SparseArray

];

aMod // MatrixForm

A = KroneckerProduct[PotentialTilde[1], PotentialTilde[1]];

B = A (1 - Unitize[A - 1/2] Unitize[A - 1]);

B // MatrixForm

It is even much more efficient to use machine precision numbers:

max = 100;

A = KroneckerProduct[PotentialTilde[1], PotentialTilde[1]];

B = A (1 - Unitize[A - 1/2] Unitize[A - 1]); // RepeatedTiming // First



NA = N[A];

NB1 = NA (1. - Unitize[NA - 0.5] Unitize[NA - 1.]); // RepeatedTiming // First



Max[Abs[B - NB]]

0.269

0.0163

0.

It is even faster (but also quite a bit more tedious, though), to operate only on the list of nonzero values by using the low-level representation of SparseArray and the following undocumented constructor:

NB2 = SparseArray[

     SparseArray @@ {

       Automatic,

       Dimensions[NA],

       NA["Background"],

       {1(*version of SparseArray implementation*),

        {

         NA["RowPointers"],

         NA["ColumnIndices"]

         },

        With[{vals = NA["NonzeroValues"]},

         vals Subtract[1., Unitize[vals - 0.5] Unitize[vals - 1.]]

         ]

        }

       }

     ]; // RepeatedTiming // First

0.00670

The advantage over modifying ArrayRules is that it leaves all array packed (they can only be backed for machine precisision numbers or machine integers; not for symbolic entries) while ArrayRules unpacks the list of nonzero positions and nonzero positions in order to generate the symbolic list of rules.

Should you wonder why SparseArray appears twice here: SparseArray@@{...} constructs the sparse array literally by the row pointers, column indices and "nonzero values" provided, even if there are actually zeroes among the "nonzero values"; wrapping this with a further SparseArray removes these "ghosts" from the internal representation, so that NB2["NonzeroValues"] returns really only the nonzero values.

a cannot be a MatrixForm:

(a = KroneckerProduct[PotentialTilde[1], 

    PotentialTilde[1]]) // MatrixForm

Convert to a regular Matrix:

Normal[a] /. {ij_ /; (NumberQ[

       ij] && ( ij != 1/2 && ij != 1)) :> 0}



(a = Normal[

     a] /. {ij_ /; (NumberQ[

          ij] && ( ij != 1/2 && ij != 1)) :> 

      0}) // MatrixForm

a



MatrixForm[a]