## Compute the basin of attraction for a steady state of an ordinary differential equation ODE

I need write a function called "compute_basin" that takes in a steady state and its k value and computes which set of initial conditions in [15,30] approach that steady state. Specifically it will return two values, "a" and "b", which describe the set [a,b] of initial conditions that approach the steady state. The function will be assessed by running

[a,b]=compute_basin(steady_state,k);

for many different values of k and the associated steady states.

Inputs:

a scalar value of a steady state

a scalar value of k (note the steady state should be correct for the k value)

Outputs:

a scalar for the lower bounds of the basin of attraction

a scalar for the upper bounds for the basin of attraction

This is what i have so far:

function [a, b] = compute_basin(Ceq, k)

% Initialize variables

a = 15; % Lower bound for the basin of attraction

b = 30; % Upper bound for the basin of attraction

threshold = 0.001; % Error threshold

% Binary search for the lower bound (a)

left = 15;

right = 30;

while left < right

mid = (left + right) / 2;

tstar = compute_time_to_steady(mid, k, Ceq);

if tstar > 0

right = mid;

else

left = mid + 0.1;

end

end

a = left;

% Binary search for the upper bound (b)

left = 15;

right = 30;

while left < right

mid = (left + right) / 2;

tstar = compute_time_to_steady(mid, k, Ceq);

if tstar > 0

left = mid;

else

right = mid – 0.1;

end

end

b = right;

end

But unfortunalty this is not correct. I have a tip saying that "There are many ways of doing this but a binary search algorithm might be helpful (Wikipedia entry for binary searchLinks to an external site.). Make sure that the [a,b] bounds are within the range of [15,30]. Note that some steady states may be unstable and so their basin of attraction will only be their exact value. In these cases a=b. It helps to solve for the upper and lower bounds separately, i.e. running the code for a and then running it again for b."

I have acces to three functions, the first is "compute_ode" that computes the derivative dC/dt for the differential equation model dcdt = 0.1 * (c-20) * (23-c) * (c-26) + k.

function dcdt = compute_ode(t,c,k)

% write the equation for the differential equation and assign it to output dcdt

dcdt = 0.1 * (c-20) * (23-c) * (c-26) + k;

end

Then i have the second "compute_states" that takes as input a value for k and then returns as output all real-valued steady states for the differential equation.

function vectCeq = compute_states(k)

% Coefficients of the polynomial equation -0.1*c^3 + 6.9*c^2 – 157.8*c + 1196 + k = 0

coeffs = [-0.1, 6.9, -157.8, 1196 + k];

% Find all roots of the polynomial

all_roots = roots(coeffs);

% Filter out complex roots

real_roots = all_roots(imag(all_roots) == 0);

% Return the real-valued steady states

vectCeq = real_roots;

end

And latsly, I have "compute_time_to_steady" that solves the differential equation and returns the earliest time it takes to get within 1% error of the steady state value:

function tstar = compute_time_to_steady(init_cond, k, steady_state)

% Set the options for the ODE solver

options = odeset(‘RelTol’, 1e-6, ‘AbsTol’, 1e-9);

% Define the ODE function as a nested function

function dcdt = ode_function(t, c)

dcdt = compute_ode(t, c, k);

end

% Initialize variables

tspan = [0 100000]; % Time range from 0 to 100,000

threshold = 0.001; % 1% error threshold

% Solve the ODE

[T, C] = ode45(@ode_function, tspan, init_cond, options);

% Iterate over the solution to find the first time it meets the criterion

for i = 1:length(T)

if abs(C(i) – steady_state) / steady_state < threshold

tstar = T(i);

return;

end

end

% If the criterion is not met within the time range, return the last time point

tstar = T(end);

end

function dcdt = compute_ode(t, c, k)

dcdt = 0.1 * (c – 20) * (23 – c) * (c – 26) + k;

end

function steady_states = compute_states(k)

% Define the steady state values based on the value of k

steady_states = [20; 23; 26];

end

I would really appriciate any help on computing the basin of attraction for a steady state!

