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

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

Answer 1

[Bjorn Gustavsson]

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

[Quanshun An]

Thanks for your answer; I modified the code as your suggestions, but the ode15s solver didn't work, neither. The command window showed: the function or the 'k1' is not defined.

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)

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,t_measured,k1,k2,k12),tq,initialconcentrations,options)

function dCdt = solver_test2(t,C,t_measured,k1q,k2q,k12q)

k1q = interp1(t_measured, k1, tq,'pchip'); % k1 vaues after interpolation

k2q = interp1(t_measured, k2, tq,'pchip'); % k2 vaues after interpolation

k12q = interp1(t_measured, k12, tq,'pchip'); % k12q vaues after interpolation

dCdt = zeros(2,1);

dCdt(1) = -k1q .* C(1);

dCdt(2) = -k2q .* C(2) + k12q .* C(1);

end

[Bjorn Gustavsson]

My bad, sloppy modification. I forgot to change name of input-parametrers in the function-definition. Is now changed.

[Quanshun An]

Hi, Bjorn. Thanks for your help.I tried this method, but the result showed that the 'tq' is not defined. I also then tried to global 'tq', but it didn't work neither. Do you know how to solve it?

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)

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,t_measured,k1,k2,k12),tq,initialconcentrations,options)

function dCdt = solver_test2(t,C,t_measured,k1,k2,k12)

k1q = interp1(t_measured, k1, tq,'pchip'); % k1 vaues after interpolation

k2q = interp1(t_measured, k2, tq,'pchip'); % k2 vaues after interpolation

k12q = interp1(t_measured, k12, tq,'pchip'); % k12q vaues after interpolation

dCdt = zeros(2,1);

dCdt(1) = -k1q .* C(1);

dCdt(2) = -k2q .* C(2) + k12q .* C(1);

end

[Bjorn Gustavsson]

Modified. Changed the name of the time that your reaction-rates are interpolated to from tq to t.

[Quanshun An]

Dear Bjorn,

Yes, it works, I do appreicate your help!

Best regards