Generalized Orr Sommerfield error returned in Eigs Function
40 views (last 30 days)
Show older comments
Good evening!
I am trying to find the eigenvalues for a generalized Orr-Sommerfeld equation. The code is the file attached above.
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 Comments
Answers (1)
Torsten
on 15 Apr 2024 at 10:28
Edited: Torsten
on 15 Apr 2024 at 10:32
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 Comments
Torsten
on 17 Apr 2024 at 17:51
Edited: Torsten
on 17 Apr 2024 at 17:53
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.
See Also
Categories
Find more on Logical in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!