Inaccurate eigenvalue/eigenfunctions approximation for 4th order differential equation

[Chi Zhang]

Hi everyone,

I have been working on a eigen problem for a 4th order differential equation and I found the results returned by eigs are way from satisfaction. The question I tried to solve is

L = chebop(@(x,u) diff(u,4), [0, 1]);

L.lbc = @(u) [diff(u,2);diff(u,3)];

L.rbc = @(u) [diff(u,2);diff(u,3)];

[V, D] = eigs(L);

norm(L*V - V*D, inf) % which returns 0.4833, while it should be (near) zero.

sum(V(:,1)*V(:,2)) % which returns 0.2506, while it should be (near) zero.

I tried the second-order version where eigs gave me a pretty good approximation

L = chebop(@(x,u) -diff(u,2), [0, 1]);

L.lbc = @(u) [diff(u,1)];

L.rbc = @(u) [diff(u,1)];

[V, D] = eigs(L);

norm(L*V - V*D, inf) % which returns 6.1555e-08.

sum(V(:,1)*V(:,2)) % which returns -1.2991e-09.

A general version of this problem can be found in, e.g. Equation (2.11) in Zuofeng Shang. Guang Cheng. "Local and global asymptotic inference in smoothing spline models." Ann. Statist. 41 (5) 2608 - 2638, October 2013. https://doi.org/10.1214/13-AOS1164 

I wonder if someone can help me with that. An example code has attached. Thanks in advance for your help.

Cheers,

Chi

[Nick Hale]

Hi Chi

Thanks for the interesting observation. Unfortunately I do not have a solution for you.

The same problem arises when discretising the problem manually (code below).

The first two eigenfunctions (spanning the null space) are far from orthogonal.  

Perhaps someone else has some insight?

Nick

n = 33;

x = chebpts(n, [0 1]);

t = chebpts(n-4, [0 1], 1);

P = barymat(t, x);

D2 = diffmat(n, 2, [0 1]);

D3 = diffmat(n, 3, [0 1]);

D4 = diffmat(n, 4, [0 1]);

% Un-comment for row-rep rather than rectangular projection:

% P = eye(n,n); P([1:2, end-1:end],:) = []; % Row replacement

A = [D2(1,:) ; D3(1,:) ; P*D4 ; D2(end,:) ; D3(end,:)];

B = [zeros(2, n) ; P ; zeros(2, n)];

[V, D] = eig(A, B);

d = diag(D);

[~, idx] = sort(abs(d), 'ascend');

d = d(idx); V = V(:,idx);

d(1:6)

%ans =
%   1.0e+04 *
%  -0.000000000750923
%   0.000000000819956
%   0.050056400486424
%   0.380353718274010
%   1.461763030034135
%   3.994379917511742

chebfun(V(:,1:4))'*chebfun(V(:,1:4))*4

%ans =
%   1.499263190538322  -0.111921196347822   0.000000071648358  -0.000000016095639
%  -0.111921196347822   1.000000366556549  -0.000000032035889   0.000000001932542
%   0.000000071648358  -0.000000032035889   1.000000020851751   0.000000060353897
%  -0.000000016095639   0.000000001932542   0.000000060353897   0.999999993844841

[Tobin Driscoll]

I am inclined to put the blame on roundoff/conditioning. A fourth-order derivative, with BCs on derivatives 2 and 3, is begging for trouble, as the norms of the matrices grow rapidly. 

Nick's code is strongly hinting that the first two modes are linear functions (zero eigenvalue); the rest of the polynomial coefficients are 5 orders of magnitude smaller. At a repeated eigenvalue, matlab will return an eigenbasis that is orthogonal in the discrete sense but not in a true Chebyshev inner product, so I'm not concerned about the upper 2x2 block of V'*V.

You might be able to get a few digits of accuracy for the first few modes, but it's unlikely you will go safely past n=100 in size in double precision.

