I am trying to solve this problem on MATLAB of which code i have shared. It involves solving eigen values of a symbolic matrix with one variable 'f' and then using those values further in calculation. It is taking whole day what should i do?

syms f
l=0.1008;
T=0.005;
Dh=0.01;
D1=0.04;
D2=0.16;
x=20*(T*pi*Dh^2/4)/(0.1008*pi*D1^2/4);
c_0=619;
k=2*pi*f/c_0;
t=(0.006+1i*k*(T+0.75*Dh))/x;
Ka=(k^2)-((1i*4*k)/(D1*t));
Kb=(k^2)-((1i*4*k*D1)/((D2^2-D1^2)*t));
Delta=[-1 0 0 0; 0 -1 0 0; 0 0  Ka^2 Ka^2-k^2; 0 0 Kb^2-k^2 Kb^2]
g=eig(Delta);
ee=vpa(g,3)
C1=[1 1 1 1];
C2=[-((g(1,1)^2+Ka^2)/(Ka^2-k^2)) -((g(2,1)^2+Ka^2)/(Ka^2-k^2)) -((g(3,1)^2+Ka^2)/(Ka^2-k^2)) -((g(4,1)^2+Ka^2)/(Ka^2-k^2))];
C3=[1/g(1,1) 1/g(2,1) 1/g(3,1) 1/g(4,1)]; 
C4=[C2(1,1)/g(1,1) C2(1,2)/g(2,1) C2(1,3)/g(3,1) C2(1,4)/g(4,1)];
A1=[C3(1,1)*exp(g(1,1)*l) C3(1,2)*exp(g(2,1)*l) C3(1,3)*exp(g(3,1)*l) C3(1,4)*exp(g(4,1)*l)];
A2=[C4(1,1)*exp(g(1,1)*l) C4(1,2)*exp(g(2,1)*l) C4(1,3)*exp(g(3,1)*l) C4(1,4)*exp(g(4,1)*l)];
A3=[-exp(g(1,1)*l)/(i*k) -exp(g(2,1)*l)/(i*k) -exp(g(3,1)*l)/(i*k) -exp(g(4,1)*l)/(i*k)];
A4=[-(exp(g(1,1)*l)*C2(1,1))/(i*k) -(exp(g(1,1)*l)*C2(1,2))/(i*k) -(exp(g(1,1)*l)*C2(1,3))/(i*k) -(exp(g(1,1)*l)*C2(1,4))/(i*k)];
E1=[C3(1,1) C3(1,2) C3(1,3) C3(1,4)];
E2=[C4(1,1) C4(1,2) C4(1,3) C4(1,4)];
E3=[1/(i*k) 1/(i*k) 1/(i*k) 1/(i*k)];
E4=[-C2(1,1)/(i*k) -C2(1,2)/(i*k) -C2(1,3)/(i*k) -C2(1,4)/(i*k)]
A=[A1; A2; A3; A4]
I=eye(4);
invA=A\I;
E=[E1; E2; E3; E4];
PP=E*invA;
vpa(PP,3)

[darova]

Maybe you should try numerical solution instead of symbolic

[Walter Roberson]

Which step is taking all day? The eig step takes roughly 20 seconds for me.

[ABHISHEK KALYANWAT]

