dae

Аргументы

initial

вектор-столбец. Он может быть равен x0 или [x0;xdot0], где x0 -- это значение состояния в начальный момент времени t0, а xdot0 -- значение производной исходного состояния или его оценка (см. ниже).

t0

вещественное число, исходный момент времени.

t

Вещественный скаляр или вектор. Указывает моменты времени для которых нужно найти решение. Заметьте, что вы можете получить решение в каждой точке шага ДАУ с помощью установки %DAEOPTIONS(2)=1 .

rtol

вещественный скаляр или вектор-столбец того же размера, что и x0, допуск относительной ошибки. Если rtol является вектором, то допуски определяются для каждой составляющей состояния.

atol

вещественный скаляр или вектор-столбец того же размера, что и x0, допуск абсолютной ошибки. Если atol является вектором, то допуски определяются для каждой составляющей состояния.

res

внешняя функция, которая вычисляет значение g(t,y,ydot). Она может быть:

функцией Scilab'а

В этом случае последовательность вызова должна быть [r,ires]=res(t,x,xdot), а res должна возвращать остаток r=g(t,x,xdot) и флаг ошибки ires.

ires = 0, если res смогла вычислить r;

ires = -1, если остаток для g(t,x,xdot) локально не определён;

ires =-2, если параметры находятся вне допустимого диапазона.

списком

Эта форма внешней функции используется для передачи параметров в функцию. Она может иметь следующий вид:

list(res,p1,p2,...)

где последовательность вызова функции res теперь

r=res(t,y,ydot,p1,p2,...)

Функция res по-прежнему возвращает значение остатка в виде функции от (t,x,xdot,x1,x2,...), а p1, p2,... являются параметрами функции.

символьной строкой

она должна ссылаться на имя подпрограммы на языке C или fortran, в предположении, что <r_name> является заданным именем.

• Последовательность вызова Fortran должна быть

  <r_name>(t, x, xdot, res, ires, rpar, ipar)

  t, x(*), xdot(*), res(*), rpar(*) имеют удвоенную точность;

  ires, ipar(*) являются целочисленными.

• Последовательность вызова C должна быть

  C2F(<r_name>)(double *t, double *x, double *xdot, double *res, integer *ires, double *rpar, integer *ipar) ,

где

• t -- текущее значение времени;

• x -- массив состояния;

• xdot -- массив производных состояния;

• res -- массив остатков;

• ires -- индикатор выполнения;

• rpar -- массив целочисленных значений параметров с плавающей запятой, которые нужны, но не могут быть установлены с помощью функции dae.

• ipar -- массив целочисленных значений параметров с плавающей запятой, которые нужны, но не могут быть установлены с помощью функции dae.

jac

Внешняя функция вычисляет значение dg/dx+cj*dg/dxdot для заданного значения параметра cj. Она может быть

функцией Scilab'а

В этом случае последовательность вызова должна быть r=jac(t,x,xdot,cj), а jac должна возвращать r=dg(t,x,xdot)/dy+cj*dg(t,x,xdot)/dxdot, где cj -- вещественный скаляр.

списком

Эта форма внешней функции используется для передачи параметров в функцию. Она может иметь следующий вид:

list(jac,p1,p2,...)

где последовательность вызова функции jac теперь

r=jac(t,x,xdot,p1,p2,...)

Функция jac по-прежнему возвращает dg/dx+cj*dg/dxdot как функцию от (t, x, xdot, cj, p1, p2,...).

символьной строкой

она должна ссылаться на имя подпрограммы на языке C или fortran, в предположении, что <j_name> является заданным именем.

• Последовательность вызова Fortran должна быть

  <j_name>(t, x, xdot, r, cj, ires, rpar, ipar)

  t, x(*), xdot(*), r(*), ci, rpar(*) имеют удвоенную точность;

  ires, ipar(*) являются целочисленными.

• Последовательность вызова C должна быть

  C2F(<j_name>)(double *t, double *x, double *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar) ,

где t, x, xdot, ires, rpar, ipar имеют аналогичное определение, что и выше, r -- массив результатов

surface

Внешняя функция вычисляет значение вектор-столбца surface(t,x) с количеством элементов ng. Каждый элемент определяет поверхность.

функцией Scilab'а

В этом случае последовательность вызова должна быть r=surface(t,x). Эта функция должна вернуть вектор с ng элементами.

списком

Эта форма внешней функции используется для передачи параметров в функцию. Она может иметь следующий вид:

list(surface,p1,p2,...)

где последовательность вызова функции surface теперь имеет вид:

r=surface(t,x,p1,p2,...)

символьной строкой

она должна ссылаться на имя подпрограммы на языке C или fortran, в предположении, что <s_name> является заданным именем.

• Последовательность вызова Fortran должна быть

  <j_name>(t, x, xdot, r, cj, ires, rpar, ipar)

  t, x(*), xdot(*), r(*), ci, rpar(*) имеют удвоенную точность;

  ires, ipar(*) являются целочисленными.

