Please note that the recommended version of Scilab is 2024.1.0. This page might be outdated.
However, this page did not exist in the previous stable version.
dassl
дифференциальное алгебраическое уравнение
Последовательность вызова
[r [,hd]]=dassl(x0,t0,t [,atol,[rtol]],res [,jac] [,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
.- 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 должна быть:
На языке 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 последовательность вызова должна быть следующего вида:
В Fortran она должна быть:
subroutine jac(t,y,yd,pd,cj,rpar,ipar) double precision t, y(*),yd(*),pd(*),cj,rpar(*) integer ipar(*)
- 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
-- индекс времени для которого требуется найти решение.
Описание
Функция dassl
интегрирует дифференциальное алгебраическое уравнение и возвращает изменение y
в заданные моменты времени.
g(t,y,ydot)=0 y(t0)=y0 and ydot(t0)=ydot0
Примеры
function [r, ires]=chemres(t, y, yd) r=[-0.04*y(1)+1d4*y(2)*y(3)-yd(1) 0.04*y(1)-1d4*y(2)*y(3)-3d7*y(2)*y(2)-yd(2) 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 y0=[1;0;0]; yd0=[-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=dassl([y0,yd0],0,t,chemres); info=list([],0,[],[],[],0,0); info(2)=1; y1=dassl([y0,yd0],0,4d10,chemres,info); y2=dassl([y0,yd0],0,4d10,chemres,chemjac,info); //Использование дополнительного аргумента для параметров //----------------------------------- function [r, ires]=chemres(t, y, yd, a, b, c) r=[-a*y(1)+b*y(2)*y(3)-yd(1) a*y(1)-b*y(2)*y(3)-c*y(2)*y(2)-yd(2) y(1)+y(2)+y(3)-1]; ires=0 endfunction function pd=chemjac(x, y, yd, cj, a, b, c) pd=[-a-cj , b*y(3) , b*y(2); a ,-b*y(3)-2*c*y(2)-cj ,-b*y(2); 1 , 1 , 1 ] endfunction y3=dassl([y0,yd0],0,t,list(chemres,0.04,1d4,3d7),list(chemjac,0.04,1d4,3d7)); //использование C-кода //------------ // - создаём C-код rescode=['void chemres(double *t, double y[], double yd[], double r[], int *ires, double rpar[], int ipar[])' ' {' ' r[0] = -0.04*y[0]+1.0e4*y[1]*y[2] -yd[0];' ' r[1] = 0.04*y[0]-1.0e4*y[1]*y[2]-3.0e7*y[1]*y[1]-yd[1];' ' r[2] = y[0]+y[1]+y[2]-1;' ' *ires = 0;' ' }']; jaccode=['void chemjac(double *t, double y[], double yd[], double pd[], double *cj, double rpar[], int ipar[])' ' {' ' /* first column*/' ' pd[0] = -0.04-*cj;' ' pd[1] = 0.04;' ' pd[2] = 1.0;' ' /* second column*/' ' pd[3] = 1.0e4*y[2];' ' pd[4] = -1.0e4*y[2]-2*3.0e7*y[1]-*cj;' ' pd[5] = 1.0;' ' /* third column*/' ' pd[6] = 1.0e4*y[1];' ' pd[7] = -1.0e4*y[1];' ' pd[8] = 1.0;' ' }']; mputl([rescode;jaccode],fullfile(TMPDIR,'mycode.c')) //создаём C-файл в директории TMPDIR // - компилируем ilib_for_link(['chemres','chemjac'],fullfile(TMPDIR,'mycode.c'),[],'c','',fullfile(TMPDIR,'loader.sce')); // - связываем его с Scilab'ом exec(fullfile(TMPDIR,'loader.sce')) // - вызов dassl y4=dassl([y0,yd0],0,t,'chemres','chemjac');
Смотрите также
- ode — программа решения обыкновенных дифференциальных уравнений
- dasrt — программа решения дифференциально-алгебраических уравнений (ДАУ) с пересечением нуля
- impl — дифференциальное алгебраическое уравнение
- fort — Fortran or C user routines call
- link — dynamic linker
- external — объект Scilab'а, внешняя функция или подпрограмма
Report an issue | ||
<< dasrt | Дифференциальное счисление, интегрирование | diff >> |