if i run the code from line 1(starts with syms f) to line 23 (starts with A4=[) , Matlab takes about 30 min to solve it and if i run the whole code , status bar remains busy all day. what could be the problem ? Is the issue with the code or the system on which I am running MATLAB ON. My sytem information are in the attached image. Code is written below.

if the problem is with code . Please suggest what can be done? @ walter ROBESON @ Darova

Thank you.

syms f
l=0.1008;
T=0.005;
Dh=0.01;
D1=0.04;
D2=0.16;
x=20*(T*pi*Dh^2/4)/(0.1008*pi*D1^2/4);
c_0=619;
k=2*pi*f/c_0;
t=(0.006+1i*k*(T+0.75*Dh))/x;
Ka=(k^2)-((1i*4*k)/(D1*t));
Kb=(k^2)-((1i*4*k*D1)/((D2^2-D1^2)*t));
Delta=[-1 0 0 0; 0 -1 0 0; 0 0  Ka^2 Ka^2-k^2; 0 0 Kb^2-k^2 Kb^2]
g=eig(Delta);
ee=vpa(g,3)
C1=[1 1 1 1];
C2=[-((g(1,1)^2+Ka^2)/(Ka^2-k^2)) -((g(2,1)^2+Ka^2)/(Ka^2-k^2)) -((g(3,1)^2+Ka^2)/(Ka^2-k^2)) -((g(4,1)^2+Ka^2)/(Ka^2-k^2))];
C3=[1/g(1,1) 1/g(2,1) 1/g(3,1) 1/g(4,1)] 
C4=[C2(1,1)/g(1,1) C2(1,2)/g(2,1) C2(1,3)/g(3,1) C2(1,4)/g(4,1)];
A1=[C3(1,1)*exp(g(1,1)*l) C3(1,2)*exp(g(2,1)*l) C3(1,3)*exp(g(3,1)*l) C3(1,4)*exp(g(4,1)*l)];
A2=[C4(1,1)*exp(g(1,1)*l) C4(1,2)*exp(g(2,1)*l) C4(1,3)*exp(g(3,1)*l) C4(1,4)*exp(g(4,1)*l)];
A3=[-exp(g(1,1)*l)/(i*k) -exp(g(2,1)*l)/(i*k) -exp(g(3,1)*l)/(i*k) -exp(g(4,1)*l)/(i*k)];
A4=[-(exp(g(1,1)*l)*C2(1,1))/(i*k) -(exp(g(1,1)*l)*C2(1,2))/(i*k) -(exp(g(1,1)*l)*C2(1,3))/(i*k) -(exp(g(1,1)*l)*C2(1,4))/(i*k)]
E1=[C3(1,1) C3(1,2) C3(1,3) C3(1,4)];
E2=[C4(1,1) C4(1,2) C4(1,3) C4(1,4)];
E3=[1/(i*k) 1/(i*k) 1/(i*k) 1/(i*k)];
E4=[-C2(1,1)/(i*k) -C2(1,2)/(i*k) -C2(1,3)/(i*k) -C2(1,4)/(i*k)]
A=[A1; A2; A3; A4]
I=eye(4);
invA=A\I;
E=[E1; E2; E3; E4];
PP=E*invA;
vpa(PP,3)

i

[Walter Roberson]

Why are you constructing invA and multiplying by E? Why not use A\E?

Note: symbolic eigenvalues are slow and the result is hard to read.

[David Goodmanson]

Hello Abhishek,

A symbolic calculation like this may well take forever and a day, so it seems worthwhile to try the doing this numerically. This means calculating the frequencies one at a time, but that can be useful to scope out the results. It seems that A, the matrix that you want to invert, has huge condition numbers for every f.

The conditon number is a crazy 4.9003e+281 at f = 8, goes infinite a bit above that frequency and has a minimum somewhere around 3e16 in the area around f = 1e-4. The Delta matrix has similar problems although not as extreme.

Anything is possible with vpa, I suppose, if you want to use enough digits and walk away from the machine for as long as it takes (go rafting on the Colorado maybe). But condition numbers like that are crackers and indicate that something is not right with the underlying calculation.

[ABHISHEK KALYANWAT]

Good Day Everyone

I tried solving the problem numerically & found final matrix 'PP' gets to infinty around f=40. It seems the problem is with the calculation. after reviewing the code , i found there may be some issue with formulation of 'Delta' matrix'.

The original Delta matrix to be solved is

Delta=[-1 0 D 0; 0 -1 0 D; D 0  Ka^2 Ka^2-k^2; 0 D Kb^2-k^2 Kb^2];

where D=d/dz (differential operator). z is direction (another variable).

Since no element of matrix 'Delta' has dependency on 'z'. I deduce MATRIX can be written equivalent to

Delta=[-1 0 0 0; 0 -1 0 0; 0 0  Ka^2 Ka^2-k^2; 0 0 Kb^2-k^2 Kb^2];

Since calculation is taking all day ,I am not sure if i am using calculas concept accurately here.

I do not know how to solve eigen vector of a matrix having differential operator like d/dz as metioned above.

Please tell me the way to solve it in its original form.

Thank you

[Walter Roberson]

where D=d/dz (differential operator)

No, if it were the differential operator then D could not appear by itself in the Δmatrix: it would have to be followed by an expression or variable. It would also be quite unlikely to be followed by or similar: if you were taking the derivative of then you would skip the + in the expression.

It looks to me as if D is a variable that is being used linearly in entries in row 3 and 4.

Possibly D is the derivative of something particular, but if so then we would have to know what it is the derivative of.

That equation does not look to me like something you would solve with eigenvalues: that is an equation for simultaneous equations. With the right hand side being all 0, you should probably be thinking about null spaces rather than eigenvalues.

[David Goodmanson]

Hi Abhishek

could you post the entire set of equations for this problem, including the definition of D? It appears that D could be d/dz, since after doing the mattix multiplication delta*y, D would be operating on the components of y.

[ABHISHEK KALYANWAT]

Thank you so much for your comments.

Here is the background of 'Delta Matrix' here notation used for this matrix is :

w=2*pi*f (f is frequency) , Everything else is constant. M1 ,M2 are 0 for this particular case.

are the zeros of the determinant of the coefficient matrix of equation (6a).

I have used 'PP' notation for matrix .

PREVIOUSLY I WAS TRYING TO CALCULATE EIGEN VECTORS TO FIND ZEROS OF THE DETERMINANT OF THE COEFFICIENT MATRIX WHICH IS TOTALLY INCORRECT. EXCUSE ME FOR MY LIMITED UNDERSTANDING OF MATHEMATICS.

  Delta=[-1 0 D 0; 0 -1 0 D; D 0  Ka^2 Ka^2-k^2; 0 D Kb^2-k^2 Kb^2]

CAN YOU TELL ME HOW CODE SHOULD BE WRITTEN TO FIND ZEROS OF THE DETERMINANT OF DELTA MATRIX.

THANK YOU.

Answer 1

[ABHISHEK KALYANWAT]

Thank you so much David ,walter and Darova for your answers. Trying to solve it numerically seems worth trying here.