Viscous_Upper.m

%--------------------------------------------------------------------------
% Bounding the minimum wave speed of the viscous FKKP-Burgers equation: 
%
%    T_t - T_xx + (u*T)_x = T*(1-T),
%   u_t - nu*u_xx + u*u_x = rho*T*(1-T),    rho,nu > 0
%
% This code produces an upper bound on the (conjectured) minimum wave speed 
% and is associated to the paper "The speed of traveling waves in a 
% FKPP-Burgers system" by Jason J. Bramburger and Christopher Henderson 
% (2020). This script is used to generate the data in Table 1.
%--------------------------------------------------------------------------

% Access YALMIP and mosek directories
addpath(genpath('YALMIP-master')) 
addpath(genpath('mosek')) 

% Clean workspace
clear all
close all
clc

format long

% Differential equation Parameters
nu = 1;
rho = 10^6;

% Bounding method parameters
d = 10; %Degree of H(T,U,V)
lambda = 1;
    
%Bisection Method
cleft = 100;
cright = 150;

while abs(cright - cleft) >= 1e-4

    cmid = 0.5*(cleft + cright);
    flag = upperbnd(cmid,rho,nu,d,lambda);

    if flag == 0
       cright = cmid; 
    else
       cleft = cmid;
    end
end

%Print Results
fprintf('The upper bound on the minimum speed for rho = %d and nu = %d is %f found using degree %d polynomials.\n',rho,cmid,d)

%%
function flag = upperbnd(c,rho,nu,d,lambda)

% SDP variables
sdpvar T U V

% Source fixed point value
u0 = c + rho - sqrt(c^2 + rho^2);
v0 = c - u0;

% Epsilon value
eps = 1e-4;

% Auxiliary function
[H, cH] = polynomial([T U V],d);

% S Procedure
d2 = d + 4;
[s1, c1] = polynomial([T U V], d2);
[s2, c2] = polynomial([T U V], d2);
[s3, c3] = polynomial([T U V], d2);
[s4, c4] = polynomial([T U], d2);
[s5, c5] = polynomial([T U], d2);

% Derivatives
dHdT = jacobian(H,T);
dHdU = jacobian(H,U);
dHdV = jacobian(H,V);

% Function replacements
HV0 = replace(H,V,0); %H(T,U,0)

%Constraints
cons = [];
cons = [cons, replace(H,[T U V], [0 0 0]) == 0]; 
cons = [cons, replace(H,[T U V], [1 1 1]) <= -eps];
cons = [cons, sos(-lambda*(dHdT*(-c*T + U*u0*T + v0*V) + dHdU*((-c*U*u0 + 0.5*(u0^2)*U^2 + rho*v0*V)/(nu*u0)) + dHdV*(T*(T-1)/v0)) - H - T*(1-T)*s1 - U*(1-U)*s2 - V*(1-V)*s3)];
cons = [cons, sos(HV0 - eps*T*(1-T) - eps*U*(1-U) - T*(1-T)*s4 - U*(1-U)*s5)];
cons = [cons, sos(s1), sos(s2), sos(s3), sos(s4), sos(s5)];

%SOS Solver
ops = sdpsettings('solver','mosek','verbose',0,'cachesolvers',1);
sol = solvesos(cons,[],ops,[cH; c1; c2; c3; c4; c5])

%Return whether solvesos failed or succeeded
flag = sol.problem;

end