How can I solve a system in complex variables ?

Hi everyone,

I am founding lot of problem in solving a system in the form A*a - B*b=0 where A and B are matrices and a and b vector columns, where inside A and B I have got my complex unknowns; the number of unknowns match the number of equations, but the problem is that the unknown are inside the matrices and I am not able to simplify the problem. To call the complex unknown I have used the sym function, and to solve the problem respect to my functions I have used the comand "solve()". I report my system below: as you can see the only unknowns are G and N1, while s B C g are terms that depends upon these unknowns. Once defined these terms I have defined my matrices and vectors (the vectors are compose by known values).

for x1=(0:0.1:1);

for x2=(1.1:0.1:2);

G = sym('G');

N1 = sym('N1');

beta_s= (mi)/ro_s;

beta_l= (-i*om*G)/ro_l;

s_s= N1/beta_s;

s_l= N1/beta_l;

B_s= sqrt(1-beta_s^2*s_s^2);

B_l= sqrt(1-beta_s^2*s_s^2);

C_s= 1-2*beta_s^2*s_s^2;

C_l= 1-2*beta_l^2*s_l^2;

g_s= exp(i*om*B_s*beta_s^-1*x1);

g_l= exp(i*om*B_s*beta_s^-1*x2);

A1=[0 
i*om*ro_s*s_s*beta_s^-1*B_s*g_s
 0 -i*om*ro_s*s_s*beta_s^-1*B_s*g_s;
 0 
i*om*ro_s*s_s*beta_s*C_s*g_s^-1
 0 
-i*om*ro_s*beta_s*C_s*g_s;
 0 
-B_s*g_s^-1 
0
 B_s*g_s;...

      0 
  -beta_s*s_s*g_s^-1
   0 
  -beta_s*s_s*g_s];

A2=[0 
i*om*ro_l*s_l*beta_l^-1*B_l*g_l
 0
 -i*om*ro_l*s_l*beta_l^-1*B_l*g_l;
 0 
i*om*ro_l*s_l*beta_l*C_l*g_l^-1 
0 
-i*om*ro_l*beta_l*C_l*g_l;
 0 
-B_l*g_l^-1 0 B_l*g_s;...

      0
   -beta_l*s_l*g_l^-1 
  0
   -beta_s*s_l*g_l];

V_1=[0; Rs1;0; Ts1];

V_2=[0;0;0;Ts2];

S= solve(A1*V_1-A2*V_2);

end

end

In the end I come out with this error

??? Error using ==> maple at 129
at offset 480, `;` unexpected

Error in ==> sym.findsym at 33
v = maple('indets', sc ,'symbol');

Error in ==> solve at 99
   vars = ['[' findsym(sym(eqns),neqns) ']'];

Error in ==> sym.solve at 49
[varargout{1:max(1,nargout)}] = solve(S{:});

Error in ==> xxx at 39
S=solve(A1*V_1-A2*V_2);

So my question is: is there a way to solve this kind of problem in Matlab? Can I obtain an explicit expression for my variable G for a certain value of x in the form G(x)= Re+Im?

Thank you very much for the help

[ChristianW]

Can you reproduce the error/problem with minimal and independent code? If yes, post the code please.

The function solve gives complex solutions:

syms x
solve('x^4 + 1 = 2*x^2 - 1')

[Walter Roberson]

Please put a breakpoint in at the solve() call, and run, and when it stops, please show us the values of A1, V_1, A2, V_2, and A1*V_1-A2*V_2

[Michele]

Dear Walter, setting a breakpoint at the solve line I obtain the following values for A1,A2,V1 and V2 A1 =

[ 0, 2821266740684990234375/1436291958794964824997567995651750363136*i*N1*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2), 0, -2821266740684990234375/1436291958794964824997567995651750363136*i*N1*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2)]

[ 0, 20000000000*i*N1*(1-947502880183141812940172165120/473751440091570878575219460881*N1^2), 0, -76923076923076928*i*(1-947502880183141812940172165120/473751440091570878575219460881*N1^2)] [ 0, -1/688296041025641*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2), 0, 1/688296041025641*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2)] [ 0, -N1, 0, -N1]

A2 =

