программа решения дифференциально-алгебраических уравнений (ДАУ) с пересечением нуля
[r,nn,[,hd]] = dasrt(x0,t0,t [,atol,[rtol]],res [,jac],ng, surf [,info] [,hd])
представляет собой либо y0 (ydot0
оценён с помощью dassl с нулём в качестве первой оценки),
либо матрицу [y0 ydot0]. g(t,y0,ydot0)
должна быть равной нулю. Если оценка ydot0 известна,
то установите info(7)=1.
вещественный вектор-столбец исходных условий.
вещественный вектор-столбец производной y по
времени в момент t0 (может быть оценкой).
вещественное число, начальный момент времени.
вещественный скаляр или вектор. Указывает моменты времени для которых
необходимо найти решение. Заметьте, что вы можете получить решение в
каждой точке шага dassl установкой info(2)=1.
вектор с двумя элементами [times num], где
times - значение момента времени, когда поверхность пересекается,
num - номер пересечённой поверхности.
вещественный скаляры или вектор-столбцы того же размера, что и
y. atol, rtol дают, соответственно, допуски абсолютной и относительной ошибки решения. Если это векторы, то допуски определены для каждого элемента y.
внешняя функция, список или строка.
Вычисляет значение g(t,y,ydot). Это может быть:
Функция Scilab.
Её последовательность вызова должна быть
[r,ires]=res(t,y,ydot) и
res должна возвращать остаток
r=g(t,y,ydot) и флаг ошибки
ires. ires = 0, если
res удалось вычислить r;
=-1, если остаток локально не определён для
(t,y,ydot); =-2, если
параметры находятся вне допустимого диапазона.
Список.
Эта форма позволяет передавать функции параметры, отличные от t,
y, ydot. Это выполняется следующим способом:
list(res,x1,x2,...)
где последовательность вызова функции res
теперь имеет вид:
r = res(t,y,ydot,x1,x2,...)
res по-прежнему возвращает
r=g(t,y,ydot как функцию от
(t,y,ydot,x1,x2,...).
![]() | Внимание:
эта форма не должна использоваться, если нет дополнительного аргумента для передачи в функцию. |
Строка.
Она должна ссылаться на имя подпрограммы на языке C или Fortran, связанной с Scilab'ом.
Последовательность вызова на языке C должна быть:
void res(double *t, double y[], double yd[], double r[],
int *ires, double rpar[], int ipar[])
На языке Fortran она должна быть:
subroutine res(t,y,yd,r,ires,rpar,ipar) double precision t, y(*),yd(*),r(*),rpar(*) integer ires,ipar(*)
Массивы rpar и ipar должны быть,
но не должны использоваться.
внешняя функция, список или строка.
Вычисляет значение dg/dy+cj*dg/dydot для заданной
величины параметра cj.
Функция Scilab.
Её последовательность вызова должна быть
r=jac(t,y,ydot,cj) и
jac должна возвращать
r=dg(t,y,ydot)/dy+cj*dg(t,y,ydot)/dydot, где
cj - вещественный скаляр.
Список.
Он должен иметь следующий вид:
list(jac,x1,x2,...)
где последовательность вызова функции
jac теперь имеет вид:
r = jac(t,y,ydot,cj,x1,x2,...)
jac по-прежнему возвращает
dg/dy+cj*dg/dydot в виде функции от
(t,y,ydot,cj,x1,x2,...).
Символьная строка.
Она должна ссылаться на имя подпрограммы на языке C или Fortran, связанной с Scilab.
В C последовательность вызова должна быть следующего вида:
void jac(double *t, double y[], double yd[], double pd[],
double *cj, double rpar[], int ipar[])
В Fortran она должна быть::
subroutine jac(t,y,yd,pd,cj,rpar,ipar) double precision t, y(*),yd(*),pd(*),cj,rpar(*) integer ipar(*)
внешняя функция, список или строка. Вычисляет
значение вектор-столбца surf(t,y) с
ng элементами. Каждый элемент определяет поверхность.
Она может быть определена как:
Функция Scilab.
Её последовательность вызова должна быть
surf(t,y).
Список.
Он должен быть следующего вида:
list(surf,x1,x2,...)
где последовательность вызова функции
surfтеперь имеет вид:
r = surf(t,y,x1,x2,...)
Символьная строка.
Она должна ссылаться на имя процедуры Fortran, связанной с Scilab'ом.
В C последовательность вызова должна быть:
void surf(int *ny, double *t, double y[], int *ng, double gout[])
В Fortran она должна быть:
subroutine surf(ny,t,y,ng,gout) double precision t, y(*),gout(*) integer ny,ng
список, содержащий 7 элементов.
Значение по умолчанию равно list([],0,[],[],[],0,0).
вещественный скаляр, который указывает максимальное время, для которого
g может выполняться, либо пустая матрица
[], если время не ограничено.
флаг, который указывает возвращать ли dassl её
вычисленные промежуточные значения (flag=1) или только
пользователь определяет значения моментов времени
(flag=0).
двухэлементный вектор, который указывает определение [ml,mu]
матрицы пределов, вычисленной с помощью jac;
r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)".
Если jac возвращает полную матрицу, то
info(3)=[].
вещественный скаляр, который указывает максимальный размер шага. Установите
info(4)=[], если ограничений нет.
вещественный скаляр, который указывает исходный размер шага.
Установите info(5)=[], если он не определён.
info(6)=1, если известно, что решение отрицательное, в
противном случае установите info(6)=0.
info(7)=1, если ydot0 является
просто оценкой; info(7)=0, если
g(t0,y0,ydot0)=0.
вещественный вектор, который позволяет хранить контекст
dassl и продолжать интегрирование.
вещественная матрица. Каждый столбец является вектором [t;x(t);xdot(t)],
где t - индекс времени для которого требуется найти решение.
Решение дифференциального уравнения в неявном виде.
g(t,y,ydot) = 0 y(t0) = y0 and ydot(t0) = ydot0
Возвращает моменты времени пересечения поверхности и номер поверхности в nn.
Подробные примеры можно найти в SCIDIR/tests/unit_tests/dassldasrt.tst
// dy/dt = ((2*log(y)+8)/t -5)*y, y(1) = 1, 1<=t<=6 // g1 = ((2*log(y)+8)/t - 5)*y // g2 = log(y) - 2.2491 y0 = 1; t = 2:6; t0 = 1; y0d = 3; atol = 1.d-6; rtol = 0; ng = 2; deff('[delta,ires] = res1(t,y,ydot)', 'ires=0; delta = ydot-((2*log(y)+8)/t-5)*y') deff('rts = gr1(t,y)', 'rts = [((2*log(y)+8)/t-5)*y;log(y)-2.2491]') [yy,nn] = dasrt([y0,y0d], t0, t, atol, rtol, res1, ng, gr1); // (Should return nn = [2.4698972 2]) | ![]() | ![]() |