ODE not evaluating correctly with time dependent parameter

I am trying to pass a time dependent constant to my ode solver however it is not yielding the right values.

Eo = 0.4;
Ei = 0;
v = -0.1;
a = 0.5;
F = 96485.33289;
R = 8.314;
T = 293.15;
tfinal = -1.5/v;
tstep = 0.1;
tspan = 0:tstep:tfinal;
En = Ei + v * tspan;
kmax = 225000;
k0 = kmax * exp(-8);

kred = k0 .* exp(((-a * F)/(R * T)) * (En - Eo));
kox  = k0 .* exp((((1-a) * F)/(R * T)) * (En - Eo));

IC = 1;

[A, B] = ode15s(@(t, y) myODE(t, y, kred, kox, En), tspan, IC);

function dydt = myODE(t, y, kred, kox, En)


kred = interp1(En, kred, t);
kox  = interp1(En, kox, t);

dydt = -kred .* y + kox;

end

A and B are being solved as singular values (0 and 1 respectively). The arrays of kred and kox are yielding the right values. What is confusing me is that I have a similar test code, which returns an array for the values A and B.

v      = 3;
tfinal = 15/v;
tStep  = 1;
tspan  = 0:tStep:tfinal;
x      = tspan * v;
xi     = 0;
xn     = xi + v*tspan;

s      =  0.5 * exp(xn - 3);
u      = -0.5 * exp(xn - 3);

IC = 1;

[A, B] = ode15s(@(t, y) myODE(t, y, s, u, xn), tspan, IC);

function dydt = myODE(t, y, s, u, xn)

s = interp1(xn, s, t);
u = interp1(xn, u, t);

dydt = -s .*y + u;

end

For this code, A and B are returning the correct values and I don't see any striking differences in how the function set up is besides the actual equations of the parameters. Any help is appreciated.

[Star Strider]

The error I get when I run you code is:

Warning: Failure at t=0.000000e+00.  Unable to meet integration tolerances without reducing the
step size below the smallest value allowed (7.905050e-323) at time t. 

Part of the problem is with the interpolations, because interpolating (actually extrapolating) outside the limits of the first argument will produce a NaN value.

However even with:

kred = interp1(En, kred, t, 'linear','extrap');
kox  = interp1(En, kox, t, 'linear','extrap');

(that would allow that extrapolation) the ‘myODE’ function quickly rails. The last finite values for ‘t’, ‘y’, ‘kred’, ‘kox’, and ‘dydt’ are:

[6.334195697423101e-02    2.153919700486809e+303    -8.007668785978458e+04     5.877317566322157e-02    1.724787555309229e+308]

The ode15s function then exits, with:

Warning: Failure at t=6.334246e-02.  Unable to meet integration tolerances without reducing the
step size below the smallest value allowed (2.220446e-16) at time t.  

I have no idea what you are doing, so I cannot advance a solution. (That is also the reason I am not posting this as an Answer, bercause it is not one.)

.

[Steven Lord]

Can you show us the mathematical form of the ODE you're trying to solve?

[Riley Stein]

@starstrider Thank you for the clarification. I am trying to model Butler-Volmer kinetics in Matlab and the differential equations require the constants kred and kox, which change with time. I recently read that time dependent constants can be passed into ODEs so that was my motivation for trying this.

[Star Strider]

Riley Stein — Your approach appears to be correct. Your implementation of the actual differential equation may not be.

[Alan Stevens]

Your En values are all negative, but you are passing positive values of t to the interp function. Are you sure En should be negative.

The corresponding terms, xn, in your other code are positive.

[Riley Stein]

Increasing the tolerance of the ode15s produced results, albeit not the right ones. Some values are being solved as negative, when they should never be negative because there is no physical meaning for it. As for the differential equations, they are from a paper and are right.

Answer 1

[Star Strider]

Create ‘kred’ and ‘kox’ as anonymous functions, and it runs:

Eo = 0.4;
Ei = 0;
v = -0.1;
a = 0.5;
F = 96485.33289;
R = 8.314;
T = 293.15;
tfinal = -1.5/v;
tstep = 0.1;
tspan = 0:tstep:tfinal;
En = Ei + v * tspan;
kmax = 225000;
k0 = kmax * exp(-8);

Enfcn = @(t) Ei + v * t;
kredfcn = @(t) k0 .* exp(((-a * F)/(R * T)) * (Enfcn(t) - Eo));
koxfcn  = @(t) k0 .* exp((((1-a) * F)/(R * T)) * (Enfcn(t) - Eo));

IC = 1;

[A, B] = ode15s(@(t, y) myODE(t, y), tspan, IC);

function dydt = myODE(t, y)


kred = kredfcn(t);
kox  = koxfcn(t);

dydt = -kred .* y + kox;

end

figure
plot(A,B)
grid

You must determine if it produces the correct result.

[Riley Stein]

I am running into the problem that kredfcn is an unrecognized function or variable in the function dydt.

[Star Strider]

Try this version:

Eo = 0.4;
Ei = 0;
v = -0.1;
a = 0.5;
F = 96485.33289;
R = 8.314;
T = 293.15;
tfinal = -1.5/v;
tstep = 0.1;
tspan = 0:tstep:tfinal;
En = Ei + v * tspan;
kmax = 225000;
k0 = kmax * exp(-8);

Enfcn = @(t) Ei + v * t;
kredfcn = @(t) k0 .* exp(((-a * F)/(R * T)) * (Enfcn(t) - Eo));
koxfcn  = @(t) k0 .* exp((((1-a) * F)/(R * T)) * (Enfcn(t) - Eo));

myODE = @(t, y) -kredfcn(t) .* y + koxfcn(t);

IC = 1;

[A, B] = ode15s(myODE, tspan, IC);

figure
plot(A,B)
grid

The ‘myODE’ function is now an anonymous function as well, so you a can run it in any script.

.

[Riley Stein]

Thank you this worked perfectly. As a curiosity question, how would this work if there were a system of differential equations?

[Star Strider]

As always, my pleasure!

(My apologies for not seeing tthe obvious relationship between ‘t’ and the functions earllier.)

It should work the same way for a system of differential equations. (It obviously depends on the differential equations and the functions, if they are not the same as they are here.)

[Riley Stein]

For future reference, the first approach also works. So long as the kred and kox parameters are passed to the function (t, y, kredfcn, koxfcn) and then called as koxfcn(t) and kredfcn(t) in the function, the script will run!