How can I solve a system in complex variables ?

8 次查看（过去 30 天）
Michele
Michele 2013-2-7
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
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2013-2-7
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
Michele 2013-2-11
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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Miroslav Balda
Miroslav Balda 2013-2-12
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!
  1 个评论
Michele
Michele 2013-2-12
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?

请先登录，再进行评论。

更多回答（3 个）

Miroslav Balda
Miroslav Balda 2013-2-7
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.
  2 个评论
Miroslav Balda
Miroslav Balda 2013-2-7
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
Michele 2013-2-11
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

请先登录，再进行评论。

Michele
Michele 2013-2-11
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.
Michele
Michele 2013-2-11
编辑：Michele 2013-2-11
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.

类别

在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by