Forward NFT for KdV

examples/fnft_kdvv_example.c

@@ -15,6 +15,7 @@
  *
  * Contributors:
  * Sander Wahls (TU Delft) 2018.
+ * Peter J. Prins (TU Delft) 2021.
  */
 
 // This example demonstrates the use of the function fnft_kdvv, which 
@@ -59,15 +60,28 @@ int main()
     // Buffer for result: Samples of the continuous spectrum
     FNFT_COMPLEX contspec[M];
 
+    // Maximum number of bound states we expect = size of the two arrays below
+    FNFT_UINT K = D;
+
+    // Buffer for result: Bound states
+    FNFT_COMPLEX bound_states[K];
+
+    // Buffer for result: Norming constants
+    FNFT_COMPLEX normconsts[K];
+
     // Default options
     fnft_kdvv_opts_t opts = fnft_kdvv_default_opts();
 
+    // Uncomment the next line to compute residues instead of norming constants
+    //opts.discspec_type = fnft_kdvv_dstype_RESIDUES;
+
     // See the header file fnft_kdvv.h for the options
 
     /** Step 3: Call fnft_nsev and check for errors **/
 
-    int ret_code = fnft_kdvv(D, q, T, M, contspec, XI, NULL, NULL, NULL,
-        &opts);
+    int ret_code = fnft_kdvv(D, q, T, M, contspec, XI, &K, bound_states,
+                             normconsts, &opts);
+
     if (ret_code != FNFT_SUCCESS) {
         printf("An error occured!\n");
         return EXIT_FAILURE;
@@ -88,5 +102,15 @@ int main()
         );
     }
 
+    printf("Discrete spectrum:\n");
+    for (FNFT_UINT i=0; i<K; i++) {
+        printf("  bound state at %g + %gI with norming constant %g + %gI\n",
+            (double)FNFT_CREAL(bound_states[i]),
+            (double)FNFT_CIMAG(bound_states[i]),
+            (double)FNFT_CREAL(normconsts[i]),
+            (double)FNFT_CIMAG(normconsts[i])
+        );
+    }
+
     return EXIT_SUCCESS;
 }

examples/mex_fnft_kdvv_example.m

@@ -14,10 +14,11 @@
 %
 % Contributors:
 % Sander Wahls (TU Delft) 2017-2018.
+% Peter J. Prins (TU Delft) 2021.
 
 % This examples demonstrates how the nonlinear Fourier transform with
 % respect to the Korteweg-de Vries equation with vanishing boundary
-% conditions can be computed using mex_fnft_nsev. The signal is a squared
+% conditions can be computed using mex_fnft_kdvv. The signal is a squared
 % sech pulse.
 
 clear all;
@@ -26,35 +27,48 @@
 %%% Setup parameters %%%
 
 
-D = 2^8;        % location of the 1st and last sample in the time domain
-T = [-10, 10];  % number of samples
-XI = [-5, 5];   % location of the 1st and last sample in the xi-domain
+D = 2^8;            % location of the 1st and last sample in the time domain
+T = [-10, 10];      % number of samples
+XI = [1e-10, 6];   % location of the 1st and last sample in the xi-domain
 
 %%% Setup the signal %%%
 
 ep_t = (T(2) - T(1)) / (D - 1);
 t = T(1):ep_t:T(2);
-q = complex(1.2*sech(t).^2);
+q = 16*sech(t).^2;
 
 %%% Compute the nonlinear Fourier transform %%%
 
-contspec = mex_fnft_kdvv(q, T, XI);
+[contspec, bound_states, norming_constants] = mex_fnft_kdvv(q, T, XI);
+% mex_fnft_kdvv has many options => run "help mex_fnft_kdvv" to learn more
 
 %%% Plot the results %%%
 
 ep_xi = (XI(2) - XI(1)) / (D - 1);
 xi = XI(1):ep_xi:XI(2);
 
 figure;
-plot(t, real(q), t, imag(q));
+plot(t, q);
 title('Time-domain');
 xlabel('t');
 ylabel('q(t)');
-legend('Real part', 'Imaginary part');
+legend('q(t)');
 
 figure;
-plot(xi, real(contspec), xi, imag(contspec));
-title('Continuous spectrum (numerically)');
+hmag=subplot(2,1,1);
+semilogy(xi, abs(contspec));
+title('Continuous spectrum');
 xlabel('\xi');
-ylabel('r(\xi)');
-legend('Real part', 'Imaginary part');
+ylabel('|r(\xi)|');
+hang=subplot(2,1,2);
+plot(xi, angle(contspec));
+xlabel('\xi');
+ylabel('\angle r(\xi)');
+linkaxes([hmag,hang],'x');
+
+figure;
+stem(imag(bound_states),real(norming_constants), 'x');
+title('Bound states and norming constants');
+xlabel('\lambda_k/j');
+ylabel('b(\lambda_k)');
+axis([0,4,-2,2]);

examples/mex_fnft_kdvv_example_2.m