[ 0, -100/688296041025641*i*N1/G^2*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2)*exp(49654294636055828125/11949949097174107789150060111986688*i*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2)), 0, 100/688296041025641*i*N1/G^2*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2)*exp(49654294636055828125/11949949097174107789150060111986688*i*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2))] [ 0, 10000000000*i*N1*(1-2*N1^2)/exp(49654294636055828125/11949949097174107789150060111986688*i*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2)), 0, -100000000000000*G*(1-2*N1^2)*exp(49654294636055828125/11949949097174107789150060111986688*i*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2))] [ 0, -1/688296041025641*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2)/exp(49654294636055828125/11949949097174107789150060111986688*i*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2)), 0, 1/688296041025641*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2)] [ 0, -N1/exp(49654294636055828125/11949949097174107789150060111986688*i*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2)), 0, -2064888123076923/5368709120000*i*N1/G*exp(49654294636055828125/11949949097174107789150060111986688*i*(473751440091570878575219460881-473751440091570906470086082560*N1^2)^(1/2))]

 V1

 0,00000000000000 + 0,00000000000000i
0,800000000000000 + 0,100000000000000i
0,00000000000000 + 0,00000000000000i
1,00000000000000 + 0,00000000000000i

V2= 0,00000000000000 + 0,00000000000000i 0,00000000000000 + 0,00000000000000i 0,00000000000000 + 0,00000000000000i 0,200000000000000 - 0,100000000000000i

Answer 1

[Miroslav Balda]

The structure of the MATLAB code without an application of the Symbolic Toolbox should be as follows:

1. All constant parameters not depending on unknown ones should be defined outside the res function, otherwise they would be assigned in each call of res, what would make computations slower.

2. The function LMFnlsq works with real parameters only. It means that also column vector p0 of the initial guess should be real, composed out of components of complex numbers.

3. The code should be kept clear and all visiblewithout hidden parts.

The following code has the right structure:

%   MAIN.m
%   ~~~~~~  
%  12.02.2013 Michele How can I solve asystem of complex variables
%
global p0 ro_s E ni mi ro_l om Ts1 Rs1 Ts2 x1 x2 beta_s V_1 V_2
%
ro_s = 2000;%   insert density of the solid in kg/m3
E = 20*10^9;%   stiffness module for the solid (= 2e10)
ni = 0.3;
mi = E/(2*(1+ni));  %   for the solid
ro_l = 1000;        %   insert density of the liquid in kg/m3
om = 10*10^6;       %   insert frequency in Hz (=10e6)
Ts1 = 1;            %   amplitude of the wave transmitted in medium 1
Rs1 = 0.8 + i*0.1;
Ts2 = 1- Rs1;
V_1 = [0; Rs1;0; Ts1];
V_2 = [0;0;0;Ts2];
x1  = 5*10^-3;
x2  = x1; 
beta_s = (mi)/ro_s;
%
%p0 = [10^6;10^6; 1; 1]; %   first guess (it gives NaNs)
p0 = [10^1;10^1; 1; 1];  %   first guess
%
[p,ssq,cnt] = LMFnlsq(@res, ones(size(p0)), 'Display',-5);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%           The solution:
p = p.*p0;
G  = complex(p(1),p(3));
N1 = complex(p(2),p(4));
%
fprintf('\n  Sum of squares    = %10.2e\n  Function evaluations = %g\n',...
  ssq,cnt);
fprintf('  G  = %12.4e + %12.4ei\n  N1 = %12.4e + %12.4ei\n',...
  p(1),p(3), p(2),p(4));

It is clear that the constant terms are evaluated outside the iteration cycle repesented by the function LMFnlsq. The supporting function res.m evaluating residuals of the system of equations has the form

 function r = res(p)
%~~~~~~~~~~~~~~~~~
global p0 ro_s E ni mi ro_l om Ts1 Rs1 Ts2 x1 x2 beta_s V_1 V_2
%
p = p.*p0;
G  = complex(p(1),p(3));
N1 = complex(p(2),p(4));
beta_l = (-i*om*G)/ro_l;
s_s = N1/beta_s;
s_l = N1/beta_l;
B_s = sqrt(1-beta_s^2*s_s^2);
B_l = sqrt(1-beta_s^2*s_s^2);
C_s = 1-2*beta_s^2*s_s^2;
C_l = 1-2*beta_l^2*s_l^2;
g_s = exp(i*om*B_s*beta_s^-1*x1);
g_l = exp(i*om*B_s*beta_s^-1*x2);
%
A1 = [0, i*om*ro_s*s_s*beta_s^-1*B_s*g_s, 0, -i*om*ro_s*s_s*beta_s^-1*B_s*g_s;
 0, i*om*ro_s*s_s*beta_s*C_s*g_s^-1, 0, -i*om*ro_s*beta_s*C_s*g_s;
 0, -B_s*g_s^-1, 0, B_s*g_s;
 0, -beta_s*s_s*g_s^-1, 0, -beta_s*s_s*g_s]; 
