# Compute steady of ODE

7 views (last 30 days)
Erik Eriksson on 28 May 2024
Commented: Sam Chak on 30 May 2024
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!
##### 2 CommentsShow NoneHide None
Torsten on 28 May 2024
Is your dcdt always a polynomial in the variable c ?
Sam Chak on 28 May 2024
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.

SAI SRUJAN on 28 May 2024
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.
I hope this helps!
Sam Chak on 30 May 2024
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.

### Categories

Find more on Programming in Help Center and File Exchange

R2023b

### Community Treasure Hunt

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

Start Hunting!