Solving simultaneous equations with big matrices using Mathematica 7

[kaizen.moto]

Hi all,

Please help me to solve my problem in Mathematica that involves simultaneous equations with large matrices. The code is provided in the attached file, Exercise1.nb

the errors say like this " General::unfl: Underflow occurred in computation. >>

General::unfl: Underflow occurred in computation. >>

General::unfl: Underflow occurred in computation. >>

General::stop: Further output of General::unfl will be suppressed during this calculation. >>

indeterminate"

Note:

1. the constants that remains during the calculation are a, q and c11. In fact the final result that Iam expecting should be interms of a q / c11.

2. it is highly recommeded to use Mathematica version 7 rather than ver 8. I don't know the particular reasons of this. I have tried to use v8 but it really slow and eventually crashed.

3. running the code using v7 should take less than a minute.

many thanks for any feedback.

 

[Bill Simpson]

Near the bottom of your notebook in this
eqn = {Ao[1] + sumU1 == 0,...}

I am just guessing, but you have never defined Ao[] or initialized it to any value. This may be the result of your cutting up a notebook to supply this example or it may be a misunderstanding.

That appears like it might be the reason for your Underflow and possibly the first of layers of problems to uncover and correct. But it is not at all clear what you think you are trying to do with this. It appears you have 34 equations x+something=0, now solve for x.

It is also possible, since you are taking 10th powers of matricies and I already see lots of variable with very tiny values that you may just be driving some values so small that they are unrepresentable and this is the cause of your underflow.

I suppose you could try littering your code with Chop[] which will drive very small values to exactly zero and see whether this results in singular matricies which have no solution. Or maybe that will just get rid of your error message.

It is not specifically an error, but near the top you define and initialize D0 and E0, then manually copy element by element these into xx and yy, then deconstruct xx and yy with subscripts into D11...D66 and E11...E66 only then to reassemble all this back into matrices, or maybe they are functions, D00 and E00.

It seems like you could do this with far less torture by
D00[a_, λ_, c11_] := Evaluate[D0-λ*IdentityMatrix[6]];
E00[a_, μ_, c11_] := Evaluate[E0-μ*IdentityMatrix[6]];

And I am not even certain why you need the Evaluate in this specific situation, other than likely it didn't work and you tried sticking in an Evaluate to see if it fixed it.

Perhaps you have good reasons why you are making some of this as complicated as you have, but doing it the way you have seems like it gives you many more opportunities to make small errors and to have the notebook be far less understandable.

If you can find a way to write notebooks that reflect high level algebra matrix operations rather than the FORTRAN method of moving individual bytes then I think this may work out much better for you. Can you turn this 15 pages of subscripting into a stunningly clear 2 pages of matrix algebra where it would be far easier to understand exactly what you are really trying to do?

 

[kaizen.moto]

Thank you for the response.

I have attached the revised notebook and frankly speaking, I have no idea on how to turn my 15 pages notebook into a 2 pages as I consider my self as a novice user to Mathematica, let alone to any other programme languages. I have never used Fortran or any other program, Mathematica is my first program.

In this revised code, I have made some simplicfications especially in finding the eigenvalues.

It seems that syntax of sub1..sub2...sub17 are recursive process but I have no idea on how to make them into a simpler code.

The problem is yet again with regard to 'underflow'. Please help me to solve this issue.

Many thanks in advance.

 

[Bill Simpson]

After you have evaluated your notebook then look at the value of eqn.

You should see that it has vast numbers of extremely small complex values.

Then look at the value of Chop[eqn].

That should be dramatically smaller because it sets all extremely small values to zero and then many zeros simply disappear. Verify that this reduced set of equations is what you are expecting.

Then evaluate