• Последовательность вызова C должна быть

  C2F(<j_name>)(double *t, double *x, double *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar) ,

где t, x, xdot, ires, rpar, ipar имеют аналогичное определение, что и выше, ng -- количество поверхностей, nx -- размерность состояния и r -- массив результатов.

rd

вектор с двумя элементами [times num], где times -- значение момента времени пересечения поверхности, num -- число пересечённых поверхностей

hd

вещественный вектор, в качестве аргумента на выходе он хранит контекст dae. Он может быть использован в качестве входного аргумента для возобновления интегрирования (горячий перезапуск).

y

вещественная матрица. Если %DAEOPTIONS(2)=1 ,то каждый столбец является вектором вида [t;x(t);xdot(t)], где t -- индекс времени для которого вычислено решение. В противном случае y -- вектор вида [x(t);xdot(t)].

Описание

Функция dae является шлюзом, построенным над функциями dassl и dasrt, разработанными для явного интегрирования дифференциальных уравнений.

g(t,x,xdot)=0
x(t0)=x0  and   xdot(t0)=xdot0

Если xdot0 не указан в исходном аргументе, то функция dae пытается вычислить его решая уравнение g(t,x0,xdot0)=0/

Если xdot0 указан в исходном аргументе, то он может быть либо совместимой производной (compatible derivative), удовлетворяющей условию g(t,x0,xdot0)=0, либо приближённым значением. В последнем случае %DAEOPTIONS(7) должен быть установлен в 1.

Конкретные примеры использования внешних функций, написанных на языке Scilab и C, представлены в modules/differential_equations/tests/unit_tests/dassldasrt.tst

Примеры

//Пример с кодом Scilab
function [r, ires]=chemres(t, y, yd)
r(1) = -0.04*y(1) + 1d4*y(2)*y(3) - yd(1);
r(2) =  0.04*y(1) - 1d4*y(2)*y(3) - 3d7*y(2)*y(2) - yd(2);
r(3) =       y(1) +     y(2)      + y(3)-1;
ires =  0;
endfunction
function pd=chemjac(x, y, yd, cj)
pd=[-0.04-cj , 1d4*y(3)               , 1d4*y(2);
0.04    ,-1d4*y(3)-2*3d7*y(2)-cj ,-1d4*y(2);
1       , 1                      , 1       ]
endfunction

x0=[1; 0; 0];
xd0=[-0.04; 0.04; 0];
t=[1.d-5:0.02:.4, 0.41:.1:4, 40, 400, 4000, 40000, 4d5, 4d6, 4d7, 4d8, 4d9, 4d10];

y=dae([x0,xd0],0,t,chemres);// возвращает запрошенные моменты времени наблюдения

%DAEOPTIONS=list([],1,[],[],[],0,0); // просит вернуть сетку точек dae
y=dae([x0,xd0],0,4d10,chemres); // без якобиана
y=dae([x0,xd0],0,4d10,chemres,chemjac); // с якобианом

//пример с кодом C (необходим C-компилятор) --------------------------------------------------
//-1- создать C-код в TMPDIR - Уравнение Вандерпола, неявная форма
code=['#include <math.h>'
'void res22(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)'
'{res[0] = yd[0] - y[1];'
' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
' '
'void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)'
'{pd[0]=*cj - 0.0;'
' pd[1]=    - (-200.0*y[0]*y[1] - 1.0);'
' pd[2]=    - 1.0;'
' pd[3]=*cj - (100.0*(1.0 - y[0]*y[0]));}'
' '
'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
'{ groot[0] = y[0];}']
mputl(code,TMPDIR+'/t22.c') 

//-2- скомпилировать и загрузить его
ilib_for_link(['res22' 'jac22' 'gr22'],'t22.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/t22loader.sce');
exec(TMPDIR+'/t22loader.sce')

//-3- запустить
rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;

//простое моделирование
t=0:0.003:300;
yy=dae([y0,y0d],t0,t,atol,rtol,'res22','jac22');
clf();plot(yy(1,:),yy(2,:))

// найти первую точку, где yy(1)=0
[yy,nn,hotd]=dae("root",[y0,y0d],t0,300,atol,rtol,'res22','jac22',ng,'gr22');
plot(yy(1,1),yy(2,1),'r+')
xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)))

// горячий перезапуск для следующей точки
t01=nn(1);[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
[yy,nn,hotd]=dae("root",[y01,y0d1],t01,300,atol,rtol,'res22','jac22',ng,'gr22',hotd);
plot(yy(1,1),yy(2,1),'r+')
xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)))

Смотрите также

• ode — программа решения обыкновенных дифференциальных уравнений
• daeoptions — установка опций для dae (программа решения дифференциальных алгоритмических уравнений)
• dassl — дифференциальное алгебраическое уравнение
• impl — дифференциальное алгебраическое уравнение
• fort — Fortran or C user routines call
• link — dynamic linker
• external — объект Scilab'а, внешняя функция или подпрограмма