Решение ODE в Scilab

Я пытаюсь провести анализ электрических цепей в Scilab, решив ODE. Но мне нужно изменить ODE в зависимости от текущего значения функции. Я реализовал решение в Scala, используя метод RK4, и оно отлично работает. Теперь я пытаюсь сделать то же самое, но с использованием стандартных функций в Scilab. И это не работает. Я попытался решить эти два ODE по отдельности, и это нормально.

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)

Итак, здесь я решаю для напряжения узла, если напряжение узла превышает верхний порог (vt), тогда используется разряд ODE, если напряжение узла ниже нижнего напряжения (vb), тогда используется заряд ODE. Я пробовал играть с %ODEOPTIONS, но безуспешно


scilab ode
person Artem Yastrebov    schedule 17.02.2016    source источник
comment
Не могли бы вы задокументировать проблемные случаи, скрипт работает отлично.   -  person Lutz Lehmann    schedule 17.02.2016

Ответы (2)


arrow_upward
1
arrow_downward

Я думаю, что это сводится к стандартной проблеме: стандартные решатели ODE имеют контроль размера шага. В отличие от фиксированного размера шага RK4. Чтобы управление размером шага работало, функция ODE должна иметь дифференцируемость по крайней мере того же порядка, что и порядок метода. Таким образом, прыжки, как в вашей функции, крайне неудачны.

Еще один момент, который следует учитывать, заключается в том, что внутренние шаги метода не всегда увеличиваются во времени, они могут перемещаться назад. См. Пример в таблице коэффициентов для Dormant-Price. Таким образом, изменения модели на основе событий, как и в вашей проблеме, могут привести к странным эффектам, если первая проблема была обойдена.

Обходной путь (не совсем решение)

Объявите state как глобальную переменную, поскольку в качестве локальной переменной значение всегда сбрасывается до глобальной константы 0.

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

arrow_upward
1
arrow_downward

вы также можете использовать опцию ode ("root" ...).

Код будет иметь вид clear state = 0 // состояние 0 заряжается, 1 разряжается 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)

Код разрядки кажется неправильным ...

person user5694329    schedule 18.02.2016