this is my sample code and i have to plot 8 by 8 matrix in my actual code, the code itself is correct but it is taking forever,any other way to plot this please help

%finding k for different values of w
syms k
c44=44;
i=sqrt(-1);
w=1:3;
y1=zeros(1,length(w));
y2=zeros(1,length(w));
s=sym(zeros(1,length(w)));
p=sym(zeros(1,length(w)));
q=sym(zeros(1,length(w)));
y=sym(zeros(1,length(w)));
for j=1:length(w)
s(j)=c44*exp(k)*w(j);
p(j)=c44/w(j);
q(j)=c44+w(j)*i/k^2;
A=[c44*s(j) 2;p(j) q(j)];
z=det(A);
y(j)=vpasolve(z,k);
y1=real(y);
y2(j)=w(j)/y1(j);
y3=y2/c44;
disp(y3)
plot(y2,w)

end

   -0.0033         0         0

   -0.0033   -0.0055         0

   -0.0033   -0.0055    0.3693

Answer 1

[Torsten]

w=1:3;
c44 = 44;
n = numel(w);
y = zeros(n,1);
y1 = zeros(n,1);
y2 = zeros(n,1);
y3 = zeros(n,1);
for j = 1:numel(w)
  y(j) = fsolve(@(k)fun(k,c44,w(j)),1,optimset('Display','none'));
  y1(j) = real(y(j));
  y2(j) = w(j)/y1(j);
  y3(j) = y2(j)/c44;
end
disp(y3)

    0.2131
    0.3015
    0.3693

plot(y2,w)

function res = fun(k,c44,w)
  s=c44*exp(k)*w;
  p=c44/w;
  q=c44+w*1i/k^2;
  A=[c44*s 2;p q];
  res=det(A);
end

Answer 2

[Walter Roberson]

There isn't really much you can do.

Your z is the sum of exponentials, and inherently involves complex numbers; there just is not reasonable way to calculate a closed-form solution. The solution, k, is going to be complex valued and approximate. There is no good logic to reduce the calculations.

The other part of your concern is that you are dealing with an 8 x 8 symbolic matrix and taking the determinant of it. The determinant of a symbolic matrix can take a fair bit of time to calculate. Here is one way that you can speed up the process:

%z=det(A);
M = sym('M', size(A));
dM = det(M);     %determinant of the replacement matrix not the original
z = subs(dM, M, A);    %substitute actual coefficients for representative coefficients

That is, you create a matrix of just variables the same size as the matrix you want to take the determinant of, take the determinant of that matrix, getting out a whole long sum like M11*M22*M34... + M14*M21*M35... which is the closed form for a determinant of a matrix that size; then you substitute the original terms for the place-holders.

Because of the way that MATLAB stores symbolic expressions internally, for larger matrices this can speed up calculation of the determinant quite a bit.

You can calculate M and dM once before the loop, just substituting in the actual coefficients within the loop.

[Odette]

Im sorry but i dont know how to do this.. also is there any way i can plot without actually finding k and taking it more like a function

[Walter Roberson]

format long g

%finding k for different values of w
syms k
c44=44;
i=sqrt(-1);
w=1:3;
y1=zeros(1,length(w));
y2=zeros(1,length(w));
s=sym(zeros(1,length(w)));
p=sym(zeros(1,length(w)));
q=sym(zeros(1,length(w)));
y=sym(zeros(1,length(w)));

sizeA  = [2 2];   %change this if you make A larger
M = sym('M', sizeA);
dM = det(M);

for j=1:length(w)
    s(j)=c44*exp(k)*w(j);
    p(j)=c44/w(j);
    q(j)=c44+w(j)*i/k^2;
    A=[c44*s(j) 2;p(j) q(j)];
    z = subs(dM, M, A);
    y(j)=vpasolve(z,k);
    y1=real(y);
    y2(j)=w(j)/y1(j);
end

y3=y2/c44;
disp(y3)

      -0.00330567344263619       -0.0055019546239557         0.369253916790257

plot(w,y2)

vpa(y(:),10)

ans = 

[Walter Roberson]

also is there any way i can plot without actually finding k and taking it more like a function

