solve 4 equations with some unknown coefficients.

I'm attempting to solve this problem. But I don't get the desired answer . Any suggestions?

I want to eliminate these coefficients A,B,C and D from the 4 equations and get one dispersion equation.

syms A B C D beta beta_y d1 alpha_y d2 sigma_ac1 sigma_ac2 omega mu;

eqn1 = A*sin(beta_y*d1) - j*C*exp(-j*beta_y*d1) + j*B*exp(j*beta_y*d1) == 0;
eqn2 = D*alpha_y*exp(-alpha_y*d2) - j*beta_y*(C*exp(-j*beta_y*d2) - B*exp(j*beta_y*d2)) == 0;
eqn3 = A*((j*beta^2)/(omega*mu)*cos(beta_y*d1) - sigma_ac1*beta_y*sin(beta_y*d1)) - (j*beta^2)/(omega*mu)*B*exp(j*beta_y*d1) - (j*beta^2)/(omega*mu)*C*exp(-j*beta_y*d1) == 0;
eqn4 = D*((j*beta^2)/(omega*mu)*exp(-alpha_y*d2) - sigma_ac2*alpha_y*exp(-alpha_y*d2)) - (j*beta^2)/(omega*mu)*B*exp(j*beta_y*d2) - (j*beta^2)/(omega*mu)*C*exp(-j*beta_y*d2) == 0;

eqns = [eqn1, eqn2, eqn3, eqn4];

vars = [A, B, C, D];

sol = solve(eqns, vars);

ASol = sol.A

ASol = 

0

BSol = sol.B

BSol = 

0

CSol = sol.C

CSol = 

0

DSol = sol.D              

DSol = 

0

[Torsten]

The equations represent a linear system in A, B, C and D.

MATLAB gives you the most obvious solution:

A = B = C = D = 0

What do you mean by "The solution is:" ? For what is the equation you write a solution ?

[Ahmad Awada]

I am sorry because i did not explain very well.

My goal is to get a dispersion equation for multi-layer structure. I got 4 equations from the electric and magnetic field equations and after that, I enforced the boundry conditions and finally, I received to these 4 equations with constant coeficients. I want to eliminate these coeficients to get one equation.

I mean, the author of the refrence paper said :We can further manipulate the 4 eqautions, to eliminate all constants A, B, C, and D. This process gives us and he wrot the last equation

How can i get this equation? any suggestions?

I hope my explanation is clear.

Thank you

[Torsten]

Maybe you should enclose the reference paper.

I cannot follow what you mean by

"I want to eliminate these coeficients to get one equation."

[Torsten]

As the page begins with "To find A, B, C, and D, we enforce ...", these constants must have been mentionned before in the article. Could you enclose this part, too ?

[Torsten]

My guess is that the equation given is a condition such that the system of equations has a nontrivial solution for A,B,C and D.

If the determinant of the linear system is 0, this will be the case. But the condition MATLAB gives is complicated; maybe you are able to simplify.

syms beta beta_y d1 alpha_y d2 sigma_ac1 sigma_ac2 omega mu real
syms A B C D
eqn1 = A*sin(beta_y*d1) - j*C*exp(-j*beta_y*d1) + j*B*exp(j*beta_y*d1) == 0;
eqn2 = D*alpha_y*exp(-alpha_y*d2) - j*beta_y*(C*exp(-j*beta_y*d2) - B*exp(j*beta_y*d2)) == 0;
eqn3 = A*((j*beta^2)/(omega*mu)*cos(beta_y*d1) - sigma_ac1*beta_y*sin(beta_y*d1)) - (j*beta^2)/(omega*mu)*B*exp(j*beta_y*d1) - (j*beta^2)/(omega*mu)*C*exp(-j*beta_y*d1) == 0;
eqn4 = D*((j*beta^2)/(omega*mu)*exp(-alpha_y*d2) - sigma_ac2*alpha_y*exp(-alpha_y*d2)) - (j*beta^2)/(omega*mu)*B*exp(j*beta_y*d2) - (j*beta^2)/(omega*mu)*C*exp(-j*beta_y*d2) == 0;
eqns = [eqn1, eqn2, eqn3, eqn4];
vars = [A, B, C, D];
[M,~] = equationsToMatrix(eqns,vars);
simplify(det(M))

ans = 

[Ahmad Awada]

Thank you for your time, I appreciate it.

[David Goodmanson]

Hi Ahmad,

Your eqn1 has a factor of j multiplying B and C. The same equation, eqn. 10 in the paper you included in a comment, does not have that factor. Could you comment on this? It's odd because it looks like including the factor of j leads to the correct answer and leaving it out does not.

[Ahmad Awada]

Hi David

you're right, it seemes a writing mistake. but my question was something else. how can i get the last equation using Matlab.

[David Goodmanson]

Hi Ahmad,

It's reasonably easy to get the result if the starting point is correct, although syms is only easily useful for part of the task. But of course it won't happen with a bad starting point. Did you just delete the image of the equations in the paper? Could you bring it back?

Take a look eqns 1 and 3, and the parallel equations 10 and 12 from the paper. The reason I'm asking about the paper is that if memory serves, eqn.10 was the same as your eqn.1 except that it did not have factors of j on the right hand side. I think (3) and (12) were the same. Not worrying about constants and signs and the terms involving sigma, looking only factors of j, I think the situation is

A =  jB+jC    (1)
jA = jB+jC    (3)

