code assistance

Hi,

I wrote the following code to solve an equation ,to find lambda that nulls the equation (Setting the equation to zero).But when i run the code,Lambda=[],i don't why.the code is:

    %%%%%%%%FIND LAMBDA AND EIGEN VALUE

    load saved_data;
          theta=pi/2;
         zeta=cos(theta);
          I=eye(n,n);
          Q=zeta*I-p*p';

 %T is a matrix(5,5)

          T=M.^(1/2)*Q*M.^(1/2);
          format short e;

 %find the eigen values

          E=eig(T);

%find the negative eigen values

          ind=find(E<0);
          G=E(ind);

%find the smallest negative eigen value

          gamma=min(abs(G));

%find lambda which is in the interval[0,-1/gamma]

          for lambda=0:-1/gamma;
              Q=zeta*I-p*p';
             W=inv(M)+lambda.*Q;
            finsym(-zm'*inv(M)*inv(W)*Q*inv(W)*inv(M)*zm);

%The equation which will be solved to find the value of lambda 
%which null the equation (Setting the equation to zero

            solve(-zm'*inv(M)*inv(W)*Q*inv(W)*inv(M)*zm);
          end

[Naz]

Use {code} button to properly post the script

[zayed]

Excuse me,but i don't know what code button,where i can find it

[Naz]

Check this out in help file: [V,D] = eig(A)

[Walter Roberson]

When you are in the editor, there is a block of 7 buttons at the top immediately under "* Body" and immediately over the box in which you can enter code. The right-most of those buttons is labeled "{} Code". For a demonstration of how to use that button effectively, please see http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup#answer_18099

[zayed]

I have found the eigen values using E=eig(T)

Answer 1

[Walter Roberson]

I do not recognize "finsym". There is an older symbolic toolbox function named "findsym", but it would be of no practical use unless you either displayed the output or used the function form of it and assigned the output to something.

The symbolic toolbox has a "solve" function. As is the case for all functions used in "function / command" equivalence, if you use the command form of it, then each (whitespace-delimited) string that follows on the command line is passed as a string input to the function. For example,

addpath C:\system x86\MATLAB

would be equivalent to calling

addpath('C:\system', 'x86\MATLAB')

Thus, invoking solve in the way you have invoked it would be equivalent to invoking

solve('''-zm''*inv(M)*inv(W)*Q*inv(W)*inv(M)*zm''');

Notice that because you used the command form, the quotes you provided become a literal part of the input.

Even if you had not provided the quotes yourself and had invoked

solve -zm'*inv(M)*inv(W)*Q*inv(W)*inv(M)*zm

then this would be equivalent to calling

solve('-zm''*inv(M)*inv(W)*Q*inv(W)*inv(M)*zm')

However because you are passing a string in to solve, solve() will only know the values of any variables which you have defined in the symbolic toolbox itself (not at the MATLAb level!). You have not defined any values at the level of the symbolic toolbox itself, so you are asking solve to solve the equation completely symbolically, which it is not able to do.

So what should you be doing? Well, you could be invoking,

solve(-zm'*inv(M)*inv(W)*Q*inv(W)*inv(M)*zm)

and then the expression will be evaluated at the MATLAB level and the result of the expression will be passed to solve(). This is generally a Good Thing To Do.

Now, I said generally a good thing to do, not always. If you go back through your code, you will see that at that point in your code, you have no undefined variables, so the expression will either evaluate to 0 or it will not evaluate to 0, and there is nothing for solve() to do.

You could be keeping track of the minimum (absolute) value over the loop, or you could be using fzero() or fsolve().

If you do keep track of the minimum value over the loop, you will find that your loop does not produce any solutions at all. You calculate gamma = min(abs(G)) so you know gamma is a non-negative real number. You then loop for lambda = 0:-1/gamma . Regardless of whether gamma is less than 1 or greater than 1, unless gamma is infinity, -1/gamma will be negative, leading to a loop request starting from 0, incrementing by 1, and ending at a negative number. That's an empty loop. (If gamma is +infinity, then the loop would be 0:0 which would be done once with lambda=0)

Please reconsider your code.

By the way, based upon your previous postings, I don't think you want

gamma = min(abs(G));

You wanted the smallest eigenvalue, and you specifically clarified that -5 is smaller than -1/2, which would correspond to

gamma = min(G);

with further clarification from you needed if G contains complex values.

If you are looking for the entry in G which has the smallest absolute value, then you need to use

[gamma, idx] = min(abs(G));
gamma = G(idx);

This could then be positive or negative; a negative value at least gives the possibility that -1/gamma would be positive (and at least 1) so your loop could proceed as written... but then one would expect that you would want to have code to deal with the possibility that the minimum value happened to be positive so -1/gamma was negative...

=====

Addition, possible modified code. Note that I did not fix the selection of eigenvalues, as you need to clarify that.

[Modified Nov 28 2011, 14:34 - p parameter, Q optimization] [Modified Nov 28 2011, 14:53 - conjugate transpose giving syntax problem]

function lambda = drive_zr17t9(p)

    format short e;

    load saved_data;

    theta = pi/2;
    zeta = cos(theta);
    I = eye(n,n);
    Q = zeta*I-p*p';

    %T is a matrix(5,5)
    Mroot = M.^(1/2);  %optimization
    T = Mroot*Q*Mroot;

    %find the eigen values

    E = eig(T);

    %find the negative eigen values -- this section still needs to be fixed

    ind = find(E<0);
    G = E(ind);

    %find the smallest negative eigen value

    gamma = min(abs(G));  %this is almost certainly wrong!

    %now solve for lambda

    bounds = sort([0, -1/gamma]);  %in case gamma is positive
    Minv = inv(M);  %optimization

    lambda = fzero(@(lambda) zr17t9(lambda,  Minv, Q, zm), bounds);  %do the searching
  end

 function r = zr17t9(lambda, Minv, Q, zm)
    Winv = inv(mInv+lambda.*Q);
    r = -ctranspose(zm)*Minv*Winv*Q*Winv*Minv*zm;
  end

[Walter Roberson]

Sorry, I didn't notice before that you were only examining negative eigenvalues. But you still have the problem that min(abs(G)) is positive and that -1 divided by that will be negative.

[zayed]

First of all i want to thank you because you clarified to me a lot of things.Regarding to the variable LAMBDA i have used (findsyn) command to define the variable,also i have used gamma = min(abs(G));because G contains complex values.But i need more explanation ,a bout invoking,solve(-zm'*inv(M)*inv(W)*Q*inv(W)*inv(M)*zm),and then the expression will be evaluated at the MATLAB level and the result of the expression will be passed to solve(),using fzero() or fsolve().which is the better and easy.What confused me that -1/gamma is -2.1162e+003(negative value).

[Walter Roberson]

findsym() only *reports* on symbolic variables that are part of an expression, and does not *define* any symbolic variables. For example if you had

z=9;

syms x y

findsym(z^3+15*x+exp(x^2)*pi)

then the findsym would reply with the symbolic variable x because that is the only symbolic variable that appears in the expression.

====

abs() of a complex number is always a non-negative number, so watch out for that.

You used G = E(find(E<0)) . I cannot locate an appropriate MATLAB reference page for that at the moment, but I do find some threads from 2001 that say that when complex numbers are compared using < or > that only the real parts are compared. If that is your intent, then it would be clearer if you added the real() call in specifically, as not many people would know that, and you would have problems if MathWorks ever changed the meaning.

===

Your expression -zm'*inv(M)*inv(W)*Q*inv(W)*inv(M)*zm has every variable defined numerically, so the expression evaluates to a double precision number. The value you are trying to seek, lambda, you have given a numeric value by the "for lambda" loop. The expression is just a number, and solve() applied to a number can be satisfied only if the number is exactly zero.

Using fzero() or fsolve() to solve the equation is, if I recall properly, a question you have already posted.

[zayed]

So how can i find the smallest negative eigen value(complex numbers).Also i do the loop because LAMBDA should be in the interval [0 ,-1/gamma].It's the question that i have been posted and i tried to write a code, but i don't have any result.If i use fzero() or fsolve(), how can i do it.

[Walter Roberson]

When you are looking for "negative" eigenvalues, do you want to examine only the real part or only the imaginary part?

The code you have posted has no LAMBDA, only lambda .

The line

for lambda = 0:-1/gamma

does NOT express that lambda should range over the interval [0,-1/gamma] . Instead, it expresses that lambda should first be set to 0, then to 1, then to 2, and so on, ending (loop body _not_ executed) the first time that lambda is greater than -1/gamma . If gamma is positive and is not infinity, then -1/gamma will be negative and since 0 is greater than any negative value, the loop body would not be executed at all.

You can do things such as

for lambda = linspace(0, -1/gamma, 100)

and that would use 100 equidistant points with endpoints 0 and -1/gamma (thus, 98 interior points.) linspace() is able to deal with the second endpoint being greater than the first endpoint: in such a situation, the list of points would be in descending order, but they would all be there.

See Answer for possible modified code (we can't format code in comments)

[zayed]

I changed the loop and it's work but it give me warning that list of equations is empty.by the way i have anew piece of information that the equation monotonically increasing and has a unique zero,

computationally efficient techniques such as the Newton method

can be used to find lambda which is real value.

[Walter Roberson]

The message about the list of equations being empty is coming from solve(), isn't it? I told you that wouldn't work...

fzero() documentation says, "The algorithm, which was originated by T. Dekker, uses a combination of bisection, secant, and inverse quadratic interpolation methods."

That is probably fairly efficient.

[Walter Roberson]

Code modified to take p as a parameter, and to optimize calculation of Q.

[Walter Roberson]

Code modified to avoid using ' operator as MATLAB was interpreting line as a command (?!)