15c challenge: accuracy of a complex equation system

[J-F Garnier]

I recently proposed a simpler solution of a problem example from the 15c Advanced Function Handbook (AFH),
solving a complex 4x4 system that represents an electric circuit driven by an AC source:
   
(See the AHF page 128 for details)

The problem can be expressed numerically by the complex matrix equation:

| (100,-350/3) (0,800/3)    (0,0)       (0,0)       |   | I1 |   | 10 |
| (0,800/3)    (1E6,-350/3) (-1E6,0)    (0,0)       |   | I2 |   | 0  |
| (0,0)        (-1E6,0)     (1E6,442/3) (0,-150)    | x | I3 | = | 0  |
| (0,0)        (0,0)        (0,-150)    (1E5,442/3) |   | I4 |   | 0  |

My solution is using a direct system solving (matrix "division") and is running only on an extended memory 15C such as the CE.
I was happy to write:

(07-29-2023 09:15 AM)J-F Garnier Wrote:  The results are in line with the AFH results obtained with the alternate method, with a minor difference in the output phase.

However, later on I wasn't satisfied by this minor difference that was after all not so minor since in the 5th place of the result.
So I run again both solutions and recorded the full results (rectangular forms):

AFH's solution:
I1 =  1.995031379 e-4  ,  4.096402192 e-3
I2 = -1.448873903 e-3  , -3.563298144 e-2
I3 = -1.454123461 e-3  , -3.563275918 e-2
I4 =  5.344584749 e-5  , -2.259869709 e-6

My solution:
I1 =  1.995210750 e-4  ,  4.096401801 e-3
I2 = -1.448866059 e-3  , -3.563300697 e-2
I3 = -1.454115617 e-3  , -3.563278473 e-2
I4 =  5.344584748 e-5  , -2.259916891 e-6

The results are indeed a bit different. To find out which is better, I compared with Free42:

Solution on Free42 (rounded to 10 digits):
I1 =  1.995795134 e-4  ,  4.096399075 e-3
I2 = -1.448833616 e-3  , -3.563298308 e-2
I3 = -1.454083170 e-3  , -3.563276083 e-2
I4 =  5.344581171 e-5  , -2.259868251 e-6

I was expecting my results to be more accurate but to my surprise, both 15C results are inaccurate !


So the challenge is:

Improve the solution on the 15c to get 9-10 correct digits for each term,
that is with an error of a few ULP on each term.
Solutions for the 41C (for instance with the Advantage ROM) are welcome too.

I will post my simple 15c solution and comments in a few days.

Note: the Free42 reference solution is obtained by using exactly the same system matrix with 10 significant digits, by doing SCI 9 RND on the system matrix before solving it.

J-F

[Werner]

There's a typo, JF: the two instances of 350/3 must actually be negative.
Cheers, Werner

[J-F Garnier]

Many thanks Werner ! It's corrected now.

J-F

[Werner]

(08-11-2023 07:56 AM)J-F Garnier Wrote:  I was expecting my results to be more accurate but to my surprise, both 15C results are inaccurate !

While in general a direct solve is somewhat more accurate (and 3x faster or so) than a multiplication by the inverse, they will both lose 4-5 digits due to the problem's condition number of ~2e4

Quote:Solutions for the 41C (for instance with the Advantage ROM) are welcome too.[/list]

not for the 42S?

Cheers, Werner

[J-F Garnier]

(08-11-2023 10:47 AM)Werner Wrote:  

Quote:Solutions for the 41C (for instance with the Advantage ROM) are welcome too.[/list]

not for the 42S?

Solutions for the 42S (or other Saturn machines) are not forbidden, just they won't be directly comparable to the 15c, due to the Saturn 12-digit accuracy.
You should then target 11-12 correct digits, and compare with Free42 rounded to 12 digits too. Please keep the 10-digit input values of the system matrix for rough comparison with the 15c.

J-F
(edit for clarification)

[John Keith]

(08-11-2023 07:56 AM)J-F Garnier Wrote:  Solution on Free42 (rounded to 10 digits):
I1 =  1.995795134 e-4  ,  4.096399075 e-3
I2 = -1.448833616 e-3  , -3.563298308 e-2
I3 = -1.454083170 e-3  , -3.563276083 e-2
I4 =  5.344581171 e-5  , -2.259868251 e-6



...

Note: the Free42 reference solution is obtained by using exactly the same system matrix with 10 significant digits, by doing SCI 9 RND on the system matrix before solving it.

The rounding itself introduces inaccuracies, although I understand the reasoning behind it. Since the values in the example are exact, the system can be solved using exact rational arithmetic. In case anyone is interested, here are the exact results rounded to 10 digits.

Code:


 ( 1.995795141E-4  ,  4.096399085E-3)   
 (-1.448833619E-3  , -3.563298308E-2)   
 (-1.454083174E-3  , -3.563276083E-2)   
 ( 5.344581171E-5  , -2.259868256E-6)

