Solução numérica de equações diferenciais .Uma maneira simples de transformar as derivadas parciais

  • View
    215

  • Download
    0

Embed Size (px)

Text of Solução numérica de equações diferenciais .Uma maneira simples de transformar as derivadas...

  • 13Soluo numrica de equaesdiferenciais parciais

    13.1Adveco pura: a onda cinemtica

    Considere a equao

    u

    t+ c u

    x= 0, u(x , 0) = g(x). (13.1)

    A sua soluo pode ser obtida pelo mtodo das caractersticas, e

    u(x , t) = g(x c t). (13.2)

    Seja ento o problema

    u

    t+ 2

    u

    x= 0, (13.3)

    u(x , 0) = 2x(1 x). (13.4)

    A condio inicial, juntamente com u(x , 1), u(x , 2) e u(x , 3) esto mostrados na figura13.1. Observe que a soluo da equao uma simples onda cinemtica.

    0.00

    0.50

    0 2 4 6 8 100

    1

    2

    3

    4

    u(x,

    t)

    tem

    po

    x

    Figura 13.1: Condio inicial da equao 13.3.

    233

  • Matemtica Aplicada 234

    Vamos adotar a notao

    uni u(x i, tn), (13.5)x i = ix , (13.6)tn = nt, (13.7)

    com

    x = L/Nx , (13.8)t = T/Nt (13.9)

    onde L, T so os tamanhos de grade no espao e no tempo, respectivamente, e Nx , Ntso os nmeros de divises no espao e no tempo.

    Uma maneira simples de transformar as derivadas parciais em diferenas finitasna equao (13.3) fazer

    u

    t

    i,n

    =un+1i unit

    +O(t), (13.10)

    u

    x

    i,n

    =uni+1 uni1x

    +O(x2). (13.11)

    Substituindo na equao (13.3), obtemos o esquema de diferenas finitas explcito:

    un+1i unit

    =c

    uni+1 uni12x

    ,

    un+1i = uni

    ct2x

    (uni+1 uni1), (13.12)

    (com c = 2 no nosso caso). Esse um esquema incondicionalmente instvel, e vaifracassar. Vamos fazer uma primeira tentativa, j conformados com o fracasso an-tecipado. Ela vai servir para desenferrujar nossas habilidades de programao demtodos de diferenas finitas.

    O programa que implementa o esquema instvel o onda1d-ins.py, mostradona listagem 13.1. Por motivos que ficaro mais claros na sequncia, ns escolhemosx = 0,01, e t = 0,0005.

    O programa gera um arquivo de sada binrio, que por sua vez lindo pelo pr-ximo programa na sequncia, surf1d-ins.py, mostrado na listagem 13.2. O nicotrabalho deste programa selecionar algumas linhas da sada de onda1d-ins.py;no caso, ns o rodamos com o comando

    surf1d-ins.py 3 250

    ,

    o que significa selecionar 3 sadas (alm da condio inicial), de 250 em 250 interva-los de tempo t. Observe que para isto ns utilizamos uma lista (v), cujos elementosso arrays.

    O resultado dos primeiros 750 intervalos de tempo de simulao mostrado nafigura 13.2. Repare como a soluo se torna rapidamente instvel. Repare tambmcomo a soluo numrica, em t = 750t = 0,375, ainda est bastante distante dostempos mostrados na soluo analtica da figura 13.1 (que vo at t = 4). Clara-mente, o esquema explcito que ns programamos jamais nos levar a uma soluonumrica satisfatria para tempos da ordem de t = 1!

  • 235 13.1 Adveco pura: a onda cinemtica

    Listagem 13.1: onda1d-ins.py Soluo de uma onda 1D com um mtodo explcitoinstvel

    1 #!/usr/bin/python2 # -*- coding: iso-8859-1 -*-3 # ----------------------------------------------------------4 # onda1d-ins resolve uma equao de onda com um mtodo5 # explcito6 #7 # uso: ./onda1d-ins.py8 # ----------------------------------------------------------9 from __future__ import print_function

    10 from __future__ import division11 fou = open('onda1d-ins.dat','wb')12 dx = 0.0113 dt = 0.000514 print('# dx = %9.4f' % dx)15 print('# dy = %9.4f' % dt)16 from numpy import zeros17 nx = int(10.0/dx) # nmero de pontos em x18 nt = int(1.0/dt) # nmero de pontos em t19 print('# nx = %9d' % nx)20 print('# nt = %9d' % nt)21 u = zeros((2,nx+1),float) # apenas 2 posies no tempo22 # so necessrias!23 def CI(x): # define a condio inicial24 if 0

  • Matemtica Aplicada 236

    Listagem 13.2: surf1d-ins.py Seleciona alguns intervalos de tempo da soluonumrica para plotagem

    1 #!/usr/bin/python2 # -*- coding: iso-8859-1 -*-3 # ----------------------------------------------------------4 # surf1d-ins.py: imprime em +1 sadas de5 # onda1d-ins a cada intervalos de tempo6 #7 # uso: ./surf1d-ins.py 8 # ----------------------------------------------------------9 from __future__ import print_function

    10 from __future__ import division11 from sys import argv12 dx = 0.0113 dt = 0.000514 print('# dx = %9.4f' % dx)15 print('# dy = %9.4f' % dt)16 nx = int(10.0/dx) # nmero de pontos em x17 print('# nx = %9d' % nx)18 m = int(argv[1]) # m sadas19 n = int(argv[2]) # a cada n intervalos de tempo20 print('# m = %9d' % m)21 print('# n = %9d' % n)22 fin = open('onda1d-ins.dat',23 'rb') # abre o arquivo com os dados24 from numpy import fromfile25 u = fromfile(fin,float,nx+1) # l a condio inicial26 v = [u] # inicializa a lista da "transposta"27 for it in range(m): # para instantes:28 for ir in range(n): # l vezes, s guarda a ltima29 u = fromfile(fin,float,nx+1)30 v.append(u) # guarda a ltima31 founam = 'surf1d-ins.dat'32 print(founam)33 fou = open(founam,'wt') # abre o arquivo de sada34 for i in range(nx+1):35 fou.write('%10.6f' % (i*dx)) # escreve o "x"36 fou.write('%10.6f' % v[0][i]) # escreve a cond inicial37 for k in range(1,m+1):38 fou.write('%10.6f' % v[k][i])# escreve o k-simo39 fou.write('\n')40 fou.close()

  • 237 13.1 Adveco pura: a onda cinemtica

    0.00

    0.50

    0 2 4 6 8 100.000

    0.125

    0.250

    0.375

    0.500

    u(x,

    t)

    tem

    pode

    sim

    ula

    o

    x

    Figura 13.2: Soluo numrica produzida por onda1d-ins.py, para t = 250t,500t e 750t.

    Por que o esquema utilizado em (13.12) fracassa? Uma forma de obter a resposta fazer uma anlise de estabilidade de von Neumann.

    A anlise de estabilidade de von Neumann consiste primeiramente em observarque, em um computador real, (13.12) jamais ser calculada com preciso infinita. Oque o computador realmente calcula um valor truncado uni . Por enquanto, ns svamos fazer essa distino de notao, entre u e u, aqui, onde ela importa. O erro detruncamento

    ni uni uni . (13.13)Note que (13.12) se aplica tanto para u quanto para u; subtraindo as equaes resul-tantes para un+1i e u

    n+1i , obtm-se a mesma equao para a evoluo de

    ni :

    n+1i = ni

    Co

    2(ni+1 ni1), (13.14)

    onde

    Co ctx

    (13.15)

    o nmero de Courant. Isso s foi possvel porque (13.12) uma equao linear emu. Mesmo para equaes no-lineares, entretanto, sempre ser possvel fazer pelomenos uma anlise local de estabilidade.

    O prximo passo da anlise de estabilidade de von Neumman escrever uma sriede Fourier para ni , na forma

    tn = nt,x i = ix ,

    ni =N/2

    l=1

    leatneikl x i , (13.16)

  • Matemtica Aplicada 238

    onde e a base dos logaritmos naturais, i =p1, N = L/x o nmero de pontos

    da discretizao em x , e L o tamanho do domnio em x .Argumentando novamente com a linearidade, desta vez de (13.14), ela vale para

    cada modo l de (13.16), donde

    lea(tn+t)eikl ix = le

    atneikl ix Co2

    leatneikl (i+1)x leatneikl (i1)x

    ; (13.17)

    eliminando o fator comum leatn+ikl ix ,

    eat = 1 Co2

    e+iklx eiklx

    = 1 iCo sen klx . (13.18)O lado direito um nmero complexo, de maneira que o lado esquerdo tambm temque ser! Como concili-los? Fazendo a = + i , e substituindo:

    e(i)t = 1 iCo sen klx;et

    cos(t) i sen(t)= 1 iCo sen klx; et cos(t) = 1, (13.19)et sen(t) = Co sen(klx). (13.20)

    As duas ltimas equaes formam um sistema no-linear nas incgnitas . O sis-tema pode ser resolvido:

    tg(t) = Cosen(klx) t = arctg

    Cosen(klx)

    .

    Note que 6= 0, donde et > 1 via (13.19), e o esquema de diferenas finitas incondicionalmente instvel.

    O mtodo de Lax Uma alternativa que produz um esquema estvel o mtodo deLax:

    un+1i =1

    2

    (uni+1+ uni1)Co(uni+1 uni1)

    . (13.21)

    Agora que ns j sabemos que esquemas numricos podem ser instveis, deve-mos fazer uma anlise de estabilidade antes de tentar implementar (13.21) nume-ricamente. Vamos a isto: utilizando novamente (13.16) e substituindo em (13.21),temos

    lea(tn+t)eikl ix =

    1

    2

    leatneikl (i+1)x + le

    atneikl (i1)x

    Co

    leatneikl (i+1)x leatneikl (i1)x

    ;

    eat =1

    2

    e+iklx + eiklx

    Co

    e+iklx eiklx

    ;

    eat = cos(klx) iCo sen(klx) (13.22)Ns podemos, claro, fazer a = i , mas h um caminho mais rpido: o truque perceber que se o fator de amplificao eat for um nmero complexo com mdulomaior que 1, o esquema ser instvel. Desejamos, portanto, que |eat 1|, o que s possvel se

    Co 1, (13.23)

  • 239 13.1 Adveco pura: a onda cinemtica

    que o critrio de estabilidade de Courant-Friedrichs-Lewy.A mgica de (13.21) que ela introduz um pouco de difuso numrica; de fato,

    podemos reescrev-la na forma

    un+1i unit

    =c uni+1 uni1

    2x+

    uni+1 2uni + uni12t

    =c uni+1 uni1

    2x+

    x2

    2t

    uni+1 2uni + uni1x2

    . (13.24)

    No custa repetir: (13.24) idntica a (13.21). Porm, comparando-a com (13.12)(nosso esquema instvel inicialmente empregado), ns vemos que ela tambm equi-valente a esta ltima, com o termo adicional

    x2/2t

    (uni+1 2uni + uni1)/x2. Oque este termo adicional significa? A resposta uma derivada numrica de ordem 2.De fato, considere as expanses em srie de Taylor

    ui+1 = ui +du

    d x

    i

    x +1

    2

    d2u

    d x2

    i

    +O(x2),

    ui1 = ui du

    d x

    i

    x +1

    2

    d2u

    d x2

    i

    +O(x2),

    e some:

    ui+1+ ui1 = 2ui +d2u

    d x2

    i

    x2+O(x2),

    d2u

    d x2

    i

    =ui+1 2ui + ui1

    x2+O(x2). (13.25)

    Portanto, a equao (13.24) ou seja: o esquema de Lax (13.21) pode ser inter-pretada tambm como uma soluo aproximada da equao de adveco-difuso

    u

    t+ c u

    c= D

    2u

    x2,

    com

    D =