Effacer les filtres
Effacer les filtres

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

34 vues (au cours des 30 derniers jours)
Erik Eriksson
Erik Eriksson le 28 Mai 2024 à 13:23
Commenté : Sam Chak le 29 Mai 2024 à 3:40
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!

Réponses (1)

Torsten
Torsten le 28 Mai 2024 à 13:26
Modifié(e) : Torsten le 28 Mai 2024 à 13:27
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.
  5 commentaires
Sam Chak
Sam Chak le 29 Mai 2024 à 3:25
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
Sam Chak le 29 Mai 2024 à 3:40
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')

Connectez-vous pour commenter.

Produits


Version

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by