The Free42 results seem to be off by no more than 10 ULP's, not significant for most real-world uses.

[J-F Garnier]

(08-11-2023 12:37 PM)John Keith Wrote:  Since the values in the example are exact, the system can be solved using exact rational arithmetic. In case anyone is interested, here are the exact results rounded to 10 digits.

Code:


 ( 1.995795141E-4  ,  4.096399085E-3 ) 
 (-1.448833619E-3  ,  -3.563298308E-2) 
 (-1.454083174E-3  ,   -.0356327608  )   
 ( 5.344581171E-5  , -2.259868256E-6 )

The Free42 results seem to be off by no more than 10 ULP's, not significant for most real-world uses.

I'm sure you understood it, but the Free42 results are not in error, they are results of a different problem.
The problem is dealing with an input system matrix with 10-digit values.
To be even more clear, that system matrix:

| (100,-116.6666667) (0,266.6666667)    (0,0)             (0,0)              |
| (0,266.6666667)    (1E6,-116.6666667) (-1E6,0)          (0,0)              |
| (0,0)              (-1E6,0)           (1E6,147.3333333) (0,-150)           |
| (0,0)              (0,0)              (0,-150)          (1E5,147.3333333)  |
assuming all the coefficients are exact.

You may try to calculate the exact solution of that problem, although the goal of my 'challenge' is to get 9-10 correct digits on the 15c or 41C.

J-F

[Werner]

mh a residual correction step still does not get 9-10 digits correct (nor a second one)

the following program will perform the solve and a residual correction, and allows for subsequent corrections:
You will need to be in big mode to run this program.

In: X:A Y:b
GSB A -> result C of solve + 1 residual correction
GSB C -> perform an additional residual correction

remark that the original matrix A is now in D, should you wish to recover it.

LBL B
STO MATRIX D
RESULT C
/
LBL C
RCL MATRIX B
STO MATRIX E
RESULT E
RCL MATRIX D
RCL MATRIX C
MATRIX 6
RCL MATRIX A
/
RCL MATRIX C
RESULT C
+
RTN

after first step:

1.995794978e-4 4.096399076e-3
-1.448833623e-3 -3.563298309e-2
-1.454083177e-3 -3.563276084e-2
5.344581171e-5 -2.259868262e-6

after GSB C:

1.995795337e-4
-1.448833607e-3
..

element C(2,1) has actually gotten worse, and we're nowhere near 9/10 digits correct.
Cheers, Werner

[J-F Garnier]

Thanks John and Werner for your interest.

I was expecting other usual math fans here to quickly find a solution (if not mine).
Remember I'm not a mathematician. My solution is what an average+ engineer can find.

I will post my (again, simple) solution on Wednesday, so last chance for your tentative!

J-F

[Albert Chan]

This just setup matrix A and B, 10-digits precision

10 DESTROY ALL @ OPTION BASE 1
20 R=4 @ COMPLEX A(R,R),B(R),X(R)
30 DATA (100,-116.6666667),(0,266.6666667),(0,0),(0,0)
31 DATA (0,266.6666667),(1E6,-116.6666667),(-1E6,0),(0,0)
32 DATA (0,0),(-1E6,0),(1E6,147.3333333),(0,-150)
33 DATA (0,0),(0,0),(0,-150),(1E5,147.3333333)
40 FOR I=1 TO R @ FOR J=1 TO R @ READ A(I,J) @ NEXT J @ NEXT I
45 DATA 10,0,0,0
50 FOR I=1 TO R @ READ B(I) @ NEXT I

Code:

>RUN
>MAT X=SYS(A,B) @ MAT DISP X
 ( 1.99579513038E-4, 4.09639907507E-3)
 (-1.44883361599E-3,-3.56329830823E-2)
 (-1.45408317044E-3,-3.56327608305E-2)
 ( 5.34458117065E-5,-2.25986825155E-6)

X(2) and X(3) are extremely close.
If we solve X(2), X(3)-X(2) instead, we get *full* 12 digits accuracy.

Code:

>FOR I=1 TO R @ A(I,2)=A(I,2)+A(I,3) @ NEXT I
>MAT X=SYS(A,B) @ X(3)=X(3)+X(2) @ MAT DISP X
 ( 1.99579513351E-4, 4.09639907503E-3)
 (-1.44883361584E-3,-3.56329830822E-2)
 (-1.45408317029E-3,-3.56327608304E-2)
 ( 5.34458117064E-5,-2.25986825132E-6)

MathPac2 Owner's Manual Supplement for HP-71 Wrote:MAT X=INV(A)*B directly computes the solution of the matrix equation A X=B instead of
inverting A then multiplying by B. It is however not exactly equivalent to MAT X=SYS(A,B)
in that it doesn't compute the residual correction that MAT..SYS does.

