Solving ODE in Scilab

Question

I am trying to do some electric circuit analysis in Scilab by solving an ODE. But I need to change an ODE depending on current value of the function. I have implemented the solution in Scala using RK4 method and it works perfectly. Now I am trying to do the same but using standard functions in Scilab. And it is not working. I have tried to solve these two ODEs seperately and it is OK.

clear
state = 0 // state 0 is charging, 1 is discharging
vb = 300.0; vt = 500.0; 
r = 100.0; rd = 10.0;
vcc = 600;
c = 48.0e-6;

function dudx = curfunc(t, uu)
    if uu < vb then state = 0
    elseif uu > vt state = 1
    end

    select state
    case 0 then // charging
        dudx = (vcc - uu) / (r * c)
    case 1 then // discharging
        dudx = - uu / (rd * c) + (vcc - uu) / (r * c)
    end
endfunction
y0 = 0
t0 = 0
t = 0:1e-6:10e-3

%ODEOPTIONS=[1, 0, 0, 1e-6, 1e-12, 2, 500, 12, 5, 0, -1, -1]
y = ode(y0, t0, t, 1e-3, 1e-6, curfunc)
clear %ODEOPTIONS
plot(t, y)

so here I am solving for a node voltage, if node voltage is exceeding top threshold (vt) then discharging ODE is used, if node voltage is going below the bottom voltage (vb) then charging ODE is used. I have tried to play with %ODEOPTIONS but no luck

Answer 1

I think it boils down to a standard problem: The standard ODE solvers have step size control. In contrast to the fixed step size of RK4. For step size control to work, the ODE function needs to have differentiability at least of the same order as the order of the method. Jumps as in your function are thus extremely unfortunate.

Another point to consider is that the internal steps of the method are not always in increasing time, they may jump backwards. See the coefficient table for Dormant-Price for an example. Thus event based model changes as in your problem may lead to strange effects if the first problem were circumvented.

---

Workaround (not really a solution)

Declare state as a global variable, since as a local variable the value gets always reset to the global constant 0.

function dudx = curfunc(t, uu)
    global state
    ...
endfunction
...
y = ode("fix",y0, t0, t, rtol, atol, curfunc)

Answer 2

you can also use the ode("root"...) option.

The code will ressemble to clear state = 0 // state 0 is charging, 1 is discharging vb = 300.0; vt = 500.0; r = 100.0; rd = 10.0; vcc = 600; c = 48.0e-6;

function dudx = charging(t, uu)
//uu<vt
  dudx = (vcc - uu) / (r * c)
endfunction
function e=chargingssurf(t,uu)
  e=uu-vt
endfunction

function dudx = discharging(t, uu)
//uu<vb
  dudx = - uu / (rd * c) + (vcc - uu) / (r * c)
endfunction
function e=dischargingssurf(t,uu)
  e=uu-vb
endfunction

y0 = 0
t0 = 0
t = 0:1e-6:10e-3

Y=[];
T=[];
[y,rt] = ode("root",y0, t0, t, 1e-3, 1e-6, charging,1,chargingssurf);
disp(rt)
k=find(t(1:$-1)<rt(1)&t(2:$)>=rt(1))
Y=[Y y];;
T=[T t(1:k) rt(1)];
[y,rt] = ode("root",y($), rd(1), t(k+1:$), 1e-3, 1e-6, discharging,1,dischargingssurf);
Y=[Y y];
T=[T t(k+1:$)];

plot(T, Y)

The discharging code seems to be wrong...