Why doesn't the ode15s solver work when I switch 'k1q','k2q','k12q' from a single value into vectors(more values)?

1 visualización (últimos 30 días)
clear;clc;close all;
t_measured = 0:5; % measured times
k1 = [0.1,0.2,0.3,0.4,0.5,0.6]; % k1 at the measured times
k2 = [0.2,0.4,0.6,0.8,1.0,1.2]; % k2 at the measured times
k12 = [0.02,0.04,0.06,0.08,0.1,0.12]; % k12 at the measured times
tq = 0:0.1:5; % times for interpolation (more points than the measured times)
k1q = interp1(t_measured, k1, tq); % k1 vaues after interpolation
k2q = interp1(t_measured, k2, tq); % k2 vaues after interpolation
k12q = interp1(t_measured, k12, tq); % k12q vaues after interpolation
% k1q = 1;
% k2q = 2;
% k12q = 3;
initialconcentrations = [1 0]; % initial concentratioins of soil and root
tol = 1e-13; % tolerance
options = odeset ('RelTol', tol,'AbsTol',[tol tol], 'NonNegative',1, 'NonNegative',2);
[t,C] = ode15s(@(t,C) solver_test2(t,C,k1q,k2q,k12q),tq,initialconcentrations,options)
function dCdt = solver_test2(t,C,k1q,k2q,k12q)
dCdt = zeros(2,1);
dCdt(1) = -k1q .* C(1);
dCdt(2) = -k2q .* C(2) + k12q .* C(1);
end

Respuesta aceptada

Bjorn Gustavsson
Bjorn Gustavsson el 16 de Abr. de 2021
Editada: Bjorn Gustavsson el 17 de Abr. de 2021
When k1q is an array you will get an error when you try to assign the product of that array with the scalar C(1) to dCdt(1) which is expected to be a scalar. Your ODE-function is expected to return the value of dCdt as a 2-by-1 array at any time it queried in ode15s. The important thing to remember when using the ODEnn-functions is that they use intricate schemes to step forward in time with tentative partial steps of different lengths in some order that we as users should consider to be "rather random" when it comes into designing our ODE-functions. When you use your arrays as input to the function there is no explicit knowledge for what values out of those that should be used at different times. To do what you want you should move the interpolation into the function solver_test2, perhaps something like this:
function dCdt = solver_test2(t,C,t_measured,k1,k2,k12)
k1q = interp1(t_measured, k1, t,'pchip'); % k1 vaues after interpolation
k2q = interp1(t_measured, k2, t,'pchip'); % k2 vaues after interpolation
k12q = interp1(t_measured, k12, t,'pchip'); % k12q vaues after interpolation
dCdt = zeros(2,1);
dCdt(1) = -k1q .* C(1);
dCdt(2) = -k2q .* C(2) + k12q .* C(1);
end
Now you will get the interpolated values of your reaction-rates (?) at any given time for which solver_test2 is called. You only have to adjust the call to ode15s accordingly:
[t,C] = ode15s(@(t,C) solver_test2(t,C,t_measured,k1,k2,k12),tq,initialconcentrations,options)
HTH
  5 comentarios
Bjorn Gustavsson
Bjorn Gustavsson el 17 de Abr. de 2021
Modified. Changed the name of the time that your reaction-rates are interpolated to from tq to t.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by