ODE not evaluating correctly with time dependent parameter

2020 7 23
1 个回答
回答已采纳
更新时间：2020 7 27
1 次查看（30 天）

0 个投票

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.

6 个评论
显示 4更早的评论 隐藏 4更早的评论

Star Strider
Star Strider 2020-7-23
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
Steven Lord 2020-7-23
Can you show us the mathematical form of the ODE you're trying to solve?
Riley Stein
Riley Stein 2020-7-23
@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
Star Strider 2020-7-23
Riley Stein — Your approach appears to be correct. Your implementation of the actual differential equation may not be.
Alan Stevens
Alan Stevens 2020-7-24
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
Riley Stein 2020-7-24
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.

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Star Strider
Star Strider 2020-7-24

0 个投票

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.

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Riley Stein
Riley Stein 2020-7-24
I am running into the problem that kredfcn is an unrecognized function or variable in the function dydt.
Star Strider
Star Strider 2020-7-24
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
Riley Stein 2020-7-25
Thank you this worked perfectly. As a curiosity question, how would this work if there were a system of differential equations?
Star Strider
Star Strider 2020-7-25
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
Riley Stein 2020-7-27
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!

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by