# Compute steady of ODE

7 views (last 30 days)

Show older comments

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 Comments

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.

### Answers (1)

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!

##### 1 Comment

Sam Chak
on 30 May 2024

### See Also

### Community Treasure Hunt

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

Start Hunting!