MathCAD Prime 9.0 "instabiities"

[FB_10795922]

Dear all, I’m experimenting since a while, with MathCAD Prime 9.0 and a specific spreadsheet, what I’d classify as “numerical instabilities”.

I’ll explain better.

I’m dealing with a complex spreadsheet; this last:

1. Reads from an Excel file the input data
2. Recalls (include command) a spreadsheet defining a function of 6 variables  

To be noticed that the above points have not created any trouble in previous spreadsheets

      3. Involves recursive formulas with matrix multiplications and inversion

     4. Involves the solution of a transcendental equation deriving from the determinant of a matrix; this last             matrix is obtained at the end of the recursive process (hinted before

    5. Is calculated using the options:

•  Calculation/Calculation Options/ Multithreading;
• TOL & CTOL @ 10-7

The spreadsheet, in itself, is not that big (around 350 Kb).

Many different issues are highlighted.

First, with plotting (preliminary to roots search).

The plotting of the transcendental equation did not work. I created instead two vectors (see Fig.1 in the attached report). It worked once (and I could see the range where the first root is); but re-opening the worksheet the vectors were recalculated and ended in errors (see Fig.1 in the attached report)

Second, with root search.

The MathCAD find and root algorithms are not working at all. I could find (chrome-extension://efaidnbmnnnibpcajpcglclefindmkaj/https://www.msme.us/2012-1-2.pdf) an algorithm “intermittently” working….I mean, the algorithm converges to a solution (reasonably close to the expected one). BUT re-opening the spreadsheet, the root is often recalculated (!); and, sometimes, a solution is not even found (!). In this last case, I re-launch the calculation changing the range for root search (1 Hz in Fig.2 of the attached report ) or simply re-writing the values (in this case 300Hz or 1Hz – Fig.2)

Third, with “contradictory” reports:

See Fig.3 of the attached report

In the Fig.3, on the left hand side, you can notice that the function e0707 is reporting an error (when calculated with Sol1_mR value); but the same calculation gives a result, on the right hand side…on the right hand side, instead, the function e0708 reports an issue….

But BEFORE (in the spreadsheet), the same functions are calculated without any problem (?!?) – see Fig.4 in the attached report.

Is anyone ever experienced this kind of issues with MathCAD (Prime 9.0 variant in particular)?

Some suggestions to alleviate these issues and/or explanations of what is happening?

Personally, I feel the issues are a mix of inexperience with some algorithms (and not “lean” spreadsheet, maybe…) and real tough calculations.

Obviously I’m available for any further question.

I’ll try to attach to attach also the spreadsheet and the needed files if really necessary…

Many thanks and best regards

[LucMeekes]

A picture may paint a thousand words. A Prime worksheet delivers ALL information.

If possible, attach the worksheet(s) and the Excel sheets needed.

Success!
L:uc

[FB_10795922]

Good afternoon, I made a test with Werner_E and it seems to work: here attached the spreadsheet and the needed (include + input) files.

Thanks and best regards

[Werner_E]

Ad Fig. 1
Its can't be seen how your function "S13" is defined, so nothing can be said about the error you experience. Is S13 really creating a square matrix based on a single frequency? What is the result of, lets say, S13(55 Hz)=...

Ad Fig. 2
Seems to be based on the same problem with function S13, but we have not enough information to say anything further

Ad Fig 3

Problem may be that you try to numerically evaluate (the equal sign at the end) a function definition. Delete the equal sign in the definition of function C7b and see if it works. Example:

Error message would be the same, no matter if a variable x is defined above or not. The error message is misleading as x is the formal function argument and it does not matter if a variable with the same name is defined or not.

As we don't see the definition of the function e0708 nor the error message thrown, nothing can be said according the error when calling this function.

Ad Fig. 4

Kind of "Yesterday the vase was still intact and standing on the ledge and today it is lying on the floor in pieces."
Obviously something has happened. We can't know what it may have been, but sure something is different now compared to the last time when the file worked OK. An automatic Windows update changing the behaviour of Prime cannot be ruled out completely, but in my opinion it is highly unlikely.  An (intentional or not) change to the Excel input sheets or in the Prime sheet seems much more likely to me.
A small change can have a big effect (such as the unintentional addition of an equals sign after a function definition, as described in #3).

Generally when debugging a worksheet you should concentrate just on the first error in the sheet as subsequent errors may be just due to it. Just go further down in the sheet after you fixed the first error.

I fully agree with Luc that whenever possible you would have to provide the worksheet and the necessary Excel sheet(s) for more help.

[FB_10795922]

Good afternoon Werner_E, i'll try to attach the files BUT I can give you already some partial answers

  S13 matrix definition:

Please have a look to S13matrix.jpg

You can see that:

1) it is recursively defined (it depends from a S12 matrix. The recursion goes back until a matrix S0 which is, by definition, a 4x4 null matrix) 

2) it depends from 4 matrices (F11, F12, F21, F22) that are submatrix of the matrix function called (include) by the spreadsheet. The original variables (6) are fixed (from data recovered by Excel) except for the frequency (omega small). On the matrices F I have very few doubts because used tenths of times and correlated successfully with data retrieved in literature.....

You can also see a calculation check (highlighted in pale blue) for the root Sol1_mR. 

Ad Fig 3

Problem may be that you try to numerically evaluate (the equal sign at the end) a function definition. Delete the equal sign in the definition of function C7b and see if it works. Example:

Just after having finished the post I simply wrote again C7b(Sol1_mR) , and "miracle": now all is ok...(see attached .jpg)...and no, apparently the = after the definition is not disturbing.... 

I've tried to attach all the file..hope it works...

Thanks & Best Regards

[Werner_E]

I was talking about the equal sign after the DEFINITION of the function, not the one after the definition of the variable D7 using a function CALL!

In your new screenshot this equal sign is not present anymore and so it works OK

Do I understand correctly that the rest of the worksheet is now OK, too?

Would be glad to hear that because even the calculation of one single matrix S13 took quite long, calculation of 601 similar matrices and their determinant for the vector X would take shear endless.

At the moment I stopped calculation and closed the worksheet and would try to open it again only if you state that there still are problems.

[FB_10795922]

Dear WERNER_E, 

other "instabilities" appeared (and some solved...) to be honest. I'm doing a step forward that seems to work....if so, I could make the spreadsheet "leaner" indeed (I have the habit to make some intermediate calculation to understand problems in advance). BUT  for my problems a spreadsheet running during all night long is not so uncommon.... If you've some general advices to tackle this kind of calculations and make them leaner and more stable (root search in particular) they are more than welcome. Best Regards

[Werner_E]

Its pretty much possible to speed up the calculations you do but because of the nested nature of your functions it would take a lot of time to to get through it and unravel the calculations. Not sure if I will find time later to give it a try.

For Prime 11, it has been announced that the development focus would be on increasing the often criticized processing speed.

I have now determined the calculation time of the single matrix S13(55 Hz) in Prime 9 and Prime 11 and Prime 11 took more than twice as long! I am on a quite slow device - time in P9 was about 47 seconds, P11 took about 104 seconds. That was extremely disappointing.
It doesn't necessarily say anything about the quality of and improvements in Prime 11, but for the calculations in your worksheet, it obviously wouldn't be worth upgrading.

[FB_10795922]

Thanks WERNER_E,

this information is precious for me. In fact, as you may guess, all the questions I put in the forum are linked among them. The link is (in a synthetic manner):

Transfer Matrix Method  for chained, beam-type structures --> ok for low natural frequencies calculations ; issues at high natural frequencies and modes calculations (the issues are: round off type liked with many matrices multiplication AND hyperbolic functions entering the terms of each matrix; more the each matrix is not well conditioned BUT you've limited field of action due to the fact that each matrix contains terms with their own unit of measure - paramount for the method) 

Possible solutions --> the recursive method you can see in my worksheet; introducing an "ad hoc" unit of measurement system (allowing to making more "homogeneous" in terms of order of magnitude each matrix constituting the transfer matrix). In the literature the first method is boasted as "the solution" (old US articles and more recent Chineses articles). But my numerical test highlighted difficulties with the natural frequencies search already at relatively low natural frequencies....the second method is only hinted in literature and it's an idea of mine but really worth to be tested.

So for the first method, maybe (in view also of your further answer here below) Prime11 is not so convenient BUT it could be worth to test the 2nd method! Now the difficulty is that it seems that I'm allowed to upgrade from 9 to 11 free of charge. BUT a Prime 11 license would be needed to have both Prime 9 and Prime 11 on your PC.....

[Werner_E]

OK, I had a closer look at your sheet now and was able to speed up the calculation of S13 by a factor larger than 600.

All that has to be done is to avoid that S13 calculates S12 twice and that S12 calculates S11 twice, ...etc.

You may use a two line program (first lines calculates S12 which is then used in the actual calculation in the second line.

Or you may decide to use an inline assignment to a local variable XX:

Using recursion we could get rid of all the function S1  to S13 and replace them for a function Snn which, additional to the frequency argument has an argument referring to the number of stage (1 .. 13).

Now that calculation speed is increased it sure makes more fun to do the plotting or play around with the root function.

While looking at the plot I wondered why the max value did not show up in plot and found that there seems to be some kind of discrepancy for the position around 247 Hz. You may want to further investigate where this stems from and if there are other outliers.

According the root-finding I don't know why both flavors of the root function fail

Not sure if an effect similar to that at 247 Hz might be the cause. The implemented algorithm tries different values of omega (not necessarily near 300 Hz) and the error message may indicate that there is at least one frequency for which S13 is returning a matrix with wrong units. I tried a few values but got a correct result every time. I may be needed to look in more detail in the defining functions but thats something I leave up to you 😉

The solve block with find also fails, but with the error "No solution found". The reason is that the magnitude of the function values is so large that tiny arguments changes effect large changes in function values and so Prime could not find a solution within the tolerance.

One way out can be to simply decrease the function values by dividing them by a large number (at least 10^4).

As you can see the solve block works OK now, but the last line shows that the result is a bit less accurate than the one you got by your calculation on the left side.

Prime 9 sheet with all the changes is attached

[Werner_E · Accepted solution]

I still don't know which frequency results in a wrong matrix S13 and is responsible for the failing of the root function.

But we can wrap the function in a try...on error.. envelope and now the first variant of the root function works OK and even with a higher accuracy than your manually derived value. But it needs much more time to find the solution compared to the solve block.

I don't understand why the second variant of the root function still fails.

It claims that the function values at the end points of the interval provided (299 Hz ... 301 Hz) are not of opposite sign which clearly is not the case.

I also tried to make the function completely unitless, but to no avail

It gets even more crazy - scaling the function, making only the argument unit-less and shifting it horizontally does the job. It does not matter in which direction the function is shifted (+1. -1).

I have no explanation for this weird behaviour!

[FB_10795922]

Dear WERNER_E,

again many thanks for all your indications: indeed I'm making a jump in MathCAD use!

About the problem you're hinting of: me also trying to find roots of the transcendental equation coming out from the recursive method I suspected

that the function is oscillating around the root with really small amplitudes...or oscillations are purely numerical: I must confess I've not all my ideas clear in my head but with transcendental equations coming from this kind of methods we can go very near to  maximum number of digits allowed by MathCAD...

[FB_10795922]

THANK YOU SO MUCH Werner_E!!

First of all, the algorithm you propose to get rid of recursion is something I was looking for since a while for many other applications beyond the specific spreadsheet you saw! 

Second: about the "ghost" root around 247 Hz (better, I use improperly Hz...as you've correctly understood it should be rad/sec) 

It's a MathCAD bug I've already seen : try to "zoom" the plot and you'll see the root disappear! In fact MathCAD is "reproducing" the asymptote of the function (247Hz is a pole of the characteristic equation det(S13(omega))=0).

Literature traces the fact that the characteristic equation coming out form this recursive method, typically presents asymptotic behaviors; they are quite disturbing for automatic search of roots.

I've retrieved an article suggesting a way to get rid of this issue. But, I must confess I've not all understood.

More, I feel it's more simple, for now, plotting and finding roots one by one (the final goal is a modal reconstruction that, normally, asks for a limited number of modes in my typical problems....).  

Third

You make me discover an option I did not know (second line of image below)....more, this is a real root very close to the expected result (around 330Hz....) 

Same for dividing with a certain amount: very good idea! I did not think we could implement so easily!

Fourth 

about: "the error message may indicate that there is at least one frequency for which S13 is returning a matrix with wrong units" I've already seen this type of message and not always a REAL problem of unit of measure....I suspect more un-ability to converge or overflow...but, as for your check, the unit of measure problem, is in my modest experience, a misleading indication, not unusual in MathCAD....

[Werner_E]

The "ghost root" at 247 Hz has nothing to do with the plot. When you evaluate values around this position you clearly see (its in the screen shot I posted) that we get a positive value for frequencies lower than 247 Hz and negative for higher frequencies.

That's what your function detS13 returns and its independent from the plot or a zoom you apply there. Sure it may be wrong values due to numerical inaccuracies or numbers going to the limit of the IEEE format used by Prime. Would be time consuming to find out, I guess.

But when zooming in makes this asymptote and sign change disappear, then its a problem with the plot - zooming in would also require to make the range used at the omega axis much denser and adapt the limit and the y axis appropriate. Primes plot facilities are very limited and a true zoom is not provided.

But with the now much speedier S13 calculation it may be a bit easier for you to further investigate.