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

​