## 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;

endI 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 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 matlab MATLAB Answers — New Questions