ODE not evaluating correctly with time dependent parameter

1 次查看（过去 30 天）
Riley Stein
Riley Stein 2020-7-23
评论： Riley Stein 2020-7-27
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更早的评论
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
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更早的评论
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by