sols = Solve[Chop[eqn], {Ao[1], Ao[2], Ao[3], Ao[4], Ao[5], Ao[6], Ao[7], Ao[8], Ao[9], Ao[10], Ao[11], Ao[12], Ao[13], Ao[14], Ao[15], Ao[16], Ao[17], Bo[1], Bo[2], Bo[3], Bo[4], Bo[5], Bo[6], Bo[7], Bo[8], Bo[9], Bo[10], Bo[11], Bo[12], Bo[13], Bo[14],
Bo[15], Bo[16], Bo[17]}]

There are still many underflows happening as it tries to solve this and I suspect those warnings actually mean that there are deep serious problems that you need to find and fix to get sensible and partly correct results out of this.

In this output you should see things like Ao[16]->0. a q/c11 which means that Ao[16] is approximately zero. Verify if those are correct in your problem or not. If they are then you can even try

sols/. 0.->0

and inspect the result. That will turn approximate zeros into exact zeros and show that about ten of your variables appear to be zero.

It appears to me, but I'm just guessing, is that the underflows and zeros mean that your system of equations is possibly singluar and it isn't finding a way to solve for Bo[15], Bo[16] and Bo[17]. Verify if that is correct or whether the solver simply failed part way through when the underflows caused the solver to fail.

 

[kaizen.moto]

After running the notebook, I agree that there are lots of smaller values of Ao and Bo interms of complex numbers. These complex numbers consists of real parts and imaginary parts. Imaginary parts are extremely small and I consider this imaginary parts as 'numerical noise' and therefore these can be neglected. I verify that these are what Iam expected.

After running with sols=Solve[Chop[eqn],{Ao[1]...Bo[17]}] for the case of m = 1 and n = 1. I still got these underflows issues. Iam expecting this is the case, all the unknowns are very close to zeros.

However, if I run the above for the case of m = 7 and n = 7, the problems of underflows disappeared. But, the values of each uknowns represent of another uknowns. For instance, Ao[1] -> 1.Ao[14], Ao[2] -> -1332.53 - 0.0530913 aq/c11 -0.000446269Ao[11]...

Iam guessing that the problem is due to the presence of aq/c11 in every equation that creating the simulataneous equations cannot be solved.

Is there a way in Mathematica to define in a such away that aq/c11 is a constant and will not be affected during the computation.

 

[kaizen.moto]

Iam just wondering if my assumptions are correct in such a way that I assume a=1,q=1 and c11=1. With these assumptions, I got the final results, i.e. all the unknowns are solved.

However, if I remove these assumptions, I got the same problems with regard to 'underflows'. These underflows have caused the solver to fail partway.

I appreciate for any feedbacks.

 

[aTn]

This usually occurs because you problem is badly scaled. In broad terms: it is ill-conditioned because the eigenvalues of your matrices have differend orders of magnitude. Try to define scaled variables by dividing your original variables by (non-zero) constants (or via linear tranformations) so that the new variables all lie in a range close to 1.

Here's an example.

Say you want to solve a linear system where x and y are your unkowns.

Let's assume that x takes values between 10^4 and 10^6 and that y takes values between 0 and 1.

