Arc Length Continuation Method: Numerical Method

I am trying plot a curve (frequency vs amplitude) that has a (sn) birfurcation, and I am trying to write a code that would use the Arc Length Continuation to get the solve the nonlinear equations and get the (frequency response) curve (include regions where two solutions exist).

The methodology is as follows

• Initialization and Initial Guess: Find an initial guess for a specific omega (parameter)
• Residual Calculation: Using a function handle to calculate residuals, based on the nonlinear equation
• Newton-Raphson Solver: Solves the nonlinear equation using the Newton-Raphson method.
• Continuation Method: Continuation method uses the arc length continuation to trace the solution curve by extrapolating new guesses and solving using the Newton-Raphson method.
• Arc Length Constraint: Ensures the next step in the continuation process is a fixed arc length from the current solution.
• Plotting: Code for plotting the FRF is included but commented out because the code is incomplete

My issue is in the function called continuation--which continunes to solve for nSteps (number of steps), by calculating the new guess, and augemeniting the Jacobain to have the new constraint equation--I believe I have the correct(?) logic but my code is incorrect.

(I know there are straight-forward ways to obtain the FRF for the function defined in Analytic Equation, the reason I am testing it out on a straight-forward equation is to make sure it is running correctly.

Any pointers appreciated

clc; clear; close all

w = [0]; % Path-finding parameters,omega--w, for initial search; starting at 0 Hz
x0 = [0;0]; % Initial guess (x(1) is the solution to the cosine and x(2) is the solution to the sine)
nSteps = 1000; % Number of steps to continue forward on the curve

% Initialize q_val
q_val = zeros(length(x0)+1,1); %this is the solution vector, it will have [x(1); x(2); w]

% Find the residuals using the initial w and the equation
equation_f = @(x, w) equation(x, w); % Function to calculate the residuals for each forcing

[sol, iter] = nlsolver(equation_f, x0, w, 10e-10); % Using Newton-Raphson to solve the equation
q_val = [sol;w]; % Store solution and corresponding w

% Perform continuation--trace the solution curve by extrapolating new guesses and solving using the Newton-Raphson method.
continuation_combined(q_val, nSteps)

% Plotting the results
% w = q_val(3,:);
% a = sqrt(q_val(1,:).^2 + q_val(2,:).^2);
% plot(w, a); grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% CONTINUATION FUNCTIONS

function q_val = continuation_combined(q, N)
q_val = q;

    for n = 1:N
        xguess = 2 * q_val(:, end) - q_val(:, end-1); % Extrapolate a new guess for the solution
        xguess = nlsolver(continuation_side(xguess), xguess); % Solve using the extrapolated guess
        q_val = [q_val, xguess]; % Append the new solution
    end

    function new_q = continuation_side(q_val)
        w = q_val(end); % Extract the current parameter value
        x = q_val(1:end-1); % Extract the state variables
        
        % Constants
        k = 1; m = 1; c = 0.02; f0 = 1; F = [f0; 0]; kc = 0.5;

        % Linear part of the equation
        A = [k - w^2 * m, -c * w; c * w, k - w^2 * m];
        
        % Nonlinear part of the equation
        B = [x(1) * (x(2)^2 + x(1)^2); x(2) * (x(2)^2 + x(1)^2)];
        
        % Residuals
        R = A * x + 0.75 * kc * B - F;
        
        % Arc length constraint
        S = 0.01;
        aa = norm(q_val(:, end) - q_val) - S;
        
        new_q = [R; aa]; % Return the combined residual and arc length constraint
    end
end

%%%%% ANALYTIC EQUATION (OBTAINED FROM HBM)
function [R] = equation(x, w)
    % Constants
    k = 1; m = 1; c = 0.02; f0 = 1; F = [f0; 0]; kc = 0.5;
    
    % Linear part of the equation
    A = [k - w^2 * m, -c * w; c * w, k - w^2 * m];
    
    % Nonlinear part of the equation
    B = [x(1) * (x(2)^2 + x(1)^2); x(2) * (x(2)^2 + x(1)^2)];
    
    % Residuals
    R = A * x + 0.75 * kc * B - F;
end

%%%%% NEWTON-RAPHSON SOLVER
function [sol, iter] = nlsolver(func, guess, w, errorbound)
    % Initializations
    norm_R0 = 1;
    iteration = 0;
    hj = 1e-8; % Step size
    x = guess;

    while norm_R0 > errorbound && iteration < 1000 % Either error bound satisfied or iteration hits max threshold for each point
        R0 = func(x, w); % Evaluate the function at x (initial conditions)
        J = zeros(length(x), length(x)); % Initialize the Jacobian matrix
        
        for r = 1:length(x)
            for c = 1:length(x)
                ej = zeros(length(x), 1);
                ej(c) = 1;
                Ri = func(x + hj * ej, w); % Evaluate the function at x + hj * ej
                J(r, c) = (Ri(r) - R0(r)) / hj; % Compute the Jacobian entry
            end
        end
        
        x = x - J \ R0; % Update x using the Jacobian and function values
        iteration = iteration + 1;
        norm_R0 = norm(R0); % Update the norm of R0
    end
    
    sol = x;
    iter = iteration;
end

Answer 1

[Umar]

Hi Cesium,

You need to modify the continuation_combined function to call the nlsolver function correctly. Instead of passing continuation_side as the first argument to nlsolver, we should pass the equation_f function, which calculates the residuals based on the nonlinear equation.

Here's the corrected code for the continuation_combined function:

function q_val = continuation_combined(q, N)

    q_val = q;
    for n = 1:N
        xguess = 2 * q_val(:, end) - q_val(:, end-1); % Extrapolate a new guess for the solution

        xguess = nlsolver(equation_f, xguess, q_val(end), 10e-10); % Solve using the extrapolated guess

        q_val = [q_val, xguess]; % Append the new solution
    end
end

With this modification, the nlsolver function will correctly solve the nonlinear equation using the provided extrapolated guess.

[Cesium Modern]

Hi Umar, thanks for the answer. I believe the function inputs should also be modified because I am getting the following error:

Unrecognized function or variable 'equation_f'.

Error in arclength>continuation_combined (line 38)

xguess = nlsolver(equation_f, xguess, q_val(end), 10e-10); % Solve using the extrapolated guess

Error in arclength (line 20)

continuation_combined(q_val,nSteps)

How would you parse in the function? feval

[Umar]

Hi Cesium,

To parse the function equation_f correctly in Matlab, you need to make sure that the function is defined or available in the current workspace or script. Check if the function equation_f is defined before calling it in the nlsolver function. You can define the function equation_f in the same script or in a separate function file and ensure it is accessible where it is being called. Verify that the function name is spelled correctly and that there are no syntax errors within the function definition. Once the function equation_f is correctly defined and accessible, the error should be resolved.

Answer 2

[Walter Roberson]

equation_f = @(x, w) equation(x, w); % Function to calculate the residuals for each forcing
[sol, iter] = nlsolver(equation_f, x0, w, 10e-10); % Using Newton-Raphson to solve the equation

That part looks valid: you are calling nlsolver passing in a function handle as the first parameter.

        xguess = nlsolver(continuation_side(xguess), xguess); % Solve using the extrapolated guess

That calls continuation_side with the parameter xguess, and expects that the result will be a function handle to pass into nlsolver. But continuation_side returns a numeric array.