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:

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

Outputs:

1. a scalar for the lower bounds of the basin of attraction
2. 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!

[Sam Chak]

Hi @Erik Eriksson

The cubic polynomial expression "" represents a downhill curve. Thus, there could be one of the following scenarios:

1. one stable equilibrium point .
2. One stable equilibrium point with two unstable equilibrium points .
3. Two stable equilibrium points with an unstable equilibrium point sandwiched in between .

I am curious to understand how there can be three basins of attraction for the computed steady states [20, 23, 26] when .

[Erik Eriksson]

Okay, how did you determine these values for k?

[Sam Chak]

Hi @Erik Eriksson

The values for k can be obtained through simulations to observe the number of stable and unstable equilibrium points. By the way, have you resolved the basin problem?

[Erik Eriksson]

Okay! No I have not, for my code i have got some hints on solving it. Here is the hint from my supervisor:

"My hint for this is to try to find the bounds separately. You are trying to find both the upper and lower bound at the exact same time. This causes problems. You know the lower bound for the basin of attraction is somewhere between 15 and the steady state. You know the upper bound for the base of traction is somewhere between 30 in the steady state. So these two problems are very similar and you solve each one separately".

Would you have any idea on how to solve this. Would be very appriciated!

[Erik Eriksson]

I managed to wolve it!

Here is the code:

function [a, b] = compute_basin(Ceq, k)
    % Initialize variables
    threshold = 0.001; % Error threshold
   
    % Binary search for the lower bound (a)
    left = 15;
    right = Ceq;
    while left < right - 0.1
        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 = max(15, left); % Ensure a is within [15, Ceq]
   
    % Binary search for the upper bound (b)
    left = Ceq;
    right = 30;
    while left < right - 0.1
        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 = min(30, right); % Ensure b is within [Ceq, 30]
    
    % Handle cases where the steady state is unstable
    if a == b
        a = Ceq;
        b = Ceq;
    end
end

Answer 1

[Nipun]

Hi Erik,

I understand that you need to write a function compute_basin to find the set of initial conditions that approach a given steady state for a differential equation. You want to find the lower and upper bounds of this set using a binary search algorithm.

Here's a revised version of your function that should help you achieve this:

function [a, b] = compute_basin(Ceq, k)
    % Initialize variables
    threshold = 0.01; % 1% error threshold for the steady state
    tol = 1e-3;       % Tolerance for binary search

    % Binary search for the lower bound (a)
    left = 15;
    right = Ceq;
    while (right - left) > tol
        mid = (left + right) / 2;
        tstar = compute_time_to_steady(mid, k, Ceq);
        if abs(tstar - Ceq) / Ceq < threshold
            right = mid;
        else
            left = mid;
        end
    end
    a = left;

    % Binary search for the upper bound (b)
    left = Ceq;
    right = 30;
    while (right - left) > tol
        mid = (left + right) / 2;
        tstar = compute_time_to_steady(mid, k, Ceq);
        if abs(tstar - Ceq) / Ceq < threshold
            left = mid;
        else
            right = mid;
        end
    end
    b = right;
end

This function uses binary search to find the lower (a) and upper (b) bounds of the basin of attraction. It adjusts the left and right bounds based on whether the solution approaches the steady state within a 1% error threshold.

Hope this helps.

Regards,

Nipun