Compute steady of ODE

First of all I have a function called "compute_ode" that computes the derivative dC/dt for the differential equation model. The function will take in three inputs and return one output. The function will be assessed by running the code: "[t,c]=ode45(@(t,c) compute_ode(t,c,k),t_range,init_cond);" for chosen values of k, t_range, and init_cond. Thus, the function should conform to the standards required by the MATLAB function ode45.

The code for that is:

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

and to call that:

t_range=[0:.1:10];
init_cond=18;
k=0; % here, you can change the k value
[t,c]=ode45(@(t,c) compute_ode(t,c,k),t_range,init_cond);
plot(t,c,'LineWidth',3)
xlabel('time','FontSize',18)
ylabel('c(t)','FontSize',18)

NOW I will write a function called "compute_states" that takes as input a value for k and then returns as output ALL real-valued steady states for the ODE. The function will be assessed by running "Ceq = compute_states(k)" for many different values of k. I also want to be sure to filter out any complex-valued steady states. I know that "roots" and "imag" might be helpful.

Input: a scalar for k

Output: a vector of steady states, where the lenght of the vector is the number of all real-valued steady states.

ODE: dC/dt = 0.1*(C-20)*(23-C)*(C-26) +k

What I have so far:

Function:

function vectCeq = compute_states(k)
    % note that vectCeq is a vector of all real-valued steady states
    

end

and to call my function:

k=0; % try for different values of k
vectCeq = compute_states(k)

I know that "roots" and "imag" commands might be helpful.

Does anyone have any suggestion or now how to solve my next step for the "compute_states"?

Thank you!

[Torsten]

Is your dcdt always a polynomial in the variable c ?

[Sam Chak]

@Erik Eriksson,

Although the ODE function 0.1*(C - 20)*(23 - C)*(C - 26) + k is a polynomial, you should aim to solve for general nonlinear functions and make certain assumptions about the system.

Despite the fact that you can numerically find the real-valued equilibrium points, you cannot tell whether those equilibrium points are stable or unstable unless you perform further analysis. Therefore, you cannot determine which value the trajectory is asymptotically converging to.

Answer 1

[SAI SRUJAN]

Hi Erik,

I understand that you are facing an issue implementing the 'compute_states' function.

To find the steady states of the given ODE, we need to solve the equation for when the derivative (dC/dt = 0). This translates to solving the polynomial equation (0.1(C-20)(23-C)(C-26) + k = 0). To do this in MATLAB, we can express the polynomial in its expanded form and then use the 'roots' function to find its roots.

Please follow the below code sample to proceed further,

function vectCeq = compute_states(k)
    % Coefficients of the polynomial in descending powers
    % The polynomial is 0.1*(C^3 - 69*C^2 + (20*23 + 23*26 + 20*26)*C - 20*23*26) + k
    % Simplified: 0.1*C^3 - 6.9*C^2 + 139.8*C - 1188 + k
    % We need to expand the polynomial and include k in the constant term
    coeffs = [0.1, -6.9, 139.8, -1188 + k];
    
    % Find the roots of the polynomial
    rootsC = roots(coeffs);
    
    % Filter only the real-valued roots
    realRoots = rootsC(imag(rootsC) == 0);
    
    % Assign the real-valued roots to the output
    vectCeq = realRoots;
end

For a comprehensive understanding of the 'roots' and 'imag' functions in MATLAB, please refer to the following documentation.

https://www.mathworks.com/help/matlab/ref/roots.html

https://www.mathworks.com/help/matlab/ref/imag.html

I hope this helps!

[Sam Chak]

Hi @Erik Eriksson

If you're implementing @SAI SRUJAN's solution in your work, kindly consider clicking 'Accept' ✔ on the answer to resolve the issue. However, I suggest refining the solution to directly utilize the 'compute_ode(t, c, k)' function, eliminating the need for manually expanding the cubic polynomial.