Thank you!I need write a function called "compute_basin" that takes in a steady state and its k value and computes which set of initial conditions in [15,30] approach that steady state. Specifically it will return two values, "a" and "b", which describe the set [a,b] of initial conditions that approach the steady state. The function will be assessed by running

[a,b]=compute_basin(steady_state,k);

for many different values of k and the associated steady states.

Inputs:

a scalar value of a steady state

a scalar value of k (note the steady state should be correct for the k value)

Outputs:

a scalar for the lower bounds of the basin of attraction

a scalar for the upper bounds for the basin of attraction

This is what i have so far:

function [a, b] = compute_basin(Ceq, k)

% Initialize variables

a = 15; % Lower bound for the basin of attraction

b = 30; % Upper bound for the basin of attraction

threshold = 0.001; % Error threshold

% Binary search for the lower bound (a)

left = 15;

right = 30;

while left < right

mid = (left + right) / 2;

tstar = compute_time_to_steady(mid, k, Ceq);

if tstar > 0

right = mid;

else

left = mid + 0.1;

end

end

a = left;

% Binary search for the upper bound (b)

left = 15;

right = 30;

while left < right

mid = (left + right) / 2;

tstar = compute_time_to_steady(mid, k, Ceq);

if tstar > 0

left = mid;

else

right = mid – 0.1;

end

end

b = right;

end

But unfortunalty this is not correct. I have a tip saying that "There are many ways of doing this but a binary search algorithm might be helpful (Wikipedia entry for binary searchLinks to an external site.). Make sure that the [a,b] bounds are within the range of [15,30]. Note that some steady states may be unstable and so their basin of attraction will only be their exact value. In these cases a=b. It helps to solve for the upper and lower bounds separately, i.e. running the code for a and then running it again for b."

I have acces to three functions, the first is "compute_ode" that computes the derivative dC/dt for the differential equation model dcdt = 0.1 * (c-20) * (23-c) * (c-26) + k.

function dcdt = compute_ode(t,c,k)

% write the equation for the differential equation and assign it to output dcdt

dcdt = 0.1 * (c-20) * (23-c) * (c-26) + k;

end

Then i have the second "compute_states" that takes as input a value for k and then returns as output all real-valued steady states for the differential equation.

function vectCeq = compute_states(k)

% Coefficients of the polynomial equation -0.1*c^3 + 6.9*c^2 – 157.8*c + 1196 + k = 0

coeffs = [-0.1, 6.9, -157.8, 1196 + k];

% Find all roots of the polynomial

all_roots = roots(coeffs);

% Filter out complex roots

real_roots = all_roots(imag(all_roots) == 0);

% Return the real-valued steady states

vectCeq = real_roots;

end

And latsly, I have "compute_time_to_steady" that solves the differential equation and returns the earliest time it takes to get within 1% error of the steady state value:

function tstar = compute_time_to_steady(init_cond, k, steady_state)

% Set the options for the ODE solver

options = odeset(‘RelTol’, 1e-6, ‘AbsTol’, 1e-9);

% Define the ODE function as a nested function

function dcdt = ode_function(t, c)

dcdt = compute_ode(t, c, k);

end

% Initialize variables

tspan = [0 100000]; % Time range from 0 to 100,000

threshold = 0.001; % 1% error threshold

% Solve the ODE

[T, C] = ode45(@ode_function, tspan, init_cond, options);

% Iterate over the solution to find the first time it meets the criterion

for i = 1:length(T)

if abs(C(i) – steady_state) / steady_state < threshold

tstar = T(i);

return;

end

end

% If the criterion is not met within the time range, return the last time point

tstar = T(end);

end

function dcdt = compute_ode(t, c, k)

dcdt = 0.1 * (c – 20) * (23 – c) * (c – 26) + k;

end

function steady_states = compute_states(k)

% Define the steady state values based on the value of k

steady_states = [20; 23; 26];

end

I would really appriciate any help on computing the basin of attraction for a steady state!

