Generalized Orr Sommerfield error returned in Eigs Function

Question

IRENE 2024-4-14，21:02

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2106786-generalized-orr-sommerfield-error-returned-in-eigs-function

编辑： Torsten 2024-4-17，17:53

orr_sommerfield_chebfun_solution.m

在 MATLAB Online 中打开

Good evening!

I am trying to find the eigenvalues for a generalized Orr-Sommerfeld equation. The code is the file attached above.

I based my code from the one readable here.

However, it returns error. I ran line by line until finding the problem line, it is this one :

%calcolo autovalori
e0 = eigs(A0,B,50,'LR');
e1 = eigs(A1,B,50,'LR');

Which returns the following errors:

% Error using lambertw
%  W0 = @(x) lambertw(0,x);
%                       ↑
% Invalid argument at position 2. Value must be numeric.
% 
% Error in orr_sommerfield_chebfun_solution>@(x)lambertw(0,x) (line 18)
% W0 = @(x) lambertw(0,x);
% 
% Error in orr_sommerfield_chebfun_solution>@(x)exp(W0(b.*x)).*(W0(b.*x)+1) (line 37)
% f0 = @(x) exp(W0(b.*x)) .* (W0(b.*x)+1);
% 
% Error in orr_sommerfield_chebfun_solution>@(x,u)((f0(x).*(l.^4.*diff(u,4)-2.*(ad.*l).^2.*diff(u,2)+ad.^4.*u))+((2.*df0(x).*l.^3.*diff(u,3))+((2.*ad.^2.*(2.*f0(x)-g0(x))+ddf0(x)).*l.^2.*diff(u,2))+((2.*ad.^2.*(df0(x)-dg0(x))).*l.*diff(u,1))+(ad.^2.*ddf0(x).*u)))./Re-1i.*ad.*(U0(x).*(l.^2.*diff(u,2)-ad.^2.*u)+ddU0(x).*u) (line 59)
% A0.op = @(x,u) (  (f0(x).*(l.^4.*diff(u,4)-2.*(ad.*l).^2.*diff(u,2)+ad.^4.*u)) + ...
% 
% Error in chebop/feval (line 190)
%     out = N.op(x, u);
% 
% Error in chebop/linearize (line 150)
%         Nu = feval(N, x, u{:});
% 
% Error in chebop/eigs (line 59)
% [L, ~, isLinear, u0] = linearize(N, N.init, [], linCheck);
% 
% Error in orr_sommerfield_chebfun_solution (line 80)
% e0 = eigs(A0,B,50,'LR');

This question was similar to my situation, in the sense that it had problems with the eigs function. I suspect that my problem comes from all the @ functions I have defined, but I can't determine what I'm doing wrong.

Does anyone has an idea for the origin of these errors? Or has a suggestion to try and find a different approach to find the eigenvalues.

Thanks!!

close all
clc
Re = 5000;                  % Reynolds number
alph = 1;                   % longitudinal Fourier parameter
G = 5.5e3;                  % costante relativa al materiale e a condiz geometriche
Bn = G./Re;                 % Bingham number 
p = Bn+1 ;                  % adimensional pressure variation
k = 62.3;                   % de kee constant
s = Bn./p;                  % yield surface
del = 1e-5;                 % delta
b = (k.*p.*(s-1))./2;       %costante ottenuta dalla trasformazione
% funzioni di Lambert
W0 = @(x) lambertw(0,x);
W1 = @(x) lambertw(-1,x);
% costanti
C0 = ((s-1)./k).*( W0((k.*p./2).*(s-1)) -1 + 1./( W0((k.*p./2).*(s-1)) ) );
C1 = ((s-1)./k).*( W1((k.*p./2).*(s-1)) -1 + 1./( W1((k.*p./2).*(s-1)) ) );
% Campi di velocità e loro derivate già trasformati
U0 = @(x) -((s-1)./(k.*b)) .* exp(W0(b.*x)) .* ( (W0(b.*x)).^2 - W0(b.*x) + 1 ) + C0 ;
U1 = @(x) -((s-1)./(k.*b)) .* exp(W1(b.*x)) .* ( (W1(b.*x)).^2 - W1(b.*x) + 1 ) + C1 ;
dU0 = @(x) (1./k) .* W0(b.*x);
dU1 = @(x) (1./k) .* W1(b.*x);
ddU0 = @(x) -(p./2) .* ( 1./(exp(W0(b.*x)).*(W0(b.*x))) );
ddU1 = @(x) -(p./2) .* ( 1./(exp(W1(b.*x)).*(W1(b.*x))) );
%funzioni per OSE
f0 = @(x) exp(W0(b.*x)) .* (W0(b.*x)+1);
f1 = @(x) exp(W1(b.*x)) .* (W1(b.*x)+1);
df0 = @(x) b./(1-s) .* (W0(b.*x)+2)./(W0(b.*x)+1);
df1 = @(x) b./(1-s) .* (W1(b.*x)+2)./(W1(b.*x)+1);
ddf0 = @(x) (b./(1-s)).^2 .* 3./( exp(W0(b.*x)).*(W0(b.*x)+1).^3 );
ddf1 = @(x) (b./(1-s)).^2 .* 3./( exp(W1(b.*x)).*(W1(b.*x)+1).^3 );
g0 = @(x) 2.*exp(W0(b.*x)) - ((2.*Bn.*k)./W0(b.*x));
g1 = @(x) 2.*exp(W1(b.*x)) - ((2.*Bn.*k)./W1(b.*x));
dg0 = @(x) b./(1-s) .* (  2./(W0(b.*x)+1) + (2.*Bn.*k)./( exp(W0(b.*x)).*(W0(b.*x)+1).*(W0(b.*x))^2 )  );
dg1 = @(x) b./(1-s) .* (  2./(W1(b.*x)+1) + (2.*Bn.*k)./( exp(W1(b.*x)).*(W1(b.*x)+1).*(W1(b.*x))^2 )  );
% operatore A
ad = alph.*del;
l = (1-s);
A0 = chebop(0,1);
A1 = chebop(0,1);
A0.op = @(x,u) (  (f0(x).*(l.^4.*diff(u,4)-2.*(ad.*l).^2.*diff(u,2)+ad.^4.*u)) + ...
    ( (2.*df0(x).*l.^3.*diff(u,3)) + ((2.*ad.^2.*(2.*f0(x)-g0(x))+ddf0(x)).*l.^2.*diff(u,2)) + ((2.*ad.^2.*(df0(x)-dg0(x))).*l.*diff(u,1)) + (ad.^2.*ddf0(x).*u) ) ...
    )./Re - 1i.*ad.*(U0(x).*(l.^2.*diff(u,2)-ad.^2.*u)+ddU0(x).*u) ;
A1.op = @(x,u) (  (f1(x).*(l.^4.*diff(u,4)-2.*(ad.*l).^2.*diff(u,2)+ad.^4.*u)) + ...
    ( (2.*df1(x).*l.^3.*diff(u,3)) + ((2.*ad.^2.*(2.*f1(x)-g1(x))+ddf1(x)).*l.^2.*diff(u,2)) + ((2.*ad.^2.*(df1(x)-dg1(x))).*l.*diff(u,1)) + (ad.^2.*ddf1(x).*u) ) ...
    )./Re - 1i.*ad.*(U1(x).*(l.^2.*diff(u,2)-ad.^2.*u)+ddU1(x).*u) ;
% operatore B, non dipende da U0 o U1!!
B = chebop(0,1);
B.op = @(x,u) (l.^2.*diff(u,2) - ad.^2.*u);
%boundary conditions [v = 0; v' = 0]
A0.lbc = [0; 0];
A0.rbc = [0; 0];
A1.lbc = [0; 0];
A1.rbc = [0; 0];
%calcolo autovalori
e0 = eigs(A0,B,50,'LR');
e1 = eigs(A1,B,50,'LR');
MS0 = 'markersize';
maxe0 = max(real(e0));
MS1 = 'markersize';
maxe1 = max(real(e1));
%plot dei grafici
figure(1);
plot(e0,'.r',MS0,14), grid on, axis([-.9 .1 -1 0]), axis square
title(sprintf('Re = %8.2f   \\lambda_r = %7.5f',Re,maxe0))
figure(2);
plot(e1,'.r',MS1,14), grid on, axis([-.9 .1 -1 0]), axis square
title(sprintf('Re = %8.2f   \\lambda_r = %7.5f',Re,maxe1))
hold off

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-4-15，10:28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2106786-generalized-orr-sommerfield-error-returned-in-eigs-function#answer_1441371

编辑：Torsten 2024-4-15，10:32

在 MATLAB Online 中打开

I suggest you comment out the command

W0 = @(x) lambertw(0,x);

like

%W0 = @(x) lambertw(0,x);

include a function for W0 as

function value = W0(x)
  x
  class(x)
  value = lambertw(x,0)
end

and check what goes wrong with your input argument x in W0.

Maybe same for W1.

It's quite an effort for us to run your code because not everyone has "chebfun" installed, I guess.

2 个评论
显示无隐藏无

IRENE 2024-4-17，17:18

在 MATLAB Online 中打开

Hello, I tried what you suggested, now I get this error messages. I don't understand what the function can do to help me, could you explain more about it?

Thanks for your help!

ans =
    'double'
% Error using lambertw
% First argument must be an integer.
% 
% Error in orr_sommerfield_chebfun_solution>W0 (line 103)
%   value = lambertw(x,0);
% 
% Error in orr_sommerfield_chebfun_solution (line 23)
% C0 = ((s-1)./k).*( W0((k.*p./2).*(s-1)) -1 + 1./( W0((k.*p./2).*(s-1)) ) );

Torsten 2024-4-17，17:51

编辑：Torsten 2024-4-17，17:53

在 MATLAB Online 中打开

I tried what you suggested, now I get this error messages.

I chose the wrong order of the arguments for W0. Use

function value = W0(x)
  x
  class(x)
  value = lambertw(0,x)
end

instead of

function value = W0(x)
  x
  class(x)
  value = lambertw(x,0)
end

I don't understand what the function can do to help me, could you explain more about it?

Form the error message

% Error using lambertw
%  W0 = @(x) lambertw(0,x);
%                       ↑
% Invalid argument at position 2. Value must be numeric.

I would conclude that lambertw is called with an input for x that is not supported by the function.

In order to find out the class and the value of this input, you will need to debug your code. Because function handles cannot be debugged, I suggested to use a function for W0 instead of a function handle and to output for which x-values the function is called and for which input for x the error finally occurs.

请先登录，再进行评论。