It would probably be better to formulate the problem for u'''' rather than for u. See (25)-(26) in http://www.chebfun.org/publications/driscoll2010.pdf, where the Orr-Sommerfeld eigenvalues are computed. That process was never automated in Chebfun.

-- Toby

--
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/0a510e65-6174-46b0-9a3a-7ac9498ac1b4n%40googlegroups.com.

[nakatsu...@gmail.com]

This is interesting and I wanted to try an algorithm I developed with Behnam Hashemi in this paper
https://epubs.siam.org/doi/pdf/10.1137/21M1445934

I'm attaching the code LSeig.m, which is needed to run the snippet below.

clear, clc
dom = [0, 1];
x   = chebfun('x', dom);
op  = @(x,u) diff(u,4);

n = 50; % degree of polynomial basis
T = chebpoly(0:n-1,dom);
dT2 = diff(T,2);
dT3 = diff(T,3);
A = op(x,T); B = T;
u0 = dT2(dom(1),:); u1 = dT3(dom(1),:);
u2 = dT2(dom(2),:); u3 = dT3(dom(2),:);
BC = [u0;u1;u2;u3]; BClam = zeros(size(BC));

tol = 1e-10; % accept solns with <tol residual
    [V,D,res] = LSeig(A,B, BC, BClam, tol,'IM');
    %[V,D,res] = LSeig(A,B, BC, BClam, tol,'qr'); % this works too
    for ii = 1:size(V,2), V(:,ii) = V(:,ii)/norm(B*V(:,ii)); end % rescale soln to unit norm

RES = A*V-B*V*D; % (some of) this should be small
disp('eigvals')
disp(diag(D).')
disp('residuals')
for ii = 1:size(RES,2), resnorms(ii) = norm(RES(:,ii)); end
disp(resnorms)
disp('boundary conditions') % some may not satisfy them
disp(vecnorm([u0*V(:,:);u1*V(:,:);u2*V(:,:);u3*V(:,:)]))

disp('ensure orthogonality of eigenvectors')
disp((T*V(:,1))'*(T*V(:,2)))  

Output:
eigvals
     0     0

residuals
     0     0

boundary conditions
     0     0

ensure orthogonality of eigenvectors
   1.5450e-18

It seems to work alright here, returning the constant and linear functions as solutions. The algorithm doesn't find any other solutions though.

Yuji

[nakatsu...@gmail.com]

Dear all, 

Behnam has explored this further, and told me that we can get the other solutions by using LSeigbc.m (attached), a variant of the method from yesterday that imposes the boundary conditions exactly (to numerical precision). Below is a MATLAB snippet. I'm also attaching a plot of the first three eigenfunctions (one of which we hadn't seen; the method computes many more).

A related class of methods are also developed in a paper with Behnam and Nick Trefethen here:
https://doi.org/10.1007/s10444-022-09994-8

Yuji

clc, clear

dom = [0, 1];

x = chebfun('x', dom);

op = @(x,u) diff(u,4);

n = 50;

tol = 1e-10;

T = chebpoly(0:n-1,dom);

dT2 = diff(T,2);

dT3 = diff(T,3);

A = op(x,T); B = T;

u0 = dT2(dom(1),:); u1 = dT3(dom(1),:);

u2 = dT2(dom(2),:); u3 = dT3(dom(2),:);

BC = [u0;u1;u2;u3]; BClam = zeros(size(BC));

[V,D,resnorms] = LSeigbc(A,B,BC,BClam,tol);

%[V,D] = LSeig(A,B,BC,BClam,tol);

ei = diag(D); [ei,ix] = sort(ei,'ascend'); V = V(:,ix); D = D(ix,ix);

disp('eigvals relative residuals')

fprintf('%15.10f \t %3.2e\n',[diag(D) resnorms].')

fprintf('\n boundary conditions\n')

fprintf('%3.2e \t %3.2e\n',[u0*V; u1*V].')

fprintf('\n orthogonality of eigenfunctions\n')

kk = 5; BB = T*V(:,1:kk); disp(norm(BB'*BB-diag(diag(BB'*BB))))

close all, eigfuns = T*V; plot(eigfuns(:,1:3)), title('three eigenfunctions'), grid on

Output:

eigvals                 relative residuals
  -0.0000000000  0.00e+00
  -0.0000000000  0.00e+00
508481.5432671806  2.08e-11
793403.1345424439  2.62e-11
1184013.5895825236  2.19e-11
1703691.0902422448  4.67e-11
2378151.6365684960  1.85e-11
3235449.0467638792  2.87e-11
4305974.9572528321  8.48e-11

 boundary conditions
0.00e+00  0.00e+00
-1.01e-09  1.95e-09
1.02e-09  -8.73e-11
-1.19e-09  3.78e-09
1.63e-09  0.00e+00
0.00e+00  7.05e-07
-2.84e-06  -2.51e-07
-2.64e-07  2.98e-07
-3.61e-06  -1.04e-06

 orthogonality of eigenfunctions
   2.9203e-11

[Chi Zhang]

Thank you Yuji! I personally try to solve it manually and to see if I can find closed forms for eigenvalues. The resulting equation I obtained for eigenvalues is (exp(lambda^(1/4) + exp(-lambda^(1/4))*cos(lambda^(1/4)) - 2 = 0. The first two non-zero positive root is about 500.5467 and 3803.537 which matches what Nick provided above. However, the associated eigenvalue in your previous thread is about 508481 and 793403. I wonder if there is something wrong with the computation.

You received this message because you are subscribed to a topic in the Google Groups "chebfun-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/chebfun-users/m0gTtq1Gug0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to chebfun-user...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/chebfun-users/56d8b91f-2ba6-4ef1-b1f0-bc92f18f1956n%40googlegroups.com.

[Behnam Hashemi]

Hi Chi,

There is nothing wrong in the computation. The approximate eigenvalues 508481 and 793403 that Yuji reported in his post should not be considered as the first two nonzero eigenvalues! Those are the first two nonzero eigenvalues whose relative residuals as computed with LSeigbc with n = 50 basis functions are less than the tolerance 10^(-10) specified in that case. Exploring different number of basis functions and different tolerances, can give you many more eigenvalues. For example, with tol = 10^(-5) LSeigbc gives the following 21 eigenvalues (n=50 as before):

eigvals                                   relative residuals
 -0.0000000000 0.00e+00
 -0.0000000000  0.00e+00

 500.5639017553  5.20e-06
3803.5370819586  1.08e-07
14617.6301303037  1.05e-09
39943.7989855812  2.83e-07
89135.4076555530  2.97e-11
173881.3154732978  3.54e-09
308208.4521062748  3.90e-11
508481.5432682041  7.14e-11

793403.1345390277  2.18e-11
1184013.5895911232  2.25e-10
1703691.0902372699  1.77e-11
2378151.6365731005  5.88e-10
3235449.0467399466  7.10e-08
4305974.9572284315  3.65e-09
5622458.8225279972  3.68e-11
7219967.9155088663  7.91e-11
9135907.3270380050  2.39e-10
11410019.9663404338  2.35e-08
14084386.5610736236  1.03e-09

The 3rd and the 4th eigenvalues in this list seem to be the first two nonzero eigenvalues. In addition, the eigenvalues 508481 and 793403 in Yuji's post appear here again; see rows 10 and 11 in the list.

By the way, the equation written in your post for lambda is not clear to me as it contains seven open parentheses but five closed ones. 

Good luck,

Behnam

[Chi Zhang]

Thank you everyone! Very helpful! I indeed obtained a very accurate eigen-systems. By the way, the equation I meant is (exp(x) + exp(-x))*cos(x) = 2 where x = lambda ^(1/4), and lambda is the eigenvalue. That is, all these non-negative roots for this equation are the eigenvalues for this particular problem.

To view this discussion on the web, visit https://groups.google.com/d/msgid/chebfun-users/29e9e36f-e54a-442a-b381-ce47663f4e5en%40googlegroups.com.