Thank you! I need write a function called "compute_basin" that takes in a steady state and its k value and computes which set of initial conditions in [15,30] approach that steady state. Specifically it will return two values, "a" and "b", which describe the set [a,b] of initial conditions that approach the steady state. The function will be assessed by running

[a,b]=compute_basin(steady_state,k);

for many different values of k and the associated steady states.

Inputs:

a scalar value of a steady state

a scalar value of k (note the steady state should be correct for the k value)

Outputs:

a scalar for the lower bounds of the basin of attraction

a scalar for the upper bounds for the basin of attraction

This is what i have so far:

function [a, b] = compute_basin(Ceq, k)

% Initialize variables

a = 15; % Lower bound for the basin of attraction

b = 30; % Upper bound for the basin of attraction

threshold = 0.001; % Error threshold

% Binary search for the lower bound (a)

left = 15;

right = 30;

while left < right

mid = (left + right) / 2;

tstar = compute_time_to_steady(mid, k, Ceq);

if tstar > 0

right = mid;

else

left = mid + 0.1;

end

end

a = left;

% Binary search for the upper bound (b)

left = 15;

right = 30;

while left < right

mid = (left + right) / 2;

tstar = compute_time_to_steady(mid, k, Ceq);

if tstar > 0

left = mid;

else

right = mid – 0.1;

end

end

b = right;

end

But unfortunalty this is not correct. I have a tip saying that "There are many ways of doing this but a binary search algorithm might be helpful (Wikipedia entry for binary searchLinks to an external site.). Make sure that the [a,b] bounds are within the range of [15,30]. Note that some steady states may be unstable and so their basin of attraction will only be their exact value. In these cases a=b. It helps to solve for the upper and lower bounds separately, i.e. running the code for a and then running it again for b."

I have acces to three functions, the first is "compute_ode" that computes the derivative dC/dt for the differential equation model dcdt = 0.1 * (c-20) * (23-c) * (c-26) + k.

function dcdt = compute_ode(t,c,k)

% write the equation for the differential equation and assign it to output dcdt

dcdt = 0.1 * (c-20) * (23-c) * (c-26) + k;

end

Then i have the second "compute_states" that takes as input a value for k and then returns as output all real-valued steady states for the differential equation.

function vectCeq = compute_states(k)

% Coefficients of the polynomial equation -0.1*c^3 + 6.9*c^2 – 157.8*c + 1196 + k = 0

coeffs = [-0.1, 6.9, -157.8, 1196 + k];

% Find all roots of the polynomial

all_roots = roots(coeffs);

% Filter out complex roots

real_roots = all_roots(imag(all_roots) == 0);

% Return the real-valued steady states

vectCeq = real_roots;

end

And latsly, I have "compute_time_to_steady" that solves the differential equation and returns the earliest time it takes to get within 1% error of the steady state value:

function tstar = compute_time_to_steady(init_cond, k, steady_state)

% Set the options for the ODE solver

options = odeset(‘RelTol’, 1e-6, ‘AbsTol’, 1e-9);

% Define the ODE function as a nested function

function dcdt = ode_function(t, c)

dcdt = compute_ode(t, c, k);

end

% Initialize variables

tspan = [0 100000]; % Time range from 0 to 100,000

threshold = 0.001; % 1% error threshold

% Solve the ODE

[T, C] = ode45(@ode_function, tspan, init_cond, options);

% Iterate over the solution to find the first time it meets the criterion

for i = 1:length(T)

if abs(C(i) – steady_state) / steady_state < threshold

tstar = T(i);

return;

end

end

% If the criterion is not met within the time range, return the last time point

tstar = T(end);

end

function dcdt = compute_ode(t, c, k)

dcdt = 0.1 * (c – 20) * (23 – c) * (c – 26) + k;

end

function steady_states = compute_states(k)

% Define the steady state values based on the value of k

steady_states = [20; 23; 26];

end

I would really appriciate any help on computing the basin of attraction for a steady state!

Thank you! basin, ode, ode45, differential equations MATLAB Answers — New Questions