Solving a set of equations numerically where solutions are complex numbers

close all
clear all
syms S11 S22 S33 S12 S13  S23;

A11= 0.2346 + 0.1345i;
A21= 0.8550 - 0.4363i;
A22= -0.0412 + 0.2931i;

B11= 0.2346 + 0.1345i;
B21= 0.8550 - 0.4363i;
B22= -0.0412 + 0.2931i;

D11= 0.2346 + 0.1345i;
D21= 0.8550 - 0.4363i;
D22= -0.0412 + 0.2931i;

c11= 0.169033453278504 + 0.127942837737115i;
c22= 0.169951074114244 + 0.148423165394202i;
c33= 0.168939745364374 + 0.142538799523638i;
c21= -0.00172796568379862 + 0.201091742020021i;
c31= -0.00152297772077853 + 0.199585730377658i;
c32= -0.00523160858486265 + 0.199472640309087i;

eqn1 = (A11.*A22.*S11 - A21.^2.*S11 - A11 + A11.*B22.*S22 + A11.*D22.*S33 - A21.^2.*B22.*S12.^2 - A21.^2.*D22.*S13.^2 + A11.*A22.*B22.*S12.^2 + A11.*A22.*D22.*S13.^2 + A11.*B22.*D22.*S23.^2 + A21.^2.*B22.*S11.*S22 + A21.^2.*D22.*S11.*S33 + A21.^2.*B22.*D22.*S11.*S23.^2 + A21.^2.*B22.*D22.*S13.^2.*S22 + A21.^2.*B22.*D22.*S12.^2.*S33 - A11.*A22.*B22.*S11.*S22 - A11.*A22.*D22.*S11.*S33 - A11.*B22.*D22.*S22.*S33 - A11.*A22.*B22.*D22.*S11.*S23.^2 - A11.*A22.*B22.*D22.*S13.^2.*S22 - A11.*A22.*B22.*D22.*S12.^2.*S33 - 2.*A21.^2.*B22.*D22.*S12.*S13.*S23 - A21.^2.*B22.*D22.*S11.*S22.*S33 + 2.*A11.*A22.*B22.*D22.*S12.*S13.*S23 + A11.*A22.*B22.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c11==0
eqn2 = -(A21.*B21.*(S12 + D22.*S13.*S23 - D22.*S12.*S33))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c21==0
eqn3 = -(A21.*D21.*(S13 + B22.*S12.*S23 - B22.*S13.*S22))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c31==0
eqn4 = (A22.*B11.*S11 - B21.^2.*S22 - B11 + B11.*B22.*S22 + B11.*D22.*S33 - A22.*B21.^2.*S12.^2 - B21.^2.*D22.*S23.^2 + A22.*B11.*B22.*S12.^2 + A22.*B11.*D22.*S13.^2 + B11.*B22.*D22.*S23.^2 + A22.*B21.^2.*S11.*S22 + B21.^2.*D22.*S22.*S33 + A22.*B21.^2.*D22.*S11.*S23.^2 + A22.*B21.^2.*D22.*S13.^2.*S22 + A22.*B21.^2.*D22.*S12.^2.*S33 - A22.*B11.*B22.*S11.*S22 - A22.*B11.*D22.*S11.*S33 - B11.*B22.*D22.*S22.*S33 - A22.*B11.*B22.*D22.*S11.*S23.^2 - A22.*B11.*B22.*D22.*S13.^2.*S22 - A22.*B11.*B22.*D22.*S12.^2.*S33 - 2.*A22.*B21.^2.*D22.*S12.*S13.*S23 - A22.*B21.^2.*D22.*S11.*S22.*S33 + 2.*A22.*B11.*B22.*D22.*S12.*S13.*S23 + A22.*B11.*B22.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c22==0
eqn5 = -(B21.*D21.*(S23 + A22.*S12.*S13 - A22.*S11.*S23))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c32==0
eqn6 = (A22.*D11.*S11 - D21.^2.*S33 - D11 + B22.*D11.*S22 + D11.*D22.*S33 - A22.*D21.^2.*S13.^2 - B22.*D21.^2.*S23.^2 + A22.*B22.*D11.*S12.^2 + A22.*D11.*D22.*S13.^2 + B22.*D11.*D22.*S23.^2 + A22.*D21.^2.*S11.*S33 + B22.*D21.^2.*S22.*S33 + A22.*B22.*D21.^2.*S11.*S23.^2 + A22.*B22.*D21.^2.*S13.^2.*S22 + A22.*B22.*D21.^2.*S12.^2.*S33 - A22.*B22.*D11.*S11.*S22 - A22.*D11.*D22.*S11.*S33 - B22.*D11.*D22.*S22.*S33 - A22.*B22.*D11.*D22.*S11.*S23.^2 - A22.*B22.*D11.*D22.*S13.^2.*S22 - A22.*B22.*D11.*D22.*S12.^2.*S33 - 2.*A22.*B22.*D21.^2.*S12.*S13.*S23 - A22.*B22.*D21.^2.*S11.*S22.*S33 + 2.*A22.*B22.*D11.*D22.*S12.*S13.*S23 + A22.*B22.*D11.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c33==0

[S11, S22, S33, S12, S13,  S23]=vpasolve([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6], [S11, S22, S33, S12, S13,  S23])

What am I doing wrong? I'm not getting any answers. Is there any other ways such as using fsolve and FminCon?

Answer 1

[John D'Errico]

Not all nasty looking, high order nonlinear systems of equations have solutions. vpasolve is a capable numerical solver, but it has limits. One issue is that when working in high precision, things take longer to do. All computations are slower. It will probably happen that vpasolve will return a high precision solution eventually.

So I tried posing this as a problem for fsolve. Working in double precision is far faster. fsolve can sometimes have issues with complex valued functions. But it often does seem to work.

Here is your code, pasted into a function to be used for fsolve.

function eqs = funs(S)
S11 = S(1);
S22 = S(2);
S33 = S(3);
S12 = S(4);
S13 = S(5);
S23 = S(6);
A11= 0.2346 + 0.1345i;
A21= 0.8550 - 0.4363i;
A22= -0.0412 + 0.2931i;
B11= 0.2346 + 0.1345i;
B21= 0.8550 - 0.4363i;
B22= -0.0412 + 0.2931i;
D11= 0.2346 + 0.1345i;
D21= 0.8550 - 0.4363i;
D22= -0.0412 + 0.2931i;
c11= 0.169033453278504 + 0.127942837737115i;
c22= 0.169951074114244 + 0.148423165394202i;
c33= 0.168939745364374 + 0.142538799523638i;
c21= -0.00172796568379862 + 0.201091742020021i;
c31= -0.00152297772077853 + 0.199585730377658i;
c32= -0.00523160858486265 + 0.199472640309087i;
eqs = [(A11.*A22.*S11 - A21.^2.*S11 - A11 + A11.*B22.*S22 + A11.*D22.*S33 - A21.^2.*B22.*S12.^2 - A21.^2.*D22.*S13.^2 + A11.*A22.*B22.*S12.^2 + A11.*A22.*D22.*S13.^2 + A11.*B22.*D22.*S23.^2 + A21.^2.*B22.*S11.*S22 + A21.^2.*D22.*S11.*S33 + A21.^2.*B22.*D22.*S11.*S23.^2 + A21.^2.*B22.*D22.*S13.^2.*S22 + A21.^2.*B22.*D22.*S12.^2.*S33 - A11.*A22.*B22.*S11.*S22 - A11.*A22.*D22.*S11.*S33 - A11.*B22.*D22.*S22.*S33 - A11.*A22.*B22.*D22.*S11.*S23.^2 - A11.*A22.*B22.*D22.*S13.^2.*S22 - A11.*A22.*B22.*D22.*S12.^2.*S33 - 2.*A21.^2.*B22.*D22.*S12.*S13.*S23 - A21.^2.*B22.*D22.*S11.*S22.*S33 + 2.*A11.*A22.*B22.*D22.*S12.*S13.*S23 + A11.*A22.*B22.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c11;
  -(A21.*B21.*(S12 + D22.*S13.*S23 - D22.*S12.*S33))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c21;
  -(A21.*D21.*(S13 + B22.*S12.*S23 - B22.*S13.*S22))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c31;
   (A22.*B11.*S11 - B21.^2.*S22 - B11 + B11.*B22.*S22 + B11.*D22.*S33 - A22.*B21.^2.*S12.^2 - B21.^2.*D22.*S23.^2 + A22.*B11.*B22.*S12.^2 + A22.*B11.*D22.*S13.^2 + B11.*B22.*D22.*S23.^2 + A22.*B21.^2.*S11.*S22 + B21.^2.*D22.*S22.*S33 + A22.*B21.^2.*D22.*S11.*S23.^2 + A22.*B21.^2.*D22.*S13.^2.*S22 + A22.*B21.^2.*D22.*S12.^2.*S33 - A22.*B11.*B22.*S11.*S22 - A22.*B11.*D22.*S11.*S33 - B11.*B22.*D22.*S22.*S33 - A22.*B11.*B22.*D22.*S11.*S23.^2 - A22.*B11.*B22.*D22.*S13.^2.*S22 - A22.*B11.*B22.*D22.*S12.^2.*S33 - 2.*A22.*B21.^2.*D22.*S12.*S13.*S23 - A22.*B21.^2.*D22.*S11.*S22.*S33 + 2.*A22.*B11.*B22.*D22.*S12.*S13.*S23 + A22.*B11.*B22.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c22;
  -(B21.*D21.*(S23 + A22.*S12.*S13 - A22.*S11.*S23))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c32;
  (A22.*D11.*S11 - D21.^2.*S33 - D11 + B22.*D11.*S22 + D11.*D22.*S33 - A22.*D21.^2.*S13.^2 - B22.*D21.^2.*S23.^2 + A22.*B22.*D11.*S12.^2 + A22.*D11.*D22.*S13.^2 + B22.*D11.*D22.*S23.^2 + A22.*D21.^2.*S11.*S33 + B22.*D21.^2.*S22.*S33 + A22.*B22.*D21.^2.*S11.*S23.^2 + A22.*B22.*D21.^2.*S13.^2.*S22 + A22.*B22.*D21.^2.*S12.^2.*S33 - A22.*B22.*D11.*S11.*S22 - A22.*D11.*D22.*S11.*S33 - B22.*D11.*D22.*S22.*S33 - A22.*B22.*D11.*D22.*S11.*S23.^2 - A22.*B22.*D11.*D22.*S13.^2.*S22 - A22.*B22.*D11.*D22.*S12.^2.*S33 - 2.*A22.*B22.*D21.^2.*S12.*S13.*S23 - A22.*B22.*D21.^2.*S11.*S22.*S33 + 2.*A22.*B22.*D11.*D22.*S12.*S13.*S23 + A22.*B22.*D11.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c33];

After saving this function on my search path, I used fsolve, as follows:

format long g
Sfinal = fsolve(@funs,(1+i)*ones(1,6))

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

Sfinal =
Columns 1 through 3
      -0.0576463452402932 -    0.0753664180207837i        -0.0744682052419572 -    0.0630375797716533i        -0.0698118920703255 -    0.0672129482711116i
Columns 4 through 6
       -0.185173231731065 +     0.109260636702373i         -0.183468124937688 +      0.10862578397362i         -0.186854054845906 +     0.104967752613526i

As a test that this is a solution in double precision, evaluate funs at that point.

funs(Sfinal)
ans =
       2.04330996567137e-12 +  8.08547673258886e-13i
       2.06829562157673e-12 +  8.28642710004601e-13i
       2.07329856408145e-12 +  8.41465785939022e-13i
       2.07303618715571e-12 +   7.6394446324457e-13i
       2.08355208086708e-12 +  8.14071032806396e-13i
       2.05693795329864e-12 +  7.88147325181399e-13i

Note that when calling fsolve, it was necessary to use a complex start point. Even then, I have seen cases where fsolve fails for complex valued problems, though I have not yet pinned down why it may fail. It did seem to work here though.

[John D'Errico]

I'll return to this later on today, to show why it is that VPASOLVE fails to finish in a reasonable time.

[John D'Errico]

Back now...

So why is vpasolve having a problem here? The answer is simple, if you appreciate how the symbolic toolbox works with numbers, how they are stored.

Create a simple number, a.

a = sym(1.23);
a
a =
123/100

So MATLAB thinks of that number as the ratio of two integers. It works splendidly well sometimes.

a^2 
ans =
15129/10000
>> a^27
ans =
267570363636059316998473870115684948492796490598096370147/1000000000000000000000000000000000000000000000000000000

So if we raise a to some high power, or do any computations with it, that number is now an integer fraction, but with LOTS of digits.

See what happens if we try to force this number to be stored as a higher precision float.

b = vpa(a)
b =
1.23

b^27
ans =
267.57036363605931699847387011568

b^127
ans =
261787648930.5263626282005527684

So now MATLAB stores any computations as floating point numbers, with around 31 decimal digits.

But, why does any of that matter to you?

eqn1 = (A11.*A22.*S11 - A21.^2.*S11 - A11 + A11.*B22.*S22 + A11.*D22.*S33 - A21.^2.*B22.*S12.^2 - A21.^2.*D22.*S13.^2 + A11.*A22.*B22.*S12.^2 + A11.*A22.*D22.*S13.^2 + A11.*B22.*D22.*S23.^2 + A21.^2.*B22.*S11.*S22 + A21.^2.*D22.*S11.*S33 + A21.^2.*B22.*D22.*S11.*S23.^2 + A21.^2.*B22.*D22.*S13.^2.*S22 + A21.^2.*B22.*D22.*S12.^2.*S33 - A11.*A22.*B22.*S11.*S22 - A11.*A22.*D22.*S11.*S33 - A11.*B22.*D22.*S22.*S33 - A11.*A22.*B22.*D22.*S11.*S23.^2 - A11.*A22.*B22.*D22.*S13.^2.*S22 - A11.*A22.*B22.*D22.*S12.^2.*S33 - 2.*A21.^2.*B22.*D22.*S12.*S13.*S23 - A21.^2.*B22.*D22.*S11.*S22.*S33 + 2.*A11.*A22.*B22.*D22.*S12.*S13.*S23 + A11.*A22.*B22.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c11;

eqn2 = -(A21.*B21.*(S12 + D22.*S13.*S23 - D22.*S12.*S33))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c21;

eqn3 = -(A21.*D21.*(S13 + B22.*S12.*S23 - B22.*S13.*S22))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c31;

eqn4 = (A22.*B11.*S11 - B21.^2.*S22 - B11 + B11.*B22.*S22 + B11.*D22.*S33 - A22.*B21.^2.*S12.^2 - B21.^2.*D22.*S23.^2 + A22.*B11.*B22.*S12.^2 + A22.*B11.*D22.*S13.^2 + B11.*B22.*D22.*S23.^2 + A22.*B21.^2.*S11.*S22 + B21.^2.*D22.*S22.*S33 + A22.*B21.^2.*D22.*S11.*S23.^2 + A22.*B21.^2.*D22.*S13.^2.*S22 + A22.*B21.^2.*D22.*S12.^2.*S33 - A22.*B11.*B22.*S11.*S22 - A22.*B11.*D22.*S11.*S33 - B11.*B22.*D22.*S22.*S33 - A22.*B11.*B22.*D22.*S11.*S23.^2 - A22.*B11.*B22.*D22.*S13.^2.*S22 - A22.*B11.*B22.*D22.*S12.^2.*S33 - 2.*A22.*B21.^2.*D22.*S12.*S13.*S23 - A22.*B21.^2.*D22.*S11.*S22.*S33 + 2.*A22.*B11.*B22.*D22.*S12.*S13.*S23 + A22.*B11.*B22.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c22;

eqn5 = -(B21.*D21.*(S23 + A22.*S12.*S13 - A22.*S11.*S23))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c32;

eqn6 = (A22.*D11.*S11 - D21.^2.*S33 - D11 + B22.*D11.*S22 + D11.*D22.*S33 - A22.*D21.^2.*S13.^2 - B22.*D21.^2.*S23.^2 + A22.*B22.*D11.*S12.^2 + A22.*D11.*D22.*S13.^2 + B22.*D11.*D22.*S23.^2 + A22.*D21.^2.*S11.*S33 + B22.*D21.^2.*S22.*S33 + A22.*B22.*D21.^2.*S11.*S23.^2 + A22.*B22.*D21.^2.*S13.^2.*S22 + A22.*B22.*D21.^2.*S12.^2.*S33 - A22.*B22.*D11.*S11.*S22 - A22.*D11.*D22.*S11.*S33 - B22.*D11.*D22.*S22.*S33 - A22.*B22.*D11.*D22.*S11.*S23.^2 - A22.*B22.*D11.*D22.*S13.^2.*S22 - A22.*B22.*D11.*D22.*S12.^2.*S33 - 2.*A22.*B22.*D21.^2.*S12.*S13.*S23 - A22.*B22.*D21.^2.*S11.*S22.*S33 + 2.*A22.*B22.*D11.*D22.*S12.*S13.*S23 + A22.*B22.*D11.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c33;

F = [eqn1;eqn2;eqn3;eqn4;eqn5;eqn6];

Now, lets write a simple Newton iteration scheme to solve this problem. What could possibly go wrong?

function Sfinal = symnewton(F,S0,vars,tol,maxiter)
J = jacobian(F,vars);
iter = 0;
err = inf;
while (iter < maxiter) && (err > tol)
  t1 = tic;
  iter = iter + 1;
  Sfinal = S0 - (subs(J,vars,S0)\subs(F,vars,S0)).';
  err = norm(double(S0 - Sfinal));
  t2 = toc(t1);

    disp([iter,t2])
    S0 = Sfinal;
  end

A really basic Newton's method there. But it will show you what happens, and why there are problems.

SNewt = symnewton(F,(1+i)*ones(1,6),[S11 S22 S33 S12 S13 S23],0,1)
            1       0.17747
SNewt =
[ - 8122605055456855530629197935866662124570267482027917800660450120385567372272060370235562840464320379719473298307414287516387150493739654321254487441644253381450181992636494636873879/4233285397525913363432932526997088221795002840574897231408923049142896645070425032678693775050665658304738725206282756786353420993747249141660799141296215152703992811933187178496000 + 3125385859109831821285351000802937656225537794602871959324823423750730975617595905496233817418541817175340830480986304005085591800079519071774117260338578296158885639259010079466759i/12699856192577740090298797580991264665385008521724691694226769147428689935211275098036081325151996974914216175618848270359060262981241747424982397423888645458111978435799561535488000, - 41147918070278363667081096829199491105426161521685145166591673284268100281403142623648940542362185053102939305632457854190912555939566203421414000057747276290067609882368612392327923/21166426987629566817164662634985441108975014202874486157044615245714483225352125163393468875253328291523693626031413783931767104968736245708303995706481075763519964059665935892480000 + 17582419410550866091125161397713261013874228847691116060875289542234399462727195710994812055446811083394886769003091729529023348227372090604032673452278714479647605484415405102053667i/63499280962888700451493987904956323326925042608623458471133845737143449676056375490180406625759984874571080878094241351795301314906208737124911987119443227290559892178997807677440000, - 207965079411716607294654716254784839390940358425651104008008762274876949111938974548934275714962954662114437129711514207244819232946711161409879/107296938885185838774090890429792768140686871142540532469896302864330331999177455056288372780497119183210174621464835439574407174152582594560000 + 28500343109474456957625055525683511759421905196308493553946738394982507302008320671053601528933261763631184427432622189910243108926414585376237i/107296938885185838774090890429792768140686871142540532469896302864330331999177455056288372780497119183210174621464835439574407174152582594560000, - 664545329241691026416398906853416867373027730681613761900362996942596241812318973064640971267593602221342419467343805896933312652743732170404869/321890816655557516322272671289378304422060613427621597409688908592990995997532365168865118341491357549630523864394506318723221522457747783680000 + 47523633026026185942466725919845745435451354644913422029138572727670099415309123880048509488473082087436646713119154179402930847494781396889069i/107296938885185838774090890429792768140686871142540532469896302864330331999177455056288372780497119183210174621464835439574407174152582594560000, - 43653975711624587315997628666562573402040774828283297342769337666514111184975993808901119086755926595022547421875641067363465321773473647361447176115329931247553717730545150696708851/21166426987629566817164662634985441108975014202874486157044615245714483225352125163393468875253328291523693626031413783931767104968736245708303995706481075763519964059665935892480000 + 27864813709732532008917126178864901128673713801042737947771870862369870034650062643302831010992244669954162399980549679685874379672135455976961186088206776145246649730627269307250979i/63499280962888700451493987904956323326925042608623458471133845737143449676056375490180406625759984874571080878094241351795301314906208737124911987119443227290559892178997807677440000, - 8757002759396037860056503865573775739429101024094919041800679081712366649407046194210042852301988097310670365087105914586961458100417963780157892872071513702982776749300381057513623/4233285397525913363432932526997088221795002840574897231408923049142896645070425032678693775050665658304738725206282756786353420993747249141660799141296215152703992811933187178496000 + 5638260673780584383003430579462646581779906168819362270328209902586118224207852663193948276256613332578124602839633076720892719758547688905183375023043343496078313803519538337490183i/12699856192577740090298797580991264665385008521724691694226769147428689935211275098036081325151996974914216175618848270359060262981241747424982397423888645458111978435799561535488000]

Although I forced it to stop after only one iteration, the numbers are immense!

vpa(SNewt,5)
ans =
[ - 1.9187 + 0.2461i, - 1.944 + 0.27689i, - 1.9382 + 0.26562i, - 2.0645 + 0.44292i, - 2.0624 + 0.43882i, - 2.0686 + 0.44396i]

Not immense in terms of the absolute magnitude, but immense in terms of the number of digits carried. Everything is stored in the form of integer fractions, but with now thousands of digits in the numerators and denominators.

If I let the iteration scheme proceed for 3 iterations, the second iteration takes roughly one second to perform. The third iteration took 31 seconds.

Those numbers are growing exponentially in the number of digits carried. So any computations take a HUGE amount of time to perform. This is why vpasolve fails to work happily.

Can we fix it? I'm probably running short of room in this comment, so I'll add another.

[John D'Errico]

This time, I'll just force the symbolic toolbox to recast the result after every iteration using vpa. I can probably do it less crudely.

function Sfinal = symnewton(F,S0,vars,tol,maxiter)
J = jacobian(F,vars);
iter = 0;
err = inf;
while (iter < maxiter) && (err > tol)
  t1 = tic;
  iter = iter + 1;
  Sfinal = S0 - (subs(J,vars,S0)\subs(F,vars,S0)).';
  Sfinal = vpa(Sfinal);
  err = norm(double(S0 - Sfinal));
  t2 = toc(t1);

    disp([iter,err,t2])
    S0 = Sfinal;
  end

Starting from a completely random start point:

SNewt = symnewton(F,rand(1,6) + i*rand(1,6),[S11 S22 S33 S12 S13 S23],1e-30,10)
            1       3.3948      0.16698
            2       2.4699       0.0915
            3       1.2335     0.088967
            4       0.4185     0.089226
            5     0.050441      0.08973
            6   0.00089176     0.093096
            7   2.7782e-07     0.094177
            8   2.6982e-14     0.090472
            9   2.5451e-28     0.095451
           10   2.2873e-40     0.087999
SNewt =
[ - 0.057646345240273125258901514562458 - 0.075366418023328922585198984460681i, - 0.074468205241965717109279556225684 - 0.06303757977421326090733019248076i, - 0.069811892070311265721740690598179 - 0.067212948273664561538665726612952i, - 0.18517323173103393455871448894831 + 0.10926063669978503724287434139558i, - 0.18346812493765326776211584597772 + 0.10862578397102584376285770050721i, - 0.18685405484588656906799822848967 + 0.10496775261093143853085791970942i]

You can see the error decreases with each iteration, and the time per iteration stays constant.

The point is, the problem with the vpasolve solution is the numbers grow without bound. It is also why I chose a different approach for high precision computations in my HPF toolbox.

[supun]

Thank you for suggesting a solution with fsolve. Also, I really like the wonderful discussion, why vpasolve is bad for my application. Thank you John.

[John D'Errico]

vpasolve SHOULD be good for your problem. There SHOULD be a way to force vpasolve to stay with a fixed/finite (even if large) precision, but I don't see one at this point.