ode23s solver: f(x)<0 even though f'(x)>0 as soon as time becomes internal variable

Hello,

I got a problem during my simulation using the ode23s solver, probably because of a very small inconstancy.

Firstly I got:

xInit = [0,0,0,0,0,0,0];
dxdt = odefun(t,x); % with dxdt = [..,..  ,duds, ...,..];
[t,x] = ode23s(odefun, [0 250e-9], xInit); % so i simulate from t = 0ns until t = 250ns

in the odefun there is a statement calculating duds:

ids = ia + ib + ir + ip;
duds = 1/C*ids;

ir is a function of uds and duds is the derivative of this uds:

if uds<0
    ir = f*exp(g*(-uds))-1;
else
    ir = 0;
end

Furthermore:

For 0<t<t0 the following statement applies without exception: iq = 0

As long as this applies the simulation runs smoothly, but at the time instance t0 there is an event triggering the start of the calculation of iq with:

iq = a*tanh(b*(t-t0))^c; % a,b,c are just coefficients and t is the simulation-time

The very moment this calculation starts, there is an error accumulating rapidly. Even though it applies that duds > 0 before and after the time instance t0 the variable uds starts to become more and more negative after t0, although it was becoming more positive before t0 when duds > 0 applies. Because ir is an exp-function of (-uds) it starts to diverge when uds becomes more negative.

I already understood that odexx is not solely taking duds to calculate uds, but it's doing further semi-random steps in order to compute higher-order derivatives. And I guess that suddenly introducing a time dependent variable does not come easily, but how do I solve this.

Discontinuities are fatal, but because tanh(t-t0) = 0 when t=t0 and because tanh is a continuous differentiable as like as its derivatives it should work, right?

Is there a way to make ode work with that, or just tell it to work with duds alone to calculate my uds of the next time instance?

Best regards,

Fabian

[Harald]

Hi,

I have a hard time following this. Could you post a single example that we can run that exposes the issue?

Best wishes,

Harald

[FJ]

Hi,

unfortunately it is hard to supply an easy example but i could share part of the code although i got to encrypt it partially. Everything of interest should be visible, though.

Main.p runs the simulation and gives a plot.Focus on the orange plot of "uds" in the upper subplot.

If you run everything as it is, it works fine. If you add in the "- irbip;" in line 18 of 'calculateDevice.m' the above discribed behaviour is to be seen.

Generally irbip= 0 until 124 ns where a calculation of it is triggered. After this point it applies that:

irbip = a*tanh(b*(t-t0))^c; % a,b,c are just coefficients and t is the simulation-time

I am aware that this piece-wise function may be a continous function but not a continuously differentiable function. So i made its numeric derivative continous, but it didn't help a thing.

I would really appreciate some help.

[Harald]

Running Main.p for me results in

Array indices must be positive integers or logical values.

Error in Main

Also, I doubt that I will be able to sufficiently investigate if I cannot inspect the call to ode23s and the function supplied to it.

I am not sure if I understand the problem at all. If a discontinuity at time 124ns constitutes the problem, you can try solving the ode for [0, 124e-9] and use the final state as an initial state for solving for [124e-9, 250e-9].

Best wishes,

Harald

[FJ]

Hello,

i deeply apologize for this mistake. The new attachment is executable. Double-checked it.

I really appreciate your help.

To give some further information:

1. There is no discontinuity any more.
2. There are varous class objects of the class device (here we call the specific device 'deevice.p') and before calling the solver there is no option to know when this event, that the variable 'id' hits zero, occurs. 'id' crossing zero triggers the start of the calculation of 'irbip'. Before that point of time irbip=0 and starting from there it follows a tanh^x ( with 1<x<3) function. So at point t0, when it gets triggert irbip = 0 applies. Because i don't know when the event happens i don't know how to split it. t0 depends on the device and many more variables. Here it happens to be around 124 ns but under other circumstances it could as well be somewhere else. You can set a conditional breakbpoint on irbip~=0 to find the first simulation point of interest.

Best regards

[Torsten]

Use the "event" option of the ODE solvers:

https://uk.mathworks.com/help/matlab/math/ode-event-location.html

[FJ]

Thanks @Torsten,

now i know how to split it and while doing it i found my mistake.

Unfortunately, splitting it didn't solve the issue.

Maybe someone else finds the following helpful:

It took me some while to find this easy bug. But ode seems to jump back and forth in time. At least in my case. This rarely happened during execution but led to problems and inconsistencies.

[Harald]

This "jumping back and forth" is expected behavior of an ode solver with step size control: if the estimated error violates the tolerances, a smaller step has to be taken. This implies evaluating the function at previous times.