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

2 views (last 30 days)
Quanshun An on 16 Apr 2021
Commented: Quanshun An on 17 Apr 2021
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

Bjorn Gustavsson on 16 Apr 2021
Edited: Bjorn Gustavsson on 17 Apr 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
Quanshun An on 17 Apr 2021
Dear Bjorn,
Yes, it works, I do appreicate your help!
Best regards