solve 4 equations with some unknown coefficients.

7 次查看（过去 30 天）
Ahmad Awada
Ahmad Awada 2023-6-10
编辑： David Goodmanson 2023-6-17
Ran in:
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
  10 个评论
显示 8更早的评论
David Goodmanson
David Goodmanson 2023-6-14
编辑：David Goodmanson 2023-6-14
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
Ahmad Awada 2023-6-15
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

请先登录，再进行评论。

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

回答（1 个）

David Goodmanson
David Goodmanson 2023-6-17
编辑：David Goodmanson 2023-6-17
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.

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by