0

I am using Matlab's backslash operator to solve a system of equations written as two matrices M1 and M2. These two matrices are square and tridiagonal, and so I have defined them as sparse. For example, with the dimensions of each being 5x5, they are defined as follows, with the values in each entry being dependent on some constant a:

N = 5;
a =  1e10;
M1 = spdiags([-a*ones(N,1)...        % Sub diagonal
              (1 + 2*a)*ones(N,1)... % Main Diagonal
              -a*ones(N,1)],...      % Super diagonal
              -1:1,N,N);
M2 = spdiags([+a*ones(N,1)...
              (1 - 2*a)*ones(N,1)...
              +a*ones(N,1)],...
              -1:1,N,N);
M_out = M1\M2;

So for example, M1 looks like the following in full form:

>> full(M1)

ans =

   1.0e+10 *

    2.0000   -1.0000         0         0         0
   -1.0000    2.0000   -1.0000         0         0
         0   -1.0000    2.0000   -1.0000         0
         0         0   -1.0000    2.0000   -1.0000
         0         0         0   -1.0000    2.0000

Now, if I examine the number of non-zero entries in the result M_out, then I can see they are all non-zero, which is fine:

>> nnz(M_out)

ans =

    25

The problem is that I also need to do this for larger values of the constant a. However, if, for example, a=1e16 instead, then the off-diagonal entries of M_out are automatically set to zero, presumably because they have become too small:

>> nnz(M_out)

ans =

    5

Is there a better way in Matlab of going about this problem of inverting sparse matrices? Or am I using the backslash operator in the wrong way?

teeeeee
  • 641
  • 5
  • 15

1 Answers1

2

If the size of your matrices doesn't grow too much, I recommend doing a full symbolic computation:

N = 5;
syms a
M1 = diag(-a*ones(N-1,1),-1) + diag((1 + 2*a)*ones(N,1),0) + diag(-a*ones(N-1,1),+1);
M2 = diag(+a*ones(N-1,1),-1) + diag((1 - 2*a)*ones(N,1),0) + diag(+a*ones(N-1,1),+1);
M_out = M1\M2;

M_num_1e10 = subs(M_out,a,1e10);
M_num_1e16 = subs(M_out,a,1e16);

vpa(M_num_1e10)
vpa(M_num_1e16)

In that case, you will need the Symbolic Math Toolbox. If you don't have it, I think you should considerer migrating to Python and work with SymPy.

EDIT:

Considering the way you defined your problem, you need extended precision for your computations. The double precision isn't enough. For example, in double precision (1e16+1) has to be rounded to (1e16), in other words (1e16+1)-(1e16) is equal to zero. So your problem starts in the main diagonal of your matrices. MATLAB only provides extended precision through its symbolic toolbox.

If you want to stick with double precision, you may extend the double precision yourself relying on the so called double-double arithmetic. I say that you will have to do it by yourself because I don't think there is a open source double-double library for MATLAB.

Zalnd
  • 207
  • 1
  • 8
  • I don't have access to the toolbox. Your answer is not particularly helpful - both your suggestions involve acquiring a new piece of software. I am looking for a way to obviate this problem in core *Matlab*, which I'm sure must be capable of handling this situation. – teeeeee Jan 13 '21 at 10:09
  • 1
    Python is open source, that's exactly why I suggested it as the second option. About the statement "I am looking for a way to obviate this problem in core Matlab, which I'm sure must be capable of handling this situation", well I strongly recommend you to read more about numerical analysis and computer number formats... – Zalnd Jan 13 '21 at 12:37
  • 1) This is part of a much larger set of codes, and I am not prepared to migrate thousands of lines of cpde over to a different language unless there is a very good reason (not simply because it is your preferred language). 2) Ofcourse I can spend the next six weeks taking a textbook from the libary and studying in detail all about numerical analysis and computer number formats, and will eventually solve this problem myself. However, the exact reason for this site existing is so that people with existing knowledge in particular areas can share it to help others in their problems. – teeeeee Jan 13 '21 at 12:48
  • 1
    But it seems you are a specialist in both topics. Let me examine your behavior. Someone answered your question. You replied that the answer wasn't helpful because you were not interested in buying anything nor prepared to migrate thousands of lines of code over to a different language not simply because it is the preferred language of the person you were arguing with. You were uncapable of thinking: "wait a minute, maybe the person that is trying to help me out know there is no magical solution to my problem and pointed out the best solution to the best of his knowledge". – Zalnd Jan 13 '21 at 13:20
  • Please don't get me wrong - I really appreciate you taking to time to offer an answer which I'm sure you really do believe provides a good solution to the best of your knowledge, so thank you for doing so. However, my opinion that your answer is not particularly helpful to me still stands. Don't take the feedback I gave to your answer as personal. The hope is that feedback to answers allows those answers to be refinied, and to eventually converge on something which the OP feels answers his/her question satisfactorily. – teeeeee Jan 13 '21 at 13:39
  • I restate, there is no magical solution for your problem. I will end my contribution here, so good luck. – Zalnd Jan 13 '21 at 14:32