@JohnS Thanks John, vv,

Here I post my recent findings about my problem.

 The idea is to replace

[M11 0] 
[  0 0]

with

[M11     0] 
[ 0  eps*I]

for some small value eps, small enough to perturb the problem only a small amount. Let N = M^(-1/2), multiply the equation on the left by

[N  0]
[ 0 c*I] 

and insert the identity matrix in the form

[N  0 ]   [N^-1      0]
[0 c*I]   [ 0  (1/c)*I]

between the matrices and the eigenvector

[X     ]

[lambda].

Multiply it all out blockwise and we get what they got, with

I22 = c^2*eps*I. 

You want I22 to be the identity matrix, so

c = 1/sqrt(eps).

c is what they called LV. So how small do we have to make eps to not affect the problem too much? Having no direct experience with this I don't know, but I would hazard a guess of least 10 times smaller than the smallest eigenvalue of M11. Here the smallest eigenvalue is 3e-11 so we can plan accordingly. It does lead to some pretty big LV numbers. (The largest eigenvalue is 4e-3, so scaling could be a problem).

There are spurious eigenvalues for the problem I have. Note that the original problem has a number of eigenvalues "at infinity", which correspond to large negative eigenvalues for the approximating problem. Those eigenvalues have size of order LV. The number of these eigenvalues "at infinity" is the dimension of lambda. So these eigenvalues should be ignored -- only the "small" eigenvalues should be kept. The trouble is that if LV becomes large, the smaller eigenvalues are subject to numerical (roundoff) "noise" of size LV.u, u = unit round-off. So taking LV = 10²⁰⁰ would likely destroy any accuracy the small eigenvalues have left. (I could be wrong; it will depend on the exact algorithm used.).

I am noticed these kinds of problems are solved very easily using Eigenvalues(K,M) command in Maple, because this command accepts both K and M matrices as its inputs, regardless how the K and M are.

But does anybody knows how this command compute the eigenvalues of generalized eigenvalue problem?

@vv Ok let prepare data you need.

M=5(I+1)(J+1),  N=10(I+1)+10(J+1)

For M=500 and N=200, here are the five known omegas:

The related matrices are as below:

Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/K21.mw .
 

K21.mw

K22.mw

K11.mw

K22.mw

M21.mw

M22.mw

M11.mw

M12.mw

Waiting your comments.

@vv Dont you have any comment how to choose LV in this method?

What other ways you know to find eigenvalues of problem in maple, that dont have the form of standard eigenvalue problem like mentioned problem?

@vv This Method are presented in some authentic papers, so I trusted them. Now I am using this method for my own problem.

for example:

http://www.sciencedirect.com/science/article/pii/S0045794916306113

@vv thanks. As you see my matrices are inconsistent (A matrix have very big and very small elements at the same time) . What do you recommend in such conditions?

@mehdibaghaee It seems the problem is because of matrix inconsistency(A matrix have very big and very small elements at the same time).

How cope with this problem?

@vv  when I use my own matrices, the results do not converge and increse.

But whats the problem with my matrices??

@vv any comment yet???

@vv I re arranged my question as above, please help.

Acer deleted my question, guide me to pose my question here.

Let me ask my question in another way. Here are a little introduction:

Doing above method on my problem in Maple, I got several omegas for several LVs(the omegas continuesly increase when LV increases). While I expect my results converge to specific omegas when LV increases.

Download Suals.mw

K21.mw

K22.mw

K11.mw

K12.mw

M21.mw

M22.mw

M11.mw

M12.mw

Where is the problem in my matrices?

Please guide and help me.

@taro using CodeGeneration is tedious for worksheet including 3000 lines of commands. Is there a fast method to translate a Maple worksheet into Matlab?

@taro yes this post was for four years ago. I wanted to know if there is new activity to develope more complete new tools in this case.

Anyway thanks

@tomleslie Ok, lets do the same process for some of my matrices.

Download unassign2-1.mw

But error persists

@vv Let me explain more.

Let the A be a square matrix of function of omega and b is vector function of  y1,y2,y3,...,yi, and X be X := [seq(seq(U[i, j], j = 0 .. JJ), i = 0 .. II), seq(seq(V[i, j], j = 0 .. JJ), i = 0 .. II), seq(seq(W[i, j], j = 0 .. JJ), i = 0 .. II), seq(seq(Phi[i, j], j = 0 .. JJ), i = 0 .. II), seq(seq(Xi[i, j], j = 0 .. JJ), i = 0 .. II)]

I initialize omega with for loop, then solve the  AX=b by X=inv(A)*b. So, X will be function of yi. Now I want to put them(Solved X) into series of equations(functions of X) and build the new CY=0 equation. Then determinant of C will be calculated. The omega that makes the determinant of C zero will be my output.

Note II and JJ are are 50.

It seems if there was a way to unassign X[i], the problem will be solved. But I dont know How is it pussible to unassign X[i].

@vv 

May matrix have the variables as below.

var1 := [seq(seq(U[i, j], j = 0 .. JJ), i = 0 .. II), seq(seq(V[i, j], j = 0 .. JJ), i = 0 .. II), seq(seq(W[i, j], j = 0 .. JJ), i = 0 .. II), seq(seq(Phi[i, j], j = 0 .. JJ), i = 0 .. II), seq(seq(Xi[i, j], j = 0 .. JJ), i = 0 .. II)]

How to unassign them?