@@ -0,0 +1,120 @@
+% This file is part of FNFT.
+%
+% FNFT is free software; you can redistribute it and/or
+% modify it under the terms of the version 2 of the GNU General
+% Public License as published by the Free Software Foundation.
+%
+% FNFT is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+% Contributors:
+% Sander Wahls (TU Delft) 2017-2018, 2020.
+% Shinivas Chimmalgi (TU Delft) 2019.
+% Peter J. Prins (TU Delft) 2021.
+
+% This example compares the accuracy-execution tradeoff of the various algorithms for
+% the Korteweg-de Vries equation with vanishing boundaries.
+
+clear;
+close all;
+clc;
+
+%%% Setup parameters %%%
+
+T = [0,16];     % location of the 1st and last sample in the time domain
+XI = [0.5,23];  % location of the 1st and last sample in the xi-domain
+
+%%% Setup the signal %%%
+
+A = sym(15);
+d = sym(0.5);
+exp_t0 = sym(3000);
+q_fun = @(t) A*sech((t-log(exp_t0))/d).^2;
+
+%%% Analytic expressions for the nonlinear Fourier transform %%%
+
+delta = sqrt(sym(A)*d^2+1/4);
+at = @(xi) 1/2 - 1i*xi*d + delta;
+bt = @(xi) 1/2 - 1i*xi*d - delta;
+ct = @(xi) 1 - 1i*xi*d;
+cgam = @(z) gamma(sym(z));
+
+R_fun = @(xi) exp_t0.^(-2i*xi) .* cgam(ct(xi)-at(xi)-bt(xi)) .* cgam(at(xi)) .* cgam(bt(xi)) ./ ...
+                       ( cgam(ct(xi)-at(xi)) .* cgam(ct(xi)-bt(xi)) .* cgam(at(xi)+bt(xi)-ct(xi)) );
+
+%%% Prepare variables to store errors and runtimes %%%
+
+error_FCF2_1 = [];
+error_FCF4_2 = [];
+error_FCF_RE2_1 = [];
+error_FCF_RE4_2 = [];
+time_FCF2_1 = [];
+time_FCF4_2 = [];
+time_FCF_RE2_1 = [];
+time_FCF_RE4_2 = [];
+
+%%% Iterate of number of samples, gather errors and runtimes for each %%%
+
+for D = 2.^(8:1:12)
+
+    fprintf('Running codes with D=%d...',D);
+    t = linspace(T(1),T(2),D);
+    q = double(q_fun(t)); % signal samples
+    R_exact = double(R_fun(linspace(XI(1),XI(2),D)));
+
+    % Compute the continuous part of the nonlinear Fourier transform
+    % numerically with different configurations, save correspoding errors
+    % and runtimes.
+
+    tic
+    R_comp = mex_fnft_kdvv(q, T, XI, 'discr_2split4B', 'skip_bs');
+    time_FCF2_1 = [time_FCF2_1,toc];
+    error_FCF2_1 = [error_FCF2_1, norm(R_exact-R_comp)/norm(R_exact)];
+
+    tic
+    R_comp = mex_fnft_kdvv(q, T, XI, 'discr_4split4B', 'skip_bs');
+    time_FCF4_2 = [time_FCF4_2,toc];
+    error_FCF4_2 = [error_FCF4_2, norm(R_exact-R_comp)/norm(R_exact)];
+
+    tic
+    R_comp = mex_fnft_kdvv(q, T, XI, 'discr_2split4B', 'RE', 'skip_bs');
+    time_FCF_RE2_1 = [time_FCF_RE2_1,toc];
+    error_FCF_RE2_1 = [error_FCF_RE2_1, norm(R_exact-R_comp)/norm(R_exact)];
+
+    tic
+    R_comp = mex_fnft_kdvv(q, T, XI, 'discr_4split4B', 'RE', 'skip_bs');
+    time_FCF_RE4_2 = [time_FCF_RE4_2,toc];
+    error_FCF_RE4_2 = [error_FCF_RE4_2, norm(R_exact-R_comp)/norm(R_exact)];
+
+    fprintf('Done.\n');
+end
+
+%%% Plot results %%%
+
+lw = 3;
+fs = 15;
+ms = 6;
+alw = 1;
+loglog(time_FCF2_1,error_FCF2_1,'-','linewidth',lw,'markersize',ms);                           
+hold on
+loglog(time_FCF4_2,error_FCF4_2,'--','linewidth',lw,'markersize',ms);
+ax = gca;
+ax.ColorOrderIndex = 1;
+loglog(time_FCF_RE2_1,error_FCF_RE2_1,':','linewidth',lw,'markersize',ms);
+loglog(time_FCF_RE4_2,error_FCF_RE4_2,'-.','linewidth',lw,'markersize',ms);
+axis tight
+
+set(gca, 'Units','normalized','FontUnits','points',...
+    'FontWeight','normal','FontSize',fs,'FontName','Times','LineWidth',alw)
+
+xlabel({'Execution Time (s)'},'Interpreter','latex','FontUnits','points','Fontsize',fs,'FontName','Times');
+ylabel({'Relative $L^2$-error in $\rho(\lambda)$'},'Interpreter','latex','FontUnits','points','Fontsize',fs,'FontName','Times')
+grid on
+grid minor
+le = legend("FCF$^{[2]}_1$","FCF$^{[4]}_2$","FCF\_RE$^{[2]}_1$","FCF\_RE$^{[4]}_2$");
+set(le,'interpreter','latex','FontUnits','points','Fontsize',fs-2,'FontName','Times','Location','SouthWest');