How to solve eigenvalue problem with eigenvalue apearing in the derivative terms

[Yang Yang]

Hello,

I’m trying to use Chebfun to solve an eigenvalue problem in atmosphere science. The ODE and its two boundary conditions are as follows:

k is the wavenumber as a known parameter, c is the unknown parameter (complex; considered as the eigenvalue). Given a certain k (say k=0.1), there would be many eigenvalues c and corresponding eigenfunctions. Among those, the most unstable mode is needed, which means finding the eigenvalue c with largest imaginary value.

Since c is the desired eigenvalue, I changed the above equation a little:

 

In this case, there are unknow parameter c appearing on the left-hand side of the equation. I would like to know how to use Chebfun to tackle problem like this?

 

The following is my “unfinished” codes....

 

k = 0.1; 

z = chebfun('z',[0 1]);

N = chebop( @(z,w) (1-k^2*(c+z)^2)*diff(w,2)-2/(c+z)*diff(w,1)-(2*k^2-c)*w, [0 1] ); 

N.lbc = 0;              % fixed end

N.rbc = 0;

[V,D] = eigs(N,10);

 

 

I would be very grateful if anyone can help me with this.

 

 Best regards,

Yang

[Nick Hale]

Hi Yang

What you have here is a nonlinear eigenvalue problem, which are a

bit more complicated than standard eigenvalue problems A*w = c*w.

Fortunately, by multiplying through by (c+z), you can turn your particular 

problem into a polynomial eigenvalue problem of the form 

         [c^3*A_3 + c^2*A_2*w + c*A_1 + A_0]*w = 0,

which simplifies matters a lot.

There's a routine polyeigs in MATLAB for solving matrix polynomial 

eigenvalue problems, but unfortunately the latest version of Chebfun

has not implemented this for differential equations. (This functionality

was in Chebfun Version 4, but it seems it was not introduced in the 

new version -- I might get a student to look at doing this if I get time.)

You could still use Chebfun to construct the discretised linear operators

(with boundary conditions) and then call the MATLAB polyeigs function. 

Alternatively you could try and use the old version of Chebfun.

Best regards

Nick

--
You received this message because you are subscribed to the Google Groups "chebfun-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to chebfun-user...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/chebfun-users/a19e1bf6-b79e-483d-b072-6de640489776n%40googlegroups.com.

[Yang Yang]

Hi Nick,

Thanks for your helpful suggestion. I've solved the problem.

Best,

Yang

[Nick Hale]

This ended up being easier than expected, so I spent a couple of hours re-implementing polyeigs.

There's a pull request here: https://github.com/chebfun/chebfun/tree/feature-polyeigs

The code below solves your problem (I think).

Best regards

Nick

PS. It would be nice to be able to specify problems such as this (and even standard/

generalised eigenvalue problems) as a single chebop with the eigenvalue included 

as a parameter in the function handle, and then use the automatic differentiation and 

parsing technologies we already have implemented in order to figure out  A, B, C, and D.

k = .1;

A = chebop(@(z,w) z*(1-k^2*z^2)*diff(w,2)- 2*diff(w) - 2*k^2*z*w);

A.bc = 'dirichlet';

B = chebop(@(z,w) (1-3*k^2*z^2)*diff(w,2)- 2*k^2*z*w);

C = chebop(@(z,w) -3*k^2*z*diff(w, 2));

D = chebop(@(z,w) -k^2*diff(w, 2));

[w, c] = polyeigs(A, B, C, D, 1, 'LI')

z = chebfun('z');

res = A(z,w) + c*B(z,w) + c^2*C(z,w) + c^3*D(z,w);

norm(res)

figure(1)

plot(z, real(w), z, imag(w), 'linewidth', 2)

xlabel('z'), legend('real w', 'imag w')

title(['c = ', num2str(c)])

To view this discussion on the web, visit https://groups.google.com/d/msgid/chebfun-users/78ad17ea-00f1-44df-a299-726699bc7aaen%40googlegroups.com.