%
A2 = [0 i*om*ro_l*s_l*beta_l^-1*B_l*g_l 0 -i*om*ro_l*s_l*beta_l^-1*B_l*g_l;
 0 i*om*ro_l*s_l*beta_l*C_l*g_l^-1 0 -i*om*ro_l*beta_l*C_l*g_l;
 0 -B_l*g_l^-1 0 B_l*g_s;
 0 -beta_l*s_l*g_l^-1 0 -beta_s*s_l*g_l];
%
rx = A1*V_1 - A2*V_2;     %   column vector of complex residuals
r = [real(rx(1));imag(rx(1)); real(rx(2));imag(rx(2))];

The statements remained unchanged but some formal changes leading to better visibility. A1 and A2 could be written better in the complex form. Thus formulated task gave the following results:

>> main
**************************************************
itr  nfJ   SUM(r^2)        x           dx
**************************************************
 0    1  3.1035e+038  1.0000e+000  0.0000e+000 
                      1.0000e+000  0.0000e+000 
                      1.0000e+000  0.0000e+000 
                      1.0000e+000  0.0000e+000 
 5    6  4.7887e+024 -1.5385e+002  7.5141e-003 
                      6.4562e-001  1.7133e-009 
                      3.0769e+003 -1.5147e-001 
                      6.5933e-001  1.7496e-009 

   8    9  3.1034e+011 -1.5385e+002  1.9139e-009 
                        6.4562e-001 -4.0313e-018 
                        3.0769e+003 -3.8455e-008 
                        6.5933e-001  6.1153e-018 

Sum of squares    =  3.10e+011
Function evaluations = 8
G  = -1.5385e+003 +  3.0769e+003i
N1 =  6.4562e+000 +  6.5933e-001i