Newly adjusted A is very stable.
Without resdiual correction, this run twice as fast, without giving up much accuracy.

Code:

>MAT X=INV(A)*B @ X(3)=X(3)+X(2) @ MAT DISP X
 ( 1.99579513352E-4, 4.09639907504E-3)
 (-1.44883361584E-3,-3.56329830822E-2)
 (-1.45408317029E-3,-3.56327608304E-2)
 ( 5.34458117064E-5,-2.25986825133E-6)

[Werner]

Well, that does it. Proof that I'm not even an average-minus engineer, as I didn't spot this and I should have.
Thanks, Albert! Yet another lesson to remember.

Cheers, Werner

[J-F Garnier]

Thanks Albert for your HP-71B solution.

My solution is based on considerations on the physical system but is not so different as we will see later.

As Werner noted, the condition number of the problem is very high, and is due to the presence of the R2=1E6 values in two rows of the matrix:



And here starts a kind of paradoxical situation:
- from a physical point of view, when R2 becomes higher and higher, it is more and more negligible in the electrical circuit,
- but from a mathematical point of view, the matrix becomes more and more ill-conditioned because the two row 2 and 3 becomes closer and closer.

This lead me to think that the mathematical description of the problem was not relevant for high values of R2.
If R2 becomes negligible , then the two central loops can be replaced by a single one (in red below).
So to get the better representation for high R2 values, I replaced one of these loops with this one in red, and keep another one of the two.
The system description becomes:



and the problem can be expressed numerically by the complex matrix equation:

| (100,-116.6666667) (0,266.6666667)    (0,0)             (0,0)              |   | I1 |   | 10 |
| (0,266.6666667)    (0,-116.6666667)   (0,147.3333333)   (0,-150)           |   | I2 |   | 0  |
| (0,0)              (-1E6,0)           (1E6,147.3333333) (0,-150)           | x | I3 | = | 0  |
| (0,0)              (0,0)              (0,-150)          (1E5,147.3333333)  |   | I4 |   | 0  |

This system can be easily solved, and the 15c results become:

I1 =  1.995795134 e-4  ,  4.096399080 e-3
I2 = -1.448833618 e-3  , -3.563298307 e-2
I3 = -1.454083172 e-3  , -3.563276083 e-2
I4 =  5.344581171 e-5  , -2.259868254 e-6

now very close to the Free42 reference solution:

I1 =  1.995795134 e-4  ,  4.096399075 e-3
I2 = -1.448833616 e-3  , -3.563298308 e-2
I3 = -1.454083170 e-3  , -3.563276083 e-2
I4 =  5.344581171 e-5  , -2.259868251 e-6

The problem can be solved on a 41C + Advantage Pac too.
Here are the corresponding results:

I1 =  1.995795132 e-4  ,  4.096399077 e-3
I2 = -1.448833617 e-3  , -3.563298309 e-2
I3 = -1.454083171 e-3  , -3.563276084 e-2
I4 =  5.344581172 e-5  , -2.259868252 e-6

that are globally as good as on the 15C, but different.
Indeed, the matrix functions of the Advantage Pac are based on code from the german CCD module, not on HP code from the 15C.

A last remark, some of you surely already noticed that all my considerations on the physical system just lead to add the row 2 and 3 of the initial system matrix, and replace one of the rows with it.
Albert added two columns and I added two rows. The math and physical points of view meet at the end!

J-F
   

[Albert Chan]

Adding 2 columns, solving I3-I2 instead of I3, is slightly superior than adding 2 rows.

Code:

| (100,-116.6666667) (0,266.6666667)  (0,0)             (0,0)             |   | I1  |   | 10 |
| (0,266.6666667)    (0,-116.6666667) (-1E6,0)          (0,0)             |   | I2  |   | 0  |
| (0,0)              (0,147.3333333)  (1E6,147.3333333) (0,-150)          | x |I3-I2| = | 0  |
| (0,0)              (0,-150)         (0,-150)          (1E5,147.3333333) |   | I4  |   | 0  |

We get a more stable matrix.

Even if I3-I2 is not super accurate, restored I3 will still be good.
This translated to more accuracy for the other currents.

Example, this is current residual corrections, difference of SYS(A,B) vs INV(A)*B

R2 = 1e6 ohms (OP example)

Code:

ulp error   I1      I2      I3      I4
Add rows    (-1,5)  (2,0)   (2,-1)  (2,-1)
Add cols    (0,-1)  (0,0)   (0,0)   (0,-1)

R2 = 1e9 ohms

Code:

ulp error   I1      I2      I3      I4
Add rows    (7,13)  (4,1)   (5,3)   (5,7)
Add cols    (-2,-1) (-1,0)  (-1,0)  (0,-1)