EXTRACT_FRC

Contents

function FRC = extract_FRC(obj, parName, parRange, ORDER)
        

Extract Forced Response Curve

This function extracts the forced response curve (FRC) for systems that may have internal resonances. The FRC computation is based on SSM computation. An appropriate SSM is constructed based on the resonant spectrum of system. The FRC is computed for the reduced system either by level-set technique or continuation method. The obtained FRC is finally mapped back to physical coordinates. By default, level-set method is used for 2D, underdamped SSMs. In any case, the continuation method is employed for higher dimensional SSMs.

FRC = EXTRACT_FRC(OBJ, PARNAME, PARRANGE, ORDER)

See also: FRC_LEVEL_SET, FRC_CONT_EP, FRC_CONT_PO

f1 = figure('Name','Norm');
if isnumeric(obj.FRCOptions.outdof)
    f2 = figure('Name',['Amplitude at DOFs ' num2str(obj.FRCOptions.outdof(:)')]);
else
    f2 = figure('Name','Amplitude at DOFs');
end
figs = [f1, f2];
colors = get(0,'defaultaxescolororder');
totalComputationTime = zeros(size(ORDER));

for j = 1:numel(ORDER)
        

    order = ORDER(j);
    startFRC = tic;
    if isempty(obj.System.spectrum)
        [~,~,~] = obj.System.linear_spectral_analysis();
    end
    lambda  = obj.System.spectrum.Lambda;
    assert(~isreal(lambda),'One or more eigenvalues must be underdamped for FRC computation using SSMs')

    % detect resonant eigenvalues in the parameter range
    switch lower(parName)
        case 'freq'
            if ~isempty(obj.System.fext)
                assert(~isempty(obj.System.fext.epsilon), 'The epsilon field is empty in the dynamical system external forcing');
            else
                assert(~isempty(obj.System.Fext.epsilon), 'The epsilon field is empty in the dynamical system external forcing');
            end
            [resLambda,resFreq] = find_eigs_in_freq_range(parRange,lambda,obj.FRCOptions.resType);

            % obtain subintervals around each resonant eigenvalue
            [parNodes, nSubint] = subdivide_freq_range(parRange, resFreq);
        case 'amp'
            assert(~isempty(obj.System.Omega), 'The Omega field is empty in the dynamical system');
            % find eigenvalue nearest to forcing frequency
            [~,idx] = min(abs(lambda - 1i*obj.System.Omega));
            resLambda = lambda(idx);
            % setup for loop
            parNodes = [parRange(1) parRange(end)];
            nSubint = 1;
    end

    FRC = cell(nSubint,1);
    for i=1:nSubint
        

        % tune subinterval
        parSubRange = parNodes(i:i+1)';
        parSubRange = tune_parameter_range(parSubRange, obj.FRCOptions.frac, i, nSubint);

        % detect modes resonant with resLambda(i)
        [resModes,mFreqs] = detect_resonant_modes(resLambda(i),lambda, obj.Options.IRtol);
        

FRC computation within the subinterval

        disp('*****************************************');
        disp(['Calculating FRC using SSM with master subspace: [' num2str(resModes(:).') ']']);

        switch obj.FRCOptions.method
            case 'level set'
                FRC{i} = obj.FRC_level_set(resModes,order,parName,parSubRange);
                plotStyle = 'circles';
            case 'continuation ep'
                % call continuation based method
                mFreqs = mFreqs(1:2:end)';
                runid  = ['freqSubint',num2str(i)];
                if j>1 % take lowest order solution as initial guess
                    sol_jminus1 = ep_read_solution('', [runid,'.ep'], 1);
                    obj.FRCOptions.p0 = sol_jminus1.p;
                    obj.FRCOptions.z0 = sol_jminus1.x;
                end
                FRC{i}   = obj.FRC_cont_ep(runid,resModes,order,mFreqs,parName,parSubRange);
                plotStyle = 'lines';
            case 'continuation po'
                runid  = ['freqSubint',num2str(i)];
                FRCi = obj.FRC_cont_po(runid,resModes,order,parSubRange);
                FRC{i}  = cat(2,FRCi{:});
                plotStyle = 'circles';
        end
        

    end
    % concatenate cell contents as struct arrays
    FRC = cat(1,FRC{:});
    totalComputationTime(j) = toc(startFRC);
        

Plot FRC in physical coordinates

    plot_FRC(FRC,obj.FRCOptions.outdof,order,parName,plotStyle, figs, colors(j,:))
        

end


for j = 1:numel(ORDER)
    disp(['Total time spent on FRC computation upto O(' num2str(ORDER(j)) ') = ' datestr(datenum(0,0,0,0,0,totalComputationTime(j)),'HH:MM:SS')])
end
        

end

function parRange = tune_parameter_range(parRange, frac, i, nSubint)
% amplify the parameter subinterval except on the first and the last nodes.
if i>1
    parRange(1) = frac(1)*parRange(1);
end
if i<nSubint
    parRange(2) = frac(2)*parRange(2);
end
end

function [resLambda, resFreq] = find_eigs_in_freq_range(Omega,lambda,resType)

switch resType
    case '2:1'  % Case of subharmonic resonance

    natFreq = imag(lambda);
    resFreqID = intersect(find(2*natFreq>Omega(1)), find(2*natFreq<Omega(end)));
    resFreq = natFreq(resFreqID);
    resLambda = lambda(resFreqID);

    assert(~isempty(resFreq),'Input frequency range should include dual multiple of natural frequency');

    case '1:1'

    natFreq = imag(lambda);
    resFreqID = intersect(find(natFreq>Omega(1)), find(natFreq<Omega(end)));
    resFreq = natFreq(resFreqID);
    resLambda = lambda(resFreqID);
    assert(~isempty(resFreq),'Input frequency range should include at least one natural frequency');
end

end

function [freqNodes, nSubint] = subdivide_freq_range(parRange,resFreq)
nSubint = numel(resFreq);
freqNodes = 0.5*(resFreq(1:end-1)+resFreq(2:end)); % center points of two adjacent resonant modes
freqNodes = [parRange(1); freqNodes; parRange(2)];
end
        

Published with MATLAB® R2023a