It is clear that the recommended process operates, nevertheless, the numerical values of results are dependent on the original data and formulas. It has been observed that the original initial guess is very far from stable solution (the resultshave been all NaNs. After changing it, the process converged. While G has been almost insensitive on the guess, the N1 has changed strongly with it.

I regard the presented process as solved.

Good luck!

[Michele]

Dear Miroslav, that's great, thank you. I agree with you on the fact that the algorithm works. The only thing that I don't understand is how can you say that the algorithm is convergent?

Answer 2

[Miroslav Balda]

Maybe that you can get a solution in a simpler way just in MATLAB without using symbolic approach. Your equation is complex dependent on a complex vector p composed out of 2 complex unknowns G and N1. The solution can be obtained by the following (pseudo)code:

p0 = [real(G0);real(N10); imag(G0); imag(N10)]; % initial estimates of unknowns
[p,ssq,cnt] = LMFnlsq(@res, ones(size(p0)), 'Display',-10);
%  LMFnlsq in FEX at  www.mathworks.com/matlabcentral/fileexchange/17534
p  = p.*p0;
%  The solution:
G  = complex(p(1),p(3))
N1 = complex(p(2),p(4))
fprintf('Sum of squares      = %10.2e\n...
   'Function evaluations  %g\n...
   'G  = %10.2e + %10.2ei\n...
   'NI = %10.2e + %10.2ei\n', ssq, cnt, p(1),p(3), p(2),p(4));

It is neccessary to build the function for evaluating differences of the given matrix equation from zero:

function r = res(p)
p  = p.*p0;
G  = complex(p(1),p(3))
N1 = complex(p(2),p(4))
 %  Evaluate a complex residual d between the equation and zero as
 %  d = A(p)*a(p)-B(p)*b(p); % and set the vector of its components: 
r  = [real(d); imag(d)]; % residuals

As is obvious, it is necessary to solve a system of 2 equations for 2 complex unknowns.

[Miroslav Balda]

More rigorous would be to write d = A(q)*a(q)-B(q)*b(q), where

q = [G;N1]; instead of d = A(p)*a(p)-B(p)*b(p) in the function res.

[Michele]

Dear Miroslav, following your advice I have written two codes (one posted as a general reply and this one as a reply to your comment): function r = res(p)

p0 = [10^6;10^6; 1; i]; %first guess

[p,ssq,cnt] = LMFnlsq(@res, ones(size(p0)), 'Display',-10);

p = p.*p0;

% The solution:

G = complex(p(1),p(3));

N1 = complex(p(2),p(4));

fprintf('Sum of squares = %10.2e\n','Function evaluations = %g\n ','G = %10.2e + %10.2ei\n ','NI = %10.2e + %10.2ei\n', ssq, cnt, p(1),p(3), p(2),p(4));

ro_s=2000; %insert density of the solid in kg/m3

E= 20*10^9; %stiffness module for the solid

ni=0.3;

mi= E/(2*(1+ni)); %for the solid

ro_l=1000; %insert density of the liquid in kg/m3

om= 10*10^6; %%insert frequency in Hz

Ts1= 1; %amplitude of the wave transmitted in medium 1

Rs1= 0.8 + i*0.1;

Ts2= 1- Rs1;

x1=5*10^-3;

x2=x1; beta_s= (mi)/ro_s;

beta_l= (-i*om*G)/ro_l;

s_s= N1/beta_s;

s_l= N1/beta_l;

B_s= sqrt(1-beta_s^2*s_s^2);

B_l= sqrt(1-beta_s^2*s_s^2);

C_s= 1-2*beta_s^2*s_s^2;

C_l= 1-2*beta_l^2*s_l^2;

g_s= exp(i*om*B_s*beta_s^-1*x1);

g_l= exp(i*om*B_s*beta_s^-1*x2);

A1=[0 i*om*ro_s*s_s*beta_s^-1*B_s*g_s 0 -i*om*ro_s*s_s*beta_s^-1*B_s*g_s; 0 i*om*ro_s*s_s*beta_s*C_s*g_s^-1 0 -i*om*ro_s*beta_s*C_s*g_s; 0 -B_s*g_s^-1 0 B_s*g_s;... 0 -beta_s*s_s*g_s^-1 0 -beta_s*s_s*g_s]; A2=[0 i*om*ro_l*s_l*beta_l^-1*B_l*g_l 0 -i*om*ro_l*s_l*beta_l^-1*B_l*g_l; 0 i*om*ro_l*s_l*beta_l*C_l*g_l^-1 0 -i*om*ro_l*beta_l*C_l*g_l; 0 -B_l*g_l^-1 0 B_l*g_s;... 0 -beta_l*s_l*g_l^-1 0 -beta_s*s_l*g_l]; V_1=[0; Rs1;0; Ts1]; V_2=[0;0;0;Ts2];

 %  Evaluate a complex residual d between the equation and zero as

d = A1(p)*V_1(p)-A2(p)*V_2(p); % and set the vector of its components:

r = [real(d); imag(d)]; % residuals

The result is in both case that the problem is not converging to a solution. In this case I obtain this error ??? Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N) to change the limit. Be aware that exceeding your available stack space can crash MATLAB and/or your computer.

Error in ==> LMFnlsq

I think I may reconsider the equations that I want to implement and before I accept your answer as the best one may I ask if you can check the code I have written to see if I modified your advice in the way you meant?

Still thank you for your help

Answer 3

[Michele]

Miroslav, sorry if I reply only now. Thank you for the answer, I will go through your advice and I will see if it does work.

Answer 4

[Michele]

Dear Miroslav, Walter and Christian, do you think that a code like that may be successfull? :

function solveeqs()

guess=[a b]% I give 2 values for my unknown to start an interpolation

options = optimset('MaxFunEvals',1000000,'MaxIter',100000);

[result,fval,exit,output]=fsolve(@eqns,guess,options);

result

fval

eqns(guess)

output % these are the output of my solve comand

end

function fcns=eqns(z)

X=z(1); % I say which are the unknowns

Y=z(2);

fcns(1)=function of X and Y %I give the functions (note that if the number of unknowns is < of the number of equations matlab will run anyway the program choosing an algorithm on his own)

fcns(2)= function of X and Y;

. . . fcns(n)... end

what I am doing is writing in an explicit way the functions instead of using matrices and with this algorithm I give a first guess for the interpolation towards the solution.