dasrt
программа решения дифференциально-алгебраических уравнений (ДАУ) с пересечением нуля
Синтаксис
[r,nn,[,hd]] = dasrt(x0,t0,t [,atol,[rtol]],res [,jac],ng, surf [,info] [,hd])
Аргументы
- x0
представляет собой либо
y0
(ydot0
оценён с помощьюdassl
с нулём в качестве первой оценки), либо матрицу[y0 ydot0]
.g(t,y0,ydot0)
должна быть равной нулю. Если оценкаydot0
известна, то установитеinfo(7)=1
.- y0
вещественный вектор-столбец исходных условий.
- ydot0
вещественный вектор-столбец производной
y
по времени в моментt0
(может быть оценкой).
- t0
вещественное число, начальный момент времени.
- t
вещественный скаляр или вектор. Указывает моменты времени для которых необходимо найти решение. Заметьте, что вы можете получить решение в каждой точке шага
dassl
установкойinfo(2)=1
.- nn
вектор с двумя элементами
[times num]
, гдеtimes
- значение момента времени, когда поверхность пересекается,num
- номер пересечённой поверхности.- atol, rtol
вещественный скаляры или вектор-столбцы того же размера, что и
y
.atol, rtol
дают, соответственно, допуски абсолютной и относительной ошибки решения. Если это векторы, то допуски определены для каждого элементаy
.- res
внешняя функция, список или строка. Вычисляет значение
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
должны быть, но не должны использоваться.
- jac
внешняя функция, список или строка. Вычисляет значение
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
внешняя функция, список или строка. Вычисляет значение вектор-столбца
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
- info
список, содержащий
7
элементов. Значение по умолчанию равноlist([],0,[],[],[],0,0)
.- info(1)
вещественный скаляр, который указывает максимальное время, для которого
g
может выполняться, либо пустая матрица[]
, если время не ограничено.- info(2)
флаг, который указывает возвращать ли
dassl
её вычисленные промежуточные значения (flag=1
) или только пользователь определяет значения моментов времени (flag=0
).- info(3)
двухэлементный вектор, который указывает определение
[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(4)=[]
, если ограничений нет.- info(5)
вещественный скаляр, который указывает исходный размер шага. Установите
info(5)=[]
, если он не определён.- info(6)
info(6)=1
, если известно, что решение отрицательное, в противном случае установитеinfo(6)=0
.- info(7)
info(7)=1
, еслиydot0
является просто оценкой;info(7)=0
, еслиg(t0,y0,ydot0)=0
.
- hd
вещественный вектор, который позволяет хранить контекст
dassl
и продолжать интегрирование.- r
вещественная матрица. Каждый столбец является вектором
[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])
Смотрите также
Report an issue | ||
<< daskr | Дифференциальное счисление, интегрирование | dassl >> |