Your determinant z is an inherently complex expression in k. If you want to plot it against k to visualize where the zeros are, then you need to plot using an independent real axes and complex axes, and a dependent axes that is the value of the determinent -- a surface plot rather than a line plot. And you have to have a good idea of where the zeros are in order to be able to get a good plot range.

Would you want to plot all length(w) surfaces in the same plot, or would you want one surface plot for each w value?

Note that if you do not find k numerically, then you will not be able to find the y and y1 and y2 values.

[Odette]

Then i think i have to find k numerically.. i tried the above program but in my actual program it shows error in y(j)=vpasolve(z,k). Unable to perform assignment because the left and right sides have a different number of elements.. please do reply again

[Walter Roberson]

In your actual program change

    y(j)=vpasolve(z,k);
    y1=real(y);
    y2(j)=w(j)/y1(j);

to

    yj = vpasolve(z,k);
    if numel(yj) ~= 1
        fprintf('note: iteration #%d generated %d results, using first or NaN\n', j, numel(yj));
    end
    if isempty(yj); yj = sym(NaN); end
    y(j) = yj(1);
    
    y1(j) = real(y(j));
    y2(j)=w(j)/y1(j);

[Odette]

Now the loop is stuck i think after first iteration .. its showing iteration #1 generated 0 results ,then won't show the next iteration it keeps running

[Walter Roberson]

Well, that is certainly possible.

vpasolve() uses a modified Newton-Raphson approach. Newton-Raphson can get stuck on asymtopes descending a long slope without realizing that you have to head up the slope in order to find a better dip on the other side. Or it can just get stuck on messy bumpy surfaces -- especially ones that contain strongly non-linear terms such as exp()

@Torsten shows the approach that would be used numerically, using fsolve() -- and note for this purpose that fzero() will not work as fzero() cannot search over complex values.

You can also take the z and do something like

zfun = matlabFunction(z);
sol = fsolve(zfun, GUESS)

and that has the potential to be much faster... but if your expression involves fine balancing of competing non-linear terms then the loss of precision from going to a numeric function can be devastating.

Sometimes you are better off dividing the input, k, into real and imaginary components, and doing a ga() or pso() search over the norm of the square of z; inside the function you would reconstruct k by k = kr + 1i*ki; if you name the parameter kri then that would be complex(kri(1), kri(2))

The norm of the square of z turns out to be real(z)^2 + imag(z)^2 .

The reason to use pso or ga is that those can search a number of different locations simultaneously -- reducing the problem of ending up stuck in a single unproductive region.

You could very likely code for Vectorized calculations, such as k=complex(kri(:,1), kri(:,2)) and the output of matlabFunction is typically vectorized already.

[Odette]

yes it contains a lot of exponential function,and i try to use fsolve -Objective function is returning undefined values at initial point. FSOLVE cannot continue-this keeps showing up even when i change the initial guess

[Torsten]

We cannot find the reason for failure without the actual code where this error message appeared.

[Odette]

c44=43;
e15=-0.48;
e11=7.61;
u11=100;
d11=2.6;
p=5700;
pe=3512;
c44e=533.34;
e11e=5.02;
csh=0.38969499199;
q=1.602*(1/10^(19));
uo=10^(20);
i=sqrt(-1);
%omega(w) is from 0 to 3
sizeA=[8 8];
M=sym('M',sizeA);
dM=det(M);


w=0:0.2:3;
n=length(w);

 %preallocation
c1=zeros(1,n);
 z1=zeros(1,n);
s=sym(zeros(1,n));
s3=sym(zeros(1,n));
b=sym(zeros(1,n));
c=sym(zeros(1,n));
y=sym(zeros(1,n));
z=sym(zeros(1,n));
s4=sym(zeros(1,n));
s5=sym(zeros(1,n));
s6=sym(zeros(1,n));
p3=sym(zeros(1,n));
p4=sym(zeros(1,n));
p5=sym(zeros(1,n));
p6=sym(zeros(1,n));
q3=sym(zeros(1,n));
q4=sym(zeros(1,n));
q5=sym(zeros(1,n));
q6=sym(zeros(1,n));
y1=zeros(1,n);
y2=zeros(1,n);
sizeA=[8 8];
M=sym('M',sizeA);
dM=det(M);
%for loop for w

 for j=1:n
 syms k
 s(j)=sqrt(1-(w(j)/k*csh)^2);