A =  B+C     (10)
jA = jB+jC   (12)

In your equations, (1) has a relative phase of j between A and (B&C). (3) has no such relative phase. A relative phase j in one of those eqns, and no such phase in the other eqn, leads to the correct result for the left hand side of the answer. In the paper, neither of the two equations has a relative phase of j between A and (B&C). This leads to an incorrect result. Is there a mistake in the paper? Am I remembering this stuff wrong?

p.s. Using the eqns

A =  B+C
A = jB+jC

(equivalent to yours, and the algebra and physics make sense) gives the correct result for the left hand side of the answer.

[Ahmad Awada]

Hi

My problem was not if the paper has a mistake or not. My question was just as follows: I have 4 equations with 4 constants, can we solve it and elimenate these constants using Matlab.

If you can help in this, i will appreciate that.

Thanks

Answer 1

[David Goodmanson]

Hi Ahmad,

If the answers don't have to match up with the paper, that makes things easier. I took four equations very similar to the ones in the question and redefined and rescaled some things. Using w instead of omega,

w^2 = 1/(mu*eps)*beta^2     % 1/(mu*eps) is the wave speed (phase velocity)                   

In eqns 3 & 4 you can replace beta^2/(mu*w) with (w*eps), then divide the equations by (w*eps) and define lambda as below to arrive at the four equations used in the code below.

Symbolic variables help to some extent. As far as the substitutions and divisions, I am not an expert with syms, but most of the time I find it faster just to use pencil and paper rather than frustratingly trying to wrestle syms into doing what you want. If someone can incorporate those steps into syms I will be glad to see that done. Toward the end of the calculation there is the division cos(th1)/sin(th1) and syms has to be able to provide cot(th1) somehow, right?

syms th1
simplify(cos(th1)/sin(th1))
ans = cos(th1)/sin(th1)

I must be missing something.

As to the calculation, if you swap eqns 2 and 3 (leaving their names the same), you have a system that looks like

[a b c 0]     [A]
[d e f 0]  x  [B]         = 0
[0 g h k]     [C]
[0 m n q]     [D]

where the letters denote nonzero quantities. Since the ABCD column vector is nonzero, the determinant of the matrix must be zero. That requirement determines ABCD within an overall constant. Take a look at eqns. 1 and 3, which do not involve D.

[a b c]     [A]
[d e f]  x  [B]         = 0             (1)
            [C]

The idea is to determine the ratio R = C/B, which does not depend on the overall scaling of the ABC column vector. There are at least three ways to do this, but I'll use calculation of the null vector, the nonzero vector that is perpendicular to both of the row vectors [a b c] and [d e f] and therefore satisfies (1). This R is the left hand side of the answer.

Now calculate R = C/B with eqns 2 and 4 and the BCD column vector, which gives the right hand side of the answer.

[g h k]     [B]
[m n q]  x  [C]         = 0
            [D]

So:

% th1 = beta_y*d1;
% th2 = beta_yy*d2;
% ayd2 = alpha_y*d2;
% lambda1 = sigma_ac1/(omega*eps);
% lambda2 = sigma_ac2/(omega*eps);
% compact notation:
% w = omega;
% ay = alpha_y;
% by = beta_y;
syms A B C D ay by th1 th2 ayd2 lambda1 lambda2
eqn1 = A*sin(th1) == B*exp(j*th1) + C*exp(-j*th1);
eqn2 = D*ay*exp(-ayd2) == j*by*B*exp(j*th2) -j*by*C*exp(-j*th2);
eqn3 = A*(cos(th1) + j*by*lambda1*sin(th1)) == j*B*exp(j*th1) - j*C*exp(-j*th1);
eqn4 = D*(j*exp(-ayd2) - ay*lambda2*exp(-ayd2)) == j*B*exp(j*th2) + j*C*exp(-j*th2);
eqns = [eqn1,eqn3];
vars = [A,B,C];
[M13,~] = equationsToMatrix(eqns,vars)
n13 = null(M13);             % 3x1 vector
R_left = n13(3)/n13(2)
eqns = [eqn2,eqn4];
vars = [B,C,D];
[M24,~] = equationsToMatrix(eqns,vars);
n24 = null(M24);             % 3x1 vector
R_right = n24(2)/n24(1)

M13 = 
[                         sin(th1),    -exp(th1*1i),   -exp(-th1*1i)]
[cos(th1) + by*lambda1*sin(th1)*1i, -exp(th1*1i)*1i, exp(-th1*1i)*1i]
 
R_left =
-(exp(th1*2i)*(cos(th1) - sin(th1)*1i + by*lambda1*sin(th1)*1i))/ ...
    (cos(th1) + sin(th1)*1i + by*lambda1*sin(th1)*1i)
 
M24 =
[-by*exp(th2*1i)*1i, by*exp(-th2*1i)*1i,                           ay*exp(-ayd2)]
[   -exp(th2*1i)*1i,   -exp(-th2*1i)*1i, - ay*lambda2*exp(-ayd2) + exp(-ayd2)*1i]
 
R_right =
-(exp(th2*2i)*(ay + ay*by*lambda2 - by*1i))/(ay - ay*by*lambda2 + by*1i)

R_left contains the factor exp(th1*2i), which if it is taken over to R_right leads to the exponential factor involving delta d.

Dividing numerator and denominator of R_left by i*sin(th1) gives the answer.

Dividing numerator and denominator of R_right by ay gives the answer.