Scilab Home page | Wiki | Bug tracker | Forge | Mailing list archives | ATOMS | File exchange
Please login or create an account
Change language to: English - Français - 日本語 - Русский

Please note that the recommended version of Scilab is 6.0.1. This page might be outdated.
See the recommended documentation of this function

Ajuda Scilab >> Equações Diferenciais > dassl

dassl

Equação diferencial algébrica

Seqüência de Chamamento

[r [,hd]]=dassl(x0,t0,t [,atol,[rtol]],res [,jac] [,info] [,hd])

Parâmetros

x0

pode ser tanto y0 (ydot0 é estimado por dassl com a primeira estimativa sendo 0) ou a matriz [y0 ydot0]. g(t,y0,ydot0) deve ser igual a zero. Se você conhece apenas uma estimativa de ydot0 faça info(7)=1

y0

vetor coluna de reais de condições iniciais

ydot0

vetor coluna de reais da derivada do tempo de y em t0 (pode ser uma estimativa).

t0

número real, é o instante 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 dassl fazendo info(2)=1.

atol, rtol

escalares reais ou vetores colunas com o mesmo tamanho que y. atol, rtol fornecem respectivamente as tolerâncias de erros absolutos e relativos da solução. Se forem vetores, as tolerâncias são especificadas para cada elemento de y.

res

função, lista ou string externos. Computa o valor de g(t,y,ydot). Pode ser :

  • Uma função do Scilab

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

  • Uma lista

    Esta forma permite passar outros parâmetros além de t,y,ydot à função. Deve ser como segue:

    list(res,x1,x2,...)

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

    r=res(t,y,ydot,x1,x2,...)

    res ainda retorna r=g(t,y,ydot) como função de (t,y,ydot,x1,x2,...).

  • Um string

    Deve se referir ao nome subrotina C ou Fortran ligada ao Scilab.

    Em C a seqüência de chamamento:

    void res(double *t, double y[], double yd[], double r[],
             int *ires, double rpar[], int ipar[])

    Em Fortran deve ser:

    subroutine res(t,y,yd,r,ires,rpar,ipar)
    double precision t, y(*),yd(*),r(*),rpar(*)
    integer ires,ipar(*)

    Os arrays rpar e ipar devem estar presentes, mas não podem ser utilizados.

jac

função, lista ou string externos. Computa o valor de dg/dy+cj*dg/dydot para um dado valor do parâmetro cj

  • Uma função do Scilab

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

  • Uma lista

    Deve ser como segue

    list(jac,x1,x2,...)

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

    r=jac(t,y,ydot,cj,x1,x2,...)

    jac ainda retorna dg/dy+cj*dg/dydot como uma função de (t,y,ydot,cj,x1,x2,...).

  • Um string

    Deve se referir ao nome de uma subrotina C ou Fortran ligada ao Scilab

    Em C, a seqüência de chamamento deve ser

    void jac(double *t, double y[], double yd[], 
             double pd[], double *cj, double rpar[], int ipar[])

    Em Fortran deve ser

    subroutine jac(t,y,yd,pd,cj,rpar,ipar)
    double precision t, y(*),yd(*),pd(*),cj,rpar(*)
    integer ipar(*)
info

lista opcional que contém 7 elementos. O valor padrão é list([],0,[],[],[],0,0);

info(1)

escalar real que fornece o tempo máximo para o qual se permite avaliar g , ou uma matriz vazia [] se não houver limites impostos sobre o tempo.

info(2)

indica se dassl retorna seus valores intermediários computados (flag=1) ou apenas os valores de pontos tempo especificados pelo usuário (flag=0).

info(3)

vetor de dois componentes que fornece a definição [ml,mu] de matriz de bandas computada por jac; r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)" .Se jac retorna uma matriz cheia, faça info(3)=[].

info(4)

escalar real que fornece o tamanho de passo máximo. Faça info(4)=[] se não houver limitação.

info(5)

escalar real que fornece o passo inicial. Faça info(5)=[] se não for especificado.

info(6)

faça info(6)=1 caso se saiba que a soulução não é negativa, em outro caso, ajuste info(6)=0.

info(7)

faça info(7)=1 se ydot0 é apenas uma estimativa, info(7)=0 se g(t0,y0,ydot0)=0.

hd

vetor de reais que permite armazenar o contexto de dassl e retomar integração

r

matriz de reais . Cada coluna é o vetor [t;x(t);xdot(t)] onde t é o índice do tempo para o qual a solução foi computada

Descrição

A função dassl integra a equação alegébrica-diferencial e retorna a evolução deyem pontos de tempo dados.

g(t,y,ydot)=0
y(t0)=y0  e   ydot(t0)=ydot0

Exemplos

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);
 
//utilizando um argumento extra para parâmetros
//-----------------------------------
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));

//utilizando código c
//------------
// - criando o código 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[])'
         ' {'
         '   /* primeira coluna*/'
         '   pd[0] = -0.04-*cj;'
         '   pd[1] =  0.04;'
         '   pd[2] =  1.0;'
         '    /* segunda coluna*/'
         '   pd[3] =  1.0e4*y[2];'
         '   pd[4] = -1.0e4*y[2]-2*3.0e7*y[1]-*cj;'
         '   pd[5] =  1.0;'
         '    /* terceira coluna*/'
         '   pd[6] =  1.0e4*y[1];'
         '   pd[7] = -1.0e4*y[1];'
         '   pd[8] =  1.0;'
         ' }'];
mputl([rescode;jaccode],fullfile(TMPDIR,'mycode.c')) //create the C file

// - compilando-o
ilib_for_link(['chemres','chemjac'],fullfile(TMPDIR,'mycode.c'),[],'c','',fullfile(TMPDIR,'loader.sce'));//compilando

// - ligando-o com o Scilab
exec(fullfile(TMPDIR,'loader.sce')) //linking incremental

// - chamada a dassl
y4=dassl([y0,yd0],0,t,'chemres','chemjac');
Scilab Enterprises
Copyright (c) 2011-2017 (Scilab Enterprises)
Copyright (c) 1989-2012 (INRIA)
Copyright (c) 1989-2007 (ENPC)
with contributors
Last updated:
Mon Oct 01 17:39:38 CEST 2012