COMPUTE_WHISKER

function [W_0, R_0] = compute_whisker(obj, order)

This function computes the invariant manifold in the autonomous system limit either a tensor- based or a multi-index based algorithm may be chosen for the computation. This function sets up the computation, the order $i$ cohomological equation for obtaining the parametrisation coefficients and the coefficients for the reduced dynamics are computed by calling cohomological_solution.

Leading-order terms

Lambda_E = obj.E.spectrum;
W_01     = obj.E.basis;

switch obj.Options.notation

    case 'tensor'
        % Initialize
        W_0 = cell(1,order);
        R_0 = cell(1,order);

        W_0{1} = sptensor(W_01);
        R_0{1} = sptensor(diag(Lambda_E));

        % *Outer resonance case* : issue warning
        if obj.resonance.outer.occurs
            prompt = 'Due to (near) outer resonance, the exisitence of the manifold is questionable and the underlying computation may suffer.';
            disp(prompt)
            disp('Attempting manifold computation')
        end

        for j = 2:order
            %recursively approximate at j-th order
            startOrderj = tic;
            [W_0{j},R_0{j},~] = cohomological_solution(obj,j,W_0,R_0,[]);
            obj.solInfo.timeEstimate(j) = toc(startOrderj);
            disp(['Manifold computation time at order ' num2str(j) ' = ' datestr(datenum(0,0,0,0,0,obj.solInfo.timeEstimate(j)),'HH:MM:SS')])
            fprintf('Estimated memory usage at order % 2i = %05.2E MB\n', j, obj.solInfo.memoryEstimate(j))
        end

        for j = 1:order
            W_0{j} = tensor_to_multi_index(W_0{j});
            R_0{j} = tensor_to_multi_index(R_0{j});
        end
        W_0 = [W_0{:}];
        R_0 = [R_0{:}];
    case 'multiindex'

        W_0 = repmat(struct('coeffs',[]),1,order);
        R_0 = repmat(struct('coeffs',[]),1,order);

        % Set up system for Multi-index notation version
        check_COMPtype(obj); % determine whether system is real or complex
        [W_0(1),R_0(1),multi_input] = coeffs_setup(obj,order);

        % *Outer resonance case* : issue warning
        if obj.resonance.outer.occurs
            prompt = 'Due to (near) outer resonance, the exisitence of the manifold is questionable and the underlying computation may suffer.';
            disp(prompt)
            disp('Attempting manifold computation')
        end

        for j = 2:order
            %recursively approximate at j-th order
            startOrderj = tic;
            [W_0(j).coeffs,R_0(j).coeffs,multi_input] = cohomological_solution(obj,j,W_0,R_0,multi_input);
            obj.solInfo.timeEstimate(j) = toc(startOrderj);
            disp(['Manifold computation time at order ' num2str(j) ' = ' datestr(datenum(0,0,0,0,0,obj.solInfo.timeEstimate(j)),'HH:MM:SS')])
            fprintf('Estimated memory usage at order % 2i = %05.2E MB\n', j, obj.solInfo.memoryEstimate(j))
        end
        switch obj.System.Options.DStype
            case 'real'
                [W_0,R_0] = coeffs_conj2lex(multi_input,order,W_0,R_0);
            case 'complex'
                [W_0] = coeffs_lex2revlex(W_0,'TaylorCoeff');
                [R_0] = coeffs_lex2revlex(R_0,'TaylorCoeff');
        end

Add multi-indices field to the coefficients field at every order

        [W_0,R_0] = coeffs_output(W_0,R_0,order);

end