a=d11*(c44*e11+e15^2);
b(j)=e11*d11*p*w(j)^2/k^2+w(j)*(c44*e11+e15^2)*i/k^2-c44*q*uo*u11/k^2;
c(j)=(e11*w(j)*i/k^2-q*uo*u11/k^2)*p*w(j)^2/k^2;
s3(j)=sqrt(1+(-b(j)-sqrt(b(j)^2-4*a*c(j)))/2*a);
s4(j)=-s3(j);
s5(j)=sqrt(1+(-b(j)+sqrt(b(j)^2-4*a*c(j)))/2*a);
s6(j)=-s5(j);
p3(j)=-e15*(s3(j)^2-1)/(c44*(s3(j)^2-1)+p*w(j)^2/k^2);
p4(j)=-e15*(s4(j)^2-1)/(c44*(s4(j)^2-1)+p*w(j)^2/k^2);
p5(j)=-e15*(s5(j)^2-1)/(c44*(s5(j)^2-1)+p*w(j)^2/k^2);
p6(j)=-e15*(s6(j)^2-1)/(c44*(s6(j)^2-1)+p*w(j)^2/k^2);
q3(j)=-uo*u11*(s3(j)^2-1)/(d11*(s3(j)^2-1)+w(j)*i/k^2);
q4(j)=-uo*u11*(s4(j)^2-1)/(d11*(s4(j)^2-1)+w(j)*i/k^2);
q5(j)=-uo*u11*(s5(j)^2-1)/(d11*(s5(j)^2-1)+w(j)*i/k^2);
q6(j)=-uo*u11*(s6(j)^2-1)/(d11*(s6(j)^2-1)+w(j)*i/k^2);

%first row
a11=e15*exp(k);a12=e15*exp(-k);
a13=(c44*p3(j)+e15)*exp(s3(j)*k);
a14=(c44*p4(j)+e15)*exp(s4(j)*k);
a15=(c44*p5(j)+e15)*exp(s5(j)*k);
a16=(c44*p6(j)+e15)*exp(s6(j)*k);a17=0;a18=0;

%second row
a21=-e11*exp(k);a22=-e11*exp(-k);
a23=(e15*p3(j)-e11)*exp(s3(j)*k);
a24=(e15*p4(j)-e11)*exp(s4(j)*k);
a25=(e15*p5(j)-e11)*exp(s5(j)*k);
a26=(e15*p6(j)-e11)*exp(s6(j)*k);a27=0;a28=0;

%third row
a31=uo*u11*exp(k);a32=uo*u11*exp(-k);
a33=exp(s3(j)*k)*(uo*u11-d11*q3(j));
a34=exp(s4(j)*k)*(uo*u11-d11*q4(j));
a35=exp(s5(j)*k)*(uo*u11-d11*q5(j));
a36=exp(s6(j)*k)*(uo*u11-d11*q6(j));
a37=0;a38=0;

%fourth row
a41=uo*u11;a42=uo*u11;a43=uo*u11*-d11*q3(j);
a44=uo*u11-d11*q4(j);
a45=uo*u11-d11*q5(j);
a46=uo*u11-d11*q6(j);a47=0;a48=0;

%fifth row
a51=e15;a52=e15;a53=c44*p3(j)+e15;
a54=c44*p4(j)+e15;
a55=c44*p5(j)+e15;
a56=c44*p6(j)+e15;a57=-c44e;a58=0;

%sixth row
a61=0;a62=0;a63=p3(j)/s3(j);a64=p4(j)/s4(j);a65=p5(j)/s5(j);a66=p6(j)/s6(j);a67=0;a68=0;

%seventh row
a71=-e11+e11e;a72=-e11+e11e;a73=(e15*p3(j)-e11)+e11e;
a74=(e15*p4(j)-e11)+e11e;
a75=(e15*p5(j)-e11)+e11e;a76=(e15*p6(j)-e11)+e11e;
a77=0;a78=e11e;

