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

22 次查看（过去 30 天）
Erik Eriksson
Erik Eriksson 2024-5-28
回答： Nipun 2024-6-3
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!
  5 个评论
显示 3更早的评论
Erik Eriksson
Erik Eriksson 2024-5-30
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
Erik Eriksson 2024-5-30
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Nipun
Nipun 2024-6-3
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

类别

在 Help Center 和 File Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by