dae

Parâmetros

initial

um vetor coluna. Pode ser igual a x0or [x0;xdot0], onde x0 é o valor de estado no tempo inicial t0 e ydot0 é o valor da derivada do estado inicial ou uma estimativa dela (ver abaixo).

t0

número real, o tempo inicial

t

escalar real ou vetor. Fornece instantes para os quais você deseja uma solução. Note que você pode obter soluções para cada ponto de passo de dae fazendo %DAEOPTIONS(2)=1.

rtol

escalar real ou vetor coluna com o mesmo tamanho que x0. A tolerância do erro relativo da solução. Se rtol for um vetor, as tolerâncias são especificadas para cada componente do estado.

atol

escalar real ou vetor coluna com o mesmo tamanho que x0. A tolerância do erro absoluto da solução. Se atol for um vetor, as tolerâncias são especificadas para cada componente do estado.

res

uma função externa (external). Computa o valor de g(t,y,ydot). Pode ser

Uma função do Scilab

Neste caso, a sua seqüência de chamamento pode ser [r,ires]=res(t,x,xdot) e res deve retornar o resíduo r=g(t,x,xdot) e o indicador de erro ires. ires = 0 se res obtiver sucesso ao computar r, =-1 se o resíduo é indefinido localmente para (t,x,xdot), =-2 se os parâmetros estão fora do intervalo admissível.

Uma lista

Esta forma é utilizada para passar parâmetros à função. Deve ser como segue:

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

Onde a seqüência de chamamento da função res é agora

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

res ainda retorna o valor residual como uma função de (t,x,xdot,x1,x2,...), e p1,p2,... são parâmetros da função.

Um string

Deve se referir ao nome subrotina C ou Fortran. Supondo que <r_name> ié o nome dado.

• A seqüência de chamamento em Fortran deve ser

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

  double precision t,x(*),xdot(*),res(*),rpar(*)

  integer ires,ipar(*)

• A seqüência de chamamento em C deve ser

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

onde

• t é o valor de tempo corrente

• x é o array de estados

• xdot é o array das derivadas dos estados

• res é o array de resíduos

• ires é o indicador de execução

• rpar é o array de valores de parâmetros em ponto flutuante, necessário, mas não pode ser definido pela função dae

• ipar é o array de valores de parâmetros inteiros, necessário, mas não pode ser definido pela função dae

jac

uma função externa (external). Computa o valor de dg/dx+cj*dg/dxdot para um dado valor do parâmetro cj. Pode ser

Uma função do Scilab

Sua seqüência de chamamento deve ser r=jac(t,x,xdot,cj) e a função jac deve retornar r=dg(t,x,xdot)/dy+cj*dg(t,x,xdot)/dxdot onde cj é um escalar real.

Uma lista

Esta forma é utilizada para passar parâmetros à função. Deve ser como segue:

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

Onde a seqüência de chamamento da função jac é agora

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

jac ainda retorna dg/dx+cj*dg/dxdot como uma função de (t,x,xdot,cj,p1,p2,...).

Um string

Deve se referir ao nome de uma subrotina C ou Fortran. Supondo que <j_name> é o nome dado,

• A seqüência de chamamento em Fortran deve ser

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

  double precision t, x(*), xdot(*), r(*), ci, rpar(*)

  integer ires, ipar(*)

• A seqüência de chamamento dem C deve ser

  C2F(<j_name>)(double *t, double *x, double *xdot, double *r, double *cj,

  integer *ires, double *rpar, integer *ipar)

onde t, x, xdot, ires, rpar, ipar são definidas semelhantemente como acima, r é o array de resultados.

surface

uma função externa (external). Computa o valor de cada vetor coluna surface(t,x) como componentes ng. Cada componente define uma superfície.

Uma função do Scilab

Sua seqüência de chamamento deve ser r=surface(t,x), esta função deve retornar um vetor com ng elementos.

Uma lista

Esta forma é utilizada para passar parâmetros à função. Deve ser como segue:

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

Onde a seqüência de chamamento da função surface é agora

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

String

Deve se referir ao nome de uma rotina C ou Fortran. Supondo que <s_name> é o nom dado,

• A seqüência de chamamento em Fortran deve ser

  <r_name>(nx, t, x, ng, r, rpar, ipar)

  double precision t, x(*), r(*), rpar(*)

  integer nx, ng,ipar(*)

• A seqüência de chamamento em C deve ser

  C2F(<r_name>)(double *t, double *x, double *xdot, double *r, double *cj,

  integer *ires, double *rpar, integer *ipar)

onde t, x, rpar, ipar são definidas semelhantemente como acima, ngé o número de superfícies, nx é a dimensão do estado e r é o array de resultados.

rd

um vetor com duas entradas [times num] times é o valor do tempo no qual a superfície é cruzada, num é o número da superfície cruzada.

hd

um vetor de reais que permite armazenar o contexto de dae. Pode ser utilizado como argumento de entrada para retomar a integração (recomeço rápido).

y

matriz de reais. Se %DAEOPTIONS(2)=1, cada coluna é o vetor [t;x(t);xdot(t)] onde t é o índice do tempo para o qual a solução foi computada. De outro modo, y é o vetor [x(t);xdot(t)].

Descrição

A função dae é uma porta construída sobre as funções dassl e dasrt designada para equações diferenciais implícitas.

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

Se xdot0 não for fornecido no argumentoinitial, a função dae tenta computá-lo resolvendo g(t,x0,xdot0)=0,

Se xdot0 for fornecido no argumentoinitial, pode ser tanto uma derivada compatível satisfazendo g(t,x0,xdot0)=0 ou um valor aproximado . No último caso %DAEOPTIONS(7) deve ser ajustado para 1.

Exemplos detalhados utilizando funções externas codificadas em C e Fortran são dados em modules/differential_equations/tests/unit_tests/dassldasrt.tst

Exemplos

//Exemplo com um código 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);// retorna os pontos de observação requisitados

%DAEOPTIONS=list([],1,[],[],[],0,0); // pede que a malha de dae seja retornada
y=dae([x0,xd0],0,4d10,chemres); // sem jacobiano
y=dae([x0,xd0],0,4d10,chemres,chemjac); // com jacobiano

//exemplo com um código C (requer-se um compilador C) --------------------------------------------------
//-1- criando os códigos C em TMPDIR - equação de Vanderpol, forma implícita
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- compilando e carregando
ilib_for_link(['res22' 'jac22' 'gr22'],'t22.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/t22loader.sce');
exec(TMPDIR+'/t22loader.sce')

//-3- executando
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;

//simulação simples
t=0:0.003:300;
yy=dae([y0,y0d],t0,t,atol,rtol,'res22','jac22');
clf();plot(yy(1,:),yy(2,:))

//achando o primeiro ponto onde 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)))

//recomeço rápido para o próximo ponto
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)))

Ver Também

• ode
• daeoptions
• dassl
• impl
• fort
• link
• external