ODE not evaluating correctly with time dependent parameter

Question

Riley Stein 2020-7-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/569785-ode-not-evaluating-correctly-with-time-dependent-parameter

评论： 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更早的评论隐藏 4更早的评论

Star Strider 2020-7-23

在 MATLAB Online 中打开

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.)

.

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 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.

请先登录，再进行评论。

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

Answer 1

Star Strider 2020-7-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/569785-ode-not-evaluating-correctly-with-time-dependent-parameter#answer_470650

在 MATLAB Online 中打开

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更早的评论

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 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!

请先登录，再进行评论。