You might want to define a scaled variable (let's denote it by x_s) via x_s = x/10^6. The new variable x_s will vary between 0.01 and 1.

Modify your equations to state them in terms of x_s and y. Try to solve the new system of equations with respect to the scaled variables.

I highly suggest you read on the subject of ill-conditionned linear systems since these things will pop up in all sorts numerical calculations, including numerical optimization (ill-conditioned Jacobians or Hessians), solving systems of nonlinear equations, solving differential equations, etc.

 

[kaizen.moto]

I have tried to scaled out the coefficients of the variables but none of them were successful. The problem is that I have 34 unknowns to be scaled over with 34 different equations. I have no idea on how to improve the equations so that they would become well-conditioned matrix. I have had also look for other references with various solving techniques but still no luck.

After running the code, I got solutions from m = 1 and n = 1 until m = 7 and n = 7. These are well conditioned matrices, I think. However the problems started from m = 9 and n = 9 where ill-conditioned started to occur.

Any response is very much appreciated

 

[Bill Simpson]

Just as a couple of experiments, try this

1: Replace your machine precision decimal approximations with some exact rational numbers that are close in value to the constants you are using.

2: Replace your machine precision decimal approximations with some say 50 digit precision numbers. And if that doesn't work then try 100 digits and repeat as necessary.

Neither of those would be approved by someone doing careful numerical analysis of every step in your problem, but it doesn't sound like that is being done anyway.

Do either of those substantially change what you are seeing?

If not then delete a BUNCH of your semicolons at the ends of lines that are hiding the results of each of your steps and carefully compare the m,n=7,7 with the m,n=9,9 and track down exactly which step is blowing up.

Perhaps with one of these someone can help point out what to try next

 

[kaizen.moto]

Thank you for your feedback.
Actually, I have no idea about the first method but I have done the second option by setting the precision of the machine to say 50 to 100^100. I did this by using the command such as Set Precision, Working Precision, etc. However, no changes are affected.

I have also tried to change the settings from edit -> preferece, but I don't know what to change from this method.
Please let me know if there any other methods about increasing the machine precision.

I have also tracked down the steps by removing the semicolons. the problems started from solving the simultaeous equations, i.e. result for RowReduce of badly conditioned matrix. I suspected this errors come from the order of coefficient, meaning high different magnitudes of constants coefficients among the 34 unknowns.

 

[Bill Simpson]

I am not certain I know how you have set your precision.

Example:

In[1]:= x=3.14158;Precision[x]
Out[2]= MachinePrecision

In[3]:= N[MachinePrecision]
Out[3]= 15.9546

In[4]:= x=3.14158000000000000000000000000000;Precision[x]
Out[5]= 32.4971

Can you see from that at least one way to set and check precision of constants?

Can you verify that the variables at the top of you notebook actually have 50 digits of precision?

Can you check the implementation details for the functions you are using to solve the simultaneous equations to see if they respect higher precision values or if they force everything to MachinePrecision? That is buried in the documentation somewhere. You want to make certain that the method you are using is not ignoring your extra precision. If it is not then you should be able to extend the precision far enough to overcome the failure.

"Elementary Numerical Computing With Mathematica" by Skeel and Keiper is very old but I find it an excellent resource when I am trying to bend Mathematica to solve a difficult precision problem. Used copies are sometimes cheap, but not so much today.

http://www.bestwebbuys.com/Elementa...-Mathematica-ISBN-9781588740533?isrc=b-search

 

[kaizen.moto]

Thanks for your feedback.
After running the codes with SetPrecision[....,50], I verify that the final output gives only up to 15.9546 precision. I mean, Mathematica forced everything to the default MachinePreicision.

I have got your recommeded book with me now by the way, its very good source of information about Mathematica for beginners like me.

Can you please let me know how to set Mathematica so that the code would run with 50 digits precision in everything?

Thanks for your kind help.

 

[Bill Simpson]

Ah, my guess about what you had done seems to have been correct.

SetPrecision is not an incantation telling Mathematica to change its mind.
SetPrecision is a function, like Sin, that returns a value with a new precision.

Either
C12 = 0.246269000000000000000000000000000000000000000000*c11
or
C12 = SetPrecision[0.246269, 50]*c11
should give you C12 with 50 digits of precision. (Actually the first one has only 48 digits because when I use 50 the number is displayed here with a blank inserted before the last two zeros, count your own digits in your notebook to verify correctness)

The first will ensure trailing digits are zero, the second will ensure trailing bits are zero. Print C12 after after each of these to see the result.

Every time you enter any numeric constant with only a few digits that will be interpreted as a MachinePrecision approximate number. Every time you multiply or add any such constant to any other expression the result will become a MachinePrecision approximate number and all the extra precision you might have had will be gone.

So check the Precision of all your variables to make certain you really have 50 digits before you start doing your matrix calculations.

 

[kaizen.moto]

Great..its working fine now...

Thank you so much..