Solve a differential equation and determine when it gets close to steady state

I have a function called "compute_ode" that computes the derivative dC/dt for the differential equation model ODE: dC/dt = 0.1*(C-20)*(23-C)*(C-26) +k. 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 function looks like this:

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 I am running it with:

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)

Then I wrote 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 differential equation. The function will be assessed by running

ceq = compute_states(k)

for different values of k. Looking like this:

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

Now I need to write a function called "compute_time_to_steady" that solves the differential equation and returns the earliest time (t) it takes to get within 1% error of the steady state value. In other words, if the computed value is y(t) at time t and the steady state is y* then the stopping criterion is abs(y(t)-y*)/y*<0.001 and the code returns the first time t that satisfies this criterion. Importantly, the function assumes a steady state exists and that the stopping criterion should be reached before t=100,000. If the stopping criterion is not satisfied by t=100,000 then the code should return -1, i.e. csol=-1. The function will be assessed by running the code:

csol = compute_time_to_steady(init_cond,k,steady_state);

for many different values of init_cond (a scalar), k (a scalar), and steady_state (a scalar).

Any suggestion on how to solve this? I have got some guisng saying that "find " might be useful.

Thanks in advance!

Answer 1

[Torsten]

You can use an event function for the ODE integrator and stop integration when abs(dy/dt) (in your case abs(0.1 * (c-20) * (23-c) * (c-26) + k)) becomes small.

https://uk.mathworks.com/help/matlab/math/ode-event-location.html

[Mike Croucher]

I also suggest moving to the new ODE framework The new solution framework for Ordinary Differential Equations (ODEs) in MATLAB R2023b » The MATLAB Blog - MATLAB & Simulink (mathworks.com)

I'm curious if the new SUNDIALS integration would be faster in this case too. Faster Ordinary Differential Equations (ODEs) solvers and Sensitivity Analysis of Parameters: Introducing SUNDIALS support in MATLAB » The MATLAB Blog - MATLAB & Simulink (mathworks.com)

[Sam Chak]

Hi @Erik Eriksson

When I modify the ODE, the function compute_states(k) no longer computes the states correctly. Could you please fix the code you wrote so that the function compute_states(k) depends on the function compute_ode(t, c, 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 - 30)*(33 - c)*(c - 36) + k;
end

k       = 0;
vectCeq = compute_states(k)

vectCeq = 3x1

26.0000 23.0000 20.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

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

[Sam Chak]

@Erik Eriksson, Can you apply this formula when c(t) get close to 23?

[Sam Chak]

Hi @Mike Croucher

Thank you for your suggestion. Did I compare the results correctly? I believe that the SUNDIALS solver generally provides fairly accurate solutions with the basic settings, although the simulation takes a slightly longer time to run. Despite the fact that the ode45 solver performs faster with the basic settings, it returns some irregularities for certain trajectories. Those trajectories should be within the expected exponential decay envelopes.

We can also observe some minor differences between the two solvers at time for the green and orange trajectories that begin with the initial values of .

%% Solve ODE using SUNDIALS solver
function solveWithSUNDIALS
    F               = ode;
    F.ODEFcn        = @(t, c, k) 0.1*(c - 20)*(23 - c)*(c - 26) + k;
    F.Parameters    = 0;
    F.Solver        = "cvodesNonstiff";
    
    c0  = 18:0.125/2:28;
    for j = 1:numel(c0)
        F.InitialValue  = c0(j);
        sol             = solve(F, 0, 6); 
        t               = sol.Time;
        c               = sol.Solution;
        figure(1), plot(t, c), hold on
    end
    grid on, ylim([18 28]), hold off
end

%% Solve ODE using ode45 solver
function solveWithODE45
    ODEFcn        = @(t, c, k) 0.1*(c - 20)*(23 - c)*(c - 26) + k;
    
    c0  = 18:0.125/2:28;
    for j = 1:numel(c0)
        [t, c]  = ode45(@(t, c) ODEFcn(t, c, 0), [0, 6], c0(j)); 
        figure(2), plot(t, c), hold on
    end
    grid on, ylim([18 28]), hold off
end

%% Test 1 with SUNDIALS
f1  = @() solveWithSUNDIALS;
T1  = timeit(f1)

T1 = 0.3345

%% Test 2 with ode45
f2  = @() solveWithODE45;
T2  = timeit(f2)

T2 = 0.2553

[Sam Chak]

Hi @Mike Croucher

I have conducted further testing to find the analytical solution for the first-order ODE subject to an initial value of . It turns out that the result, represented by the green trajectory in Figure 2, returned by the ode45 solver was more accurate than the SUNDIALS solver.

syms c(t)
eqn     = diff(c,t) == 0.1*(c - 20)*(23 - c)*(c - 26);
cond    = c(0)==23.0625;
cSol(t) = dsolve(eqn, cond)

cSol(t) = 

fplot(cSol, [0 6]), grid on, xlabel t, ylabel c(t)
xline( 4, '-.', 't = 4');
yline(25, '-.', 'Marker');
title('Result generated by Analytical Solution')