%eigth row
a81=0;a82=0;a83=0;a84=0;a85=0;a86=0;a87=0;a88=-1;

A=[a11 a12 a13 a14 a15 a16 a17 a18 ;
    a21 a22 a23 a24 a25 a26 a27 a28;
    a31 a32 a33 a34 a35 a36 a37 a38 ;
    a41 a42 a43 a44 a45 a46 a47 a48;
    a51 a52 a53 a54 a55 a56 a57 a58;
    a61 a62 a63 a64 a65 a66 a67 a68;
    a71 a72 a73 a74 a75 a76 a77 a78;
    a81 a82 a83 a84 a85 a86 a87 a88];
z=subs(dM,M,A);
zfun=matlabFunction(z);

y= fsolve(zfun,1);
        
  y1(j) = real(y(j));
    y2(j)=w(j)/y1(j);
 end

Error using fsolve
Objective function is returning undefined values at initial point. FSOLVE cannot continue.

 c2=y2/csh;
plot(c2,w)

 
 

[Odette]

this is my actual code i know it is very lengthy but please reply

[Walter Roberson]

At least for the first iteration, A is not full rank, so its determinant z comes out equal to 0 (but written in a complicated way that leads to a division by 0 if you use it naively)

c44=43;
e15=-0.48;
e11=7.61;
u11=100;
d11=2.6;
p=5700;
pe=3512;
c44e=533.34;
e11e=5.02;
csh=0.38969499199;
q=1.602*(1/10^(19));
uo=10^(20);
i=sqrt(-1);
%omega(w) is from 0 to 3
sizeA=[8 8];
M=sym('M',sizeA);
dM=det(M);


w=0:0.2:3;
n=length(w);

 %preallocation
c1=zeros(1,n);
 z1=zeros(1,n);
s=sym(zeros(1,n));
s3=sym(zeros(1,n));
b=sym(zeros(1,n));
c=sym(zeros(1,n));
y=sym(zeros(1,n));
z=sym(zeros(1,n));
s4=sym(zeros(1,n));
s5=sym(zeros(1,n));
s6=sym(zeros(1,n));
p3=sym(zeros(1,n));
p4=sym(zeros(1,n));
p5=sym(zeros(1,n));
p6=sym(zeros(1,n));
q3=sym(zeros(1,n));
q4=sym(zeros(1,n));
q5=sym(zeros(1,n));
q6=sym(zeros(1,n));
y1=zeros(1,n);
y2=zeros(1,n);
sizeA=[8 8];
M=sym('M',sizeA);
dM=det(M);
%for loop for w

 for j=1:n
 syms k
 s(j)=sqrt(1-(w(j)/k*csh)^2);
a=d11*(c44*e11+e15^2);
b(j)=e11*d11*p*w(j)^2/k^2+w(j)*(c44*e11+e15^2)*i/k^2-c44*q*uo*u11/k^2;
c(j)=(e11*w(j)*i/k^2-q*uo*u11/k^2)*p*w(j)^2/k^2;
s3(j)=sqrt(1+(-b(j)-sqrt(b(j)^2-4*a*c(j)))/2*a);
s4(j)=-s3(j);
s5(j)=sqrt(1+(-b(j)+sqrt(b(j)^2-4*a*c(j)))/2*a);
s6(j)=-s5(j);
p3(j)=-e15*(s3(j)^2-1)/(c44*(s3(j)^2-1)+p*w(j)^2/k^2);
p4(j)=-e15*(s4(j)^2-1)/(c44*(s4(j)^2-1)+p*w(j)^2/k^2);
p5(j)=-e15*(s5(j)^2-1)/(c44*(s5(j)^2-1)+p*w(j)^2/k^2);
p6(j)=-e15*(s6(j)^2-1)/(c44*(s6(j)^2-1)+p*w(j)^2/k^2);
q3(j)=-uo*u11*(s3(j)^2-1)/(d11*(s3(j)^2-1)+w(j)*i/k^2);
q4(j)=-uo*u11*(s4(j)^2-1)/(d11*(s4(j)^2-1)+w(j)*i/k^2);
q5(j)=-uo*u11*(s5(j)^2-1)/(d11*(s5(j)^2-1)+w(j)*i/k^2);
q6(j)=-uo*u11*(s6(j)^2-1)/(d11*(s6(j)^2-1)+w(j)*i/k^2);

