Compute steady of ODE
6 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Sam Chak
el 28 de Mayo de 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.
Respuestas (1)
SAI SRUJAN
el 28 de Mayo de 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 comentario
Sam Chak
el 30 de Mayo de 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.
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!