%first row
a11=e15*exp(k);a12=e15*exp(-k);
a13=(c44*p3(j)+e15)*exp(s3(j)*k);
a14=(c44*p4(j)+e15)*exp(s4(j)*k);
a15=(c44*p5(j)+e15)*exp(s5(j)*k);
a16=(c44*p6(j)+e15)*exp(s6(j)*k);a17=0;a18=0;

%second row
a21=-e11*exp(k);a22=-e11*exp(-k);
a23=(e15*p3(j)-e11)*exp(s3(j)*k);
a24=(e15*p4(j)-e11)*exp(s4(j)*k);
a25=(e15*p5(j)-e11)*exp(s5(j)*k);
a26=(e15*p6(j)-e11)*exp(s6(j)*k);a27=0;a28=0;

%third row
a31=uo*u11*exp(k);a32=uo*u11*exp(-k);
a33=exp(s3(j)*k)*(uo*u11-d11*q3(j));
a34=exp(s4(j)*k)*(uo*u11-d11*q4(j));
a35=exp(s5(j)*k)*(uo*u11-d11*q5(j));
a36=exp(s6(j)*k)*(uo*u11-d11*q6(j));
a37=0;a38=0;

%fourth row
a41=uo*u11;a42=uo*u11;a43=uo*u11*-d11*q3(j);
a44=uo*u11-d11*q4(j);
a45=uo*u11-d11*q5(j);
a46=uo*u11-d11*q6(j);a47=0;a48=0;

%fifth row
a51=e15;a52=e15;a53=c44*p3(j)+e15;
a54=c44*p4(j)+e15;
a55=c44*p5(j)+e15;
a56=c44*p6(j)+e15;a57=-c44e;a58=0;

%sixth row
a61=0;a62=0;a63=p3(j)/s3(j);a64=p4(j)/s4(j);a65=p5(j)/s5(j);a66=p6(j)/s6(j);a67=0;a68=0;

%seventh row
a71=-e11+e11e;a72=-e11+e11e;a73=(e15*p3(j)-e11)+e11e;
a74=(e15*p4(j)-e11)+e11e;
a75=(e15*p5(j)-e11)+e11e;a76=(e15*p6(j)-e11)+e11e;
a77=0;a78=e11e;

%eigth row
a81=0;a82=0;a83=0;a84=0;a85=0;a86=0;a87=0;a88=-1;

A=[a11 a12 a13 a14 a15 a16 a17 a18 ;
    a21 a22 a23 a24 a25 a26 a27 a28;
    a31 a32 a33 a34 a35 a36 a37 a38 ;
    a41 a42 a43 a44 a45 a46 a47 a48;
    a51 a52 a53 a54 a55 a56 a57 a58;
    a61 a62 a63 a64 a65 a66 a67 a68;
    a71 a72 a73 a74 a75 a76 a77 a78;
    a81 a82 a83 a84 a85 a86 a87 a88]

rank(A)

z = subs(dM,M,A);

simplify(z)

zfun=matlabFunction(z);

   y0 = 1;
   zfun(y0)  %what are we getting back?
   limit(z, k, y0)
   
   y = fsolve(zfun, y0);
        
  y1(j) = real(y(j));
    y2(j)=w(j)/y1(j);
 end

A = 

ans = 7

ans = 

0

ans = NaN

Error using symengine
First argument is not defined under the assumption that 'k' has the property 'Dom::Interval(9999999/10000000, 1) union Dom::Interval(1, 10000001/10000000)'.

Error in mupadengine/evalin_internal

Error in mupadengine/fevalHelper

Error in mupadengine/feval_internal

Error in sym/limit (line 62)
rSym = feval_internal(symengine, 'symobj::map', args{1}.s, 'symobj::limit', args{2}.s, args{3}.s, dir);

 c2=y2/csh;
plot(c2,w)