Notebook

Equação de adveçcão¶

Como primeiro exemplo de simulação numérica, vamos considerar a equação de advecção

$$ u_t + au_x = 0. $$

In [1]:

using DifferentialEquations
using Plots

Solução analítica da equação de advecção¶

Vamos começar estudando a equação de advecção em uma dimensão:

$$ u_t + a u_x = 0 $$

Vamos considerar o domínio $x\in \mathbb{R}$ e uma condição inicial $u_0(x)$.
A solução dessa equação é simples:

$$ u(t,x) = u_0(x-at). $$

De fato,

$$ u_t(t,x) = -au_0'(x-at), \quad u_x = u_0'(x-at), \quad \Longrightarrow u_t + au_x = -au_0'(x-at) + a u_0'(x-at) = 0. $$

Portanto, a solução é o deslocamento da condição inicial $u_0$ com velocidade constante $a$.
Em aplicacões, $u_0(x)$ representa alguma quantidade no ponto $x$, digamos (densidade de) massa, calor, etc., que é então transportada com velocidade $a$.
Caso $a>0$, a solução se desloca para a direita.
Caso $a<0$, a solucão se desloca para a esquerda.

Condições de contorno em domínios limitados¶

Em domínios limitados, é preciso incluir condições de contorno.
Uma das condições mais simples é assumir condições periódicas. Nesse caso, é como se o domínio fosse a reta toda e a condição inicial e a solução fossem copiadas infinitas vezes para os dois lados, de maneira periódica.
Outras condições clássicas são a de Dirichlet, $u(t, \tilde x)=f(t)$, e a de Neumann, $u_x(t, \tilde x)=g(t)$, onde $\tilde x$ é um ponto no bordo.
Mas nesses casos, como a condição inicial é transportada para um dos dois lados, dependento do sinal de $a$, só podemos impor a condição de contorno "no lado de origem", para evitar conflitos.
Por exemplo, se o dominio é o intervalo $I=[0, L]$ e $a>0$, então só podemos impor $u(t,0)=f(t)$ ou $u_x(t,0)=g(t)$. Se $a<0$, então só podemos impor $u(t,L)=f(t)$ ou $u_x(t,L)=g(t)$.
No caso em que $f=0$, temos a condição de Dirichlet homogêna.
No caso em que $g=0$, temos a condição de Neumann homogênea.

Simulação de deslocamento para a direita¶

Vamos assumir, aqui, condições de contorno periódicas.
O intervalo espacial é $I=[0,2\pi]$, com condição inicial $u_0(x) = \sin(x)$.
No caso de deslocamento para a direita, temos $a>0$.
Nesse caso, vamos usar diferenças finitas "para frente".

$$ \partial_x u(x) \approx \frac{u(x+h) - u(x)}{h} $$

O uso de diferenças centradas ou para trás, nesse caso, nos leva a métodos instáveis.

In [2]:

"""
Operador de diferenças finitas para trás, com condições de contorno periódicas
"""
function δ⁻(u::Vector{Float64}, h::Float64, ::Val{:per})
    du = zero(u)
    for j = 2:length(u)
        du[j] = (u[j] - u[j-1])/h
    end
    du[1] = (u[1] - u[end])/h
    return du
end

Out[2]:

δ⁻

In [3]:

function dudt_advection⁻!(dudt::Vector{Float64}, u::Vector{Float64},
        p::Vector{Float64}, t::Float64)
    a, h = p
    du = δ⁻(u, h, Val(:per))
    dudt .= - a * du
    return nothing
end

Out[3]:

dudt_advection⁻! (generic function with 1 method)

In [4]:

a = 0.5 # # velocidade de propagação
L = 2π # comprimento do intervalo [0,L]
N = 60 # número de pontos da malha
h = L/(N-1) # comprimento de cada partição na malha
x = range(0.0, L, length=N) # discretização do intervalo [0,L] com N pontos, incluindo os extremos
u₀ = sin.(x) # condição inicial, com N pontos (incluindo os dois extremos)
p = [a, h] # parâmetros
Tf = L/a # tempo final (volta completa)
τ = 0.1 # intervalo de tempo para guardar a solução
tspan = (0.0,Tf) # intervalo de tempo
prob = ODEProblem(dudt_advection⁻!, u₀, tspan, p, saveat = τ)
nothing

In [5]:

sol = solve(prob, Tsit5())
sol.retcode

Out[5]:

:Success

In [6]:

sol

Out[6]:

retcode: Success
Interpolation: 1st order linear
t: 127-element Vector{Float64}:
  0.0
  0.1
  0.2
  0.30000000000000004
  0.4
  0.5
  0.6
  0.7000000000000001
  0.8
  0.9
  1.0
  1.1
  1.2000000000000002
  ⋮
 11.5
 11.6
 11.700000000000001
 11.8
 11.9
 12.0
 12.1
 12.200000000000001
 12.3
 12.4
 12.5
 12.566370614359172
u: 127-element Vector{Vector{Float64}}:
 [0.0, 0.1062934856473654, 0.21138262362962432, 0.31407671202194876, 0.41321218576837815, 0.5076658003388399, 0.5963673585385013, 0.6783118362696161, 0.7525707698561385, 0.818302775908169  …  -0.8183027759081689, -0.7525707698561385, -0.678311836269616, -0.5963673585385014, -0.5076658003388399, -0.4132121857683782, -0.3140767120219489, -0.21138262362962446, -0.1062934856473656, -2.4492935982947064e-16]
 [-0.010057738237594616, 0.06501077026477721, 0.163223305414004, 0.2657657931936173, 0.36630895192625346, 0.46282364675910964, 0.5541060018585185, 0.63911103850578, 0.7168747543854092, 0.7865159965125561  …  -0.8437100125058529, -0.7823996938086538, -0.7122244739495479, -0.6339794662101336, -0.5485512177932172, -0.45690766488875856, -0.3600871655700709, -0.25918673478194343, -0.15534961472297223, -0.049752321454633215]
 [-0.03491395517909535, 0.032218432474368085, 0.11969124185206444, 0.2183941765271466, 0.3190127507882637, 0.4171101673112567, 0.5106947341982397, 0.5985272095014393, 0.6795828289036199, 0.7529390619287905  …  -0.8668811992358952, -0.8101218009089822, -0.7441833988695238, -0.669813101560275, -0.5878535541666131, -0.49923339110988135, -0.4049567142483335, -0.3060917160018589, -0.20375857630503935, -0.09911677052088785]
 [-0.06885415877009128, 0.0004261340836273676, 0.0805296841073726, 0.17339504970534578, 0.27213834813983, 0.3708600471578879, 0.4662983918605026, 0.5566784832932635, 0.6407979346460463, 0.7176653235382806  …  -0.8877704920630689, -0.8356792274505559, -0.7741193832705516, -0.7037884571556812, -0.625483326601638, -0.5400912200267426, -0.4485796641194421, -0.3519855213743294, -0.2514032420257676, -0.14797246348941526]
 [-0.10825998537024041, -0.03336287717760092, 0.043774944564390905, 0.131134291134544, 0.22652497314434966, 0.3246089940667689, 0.42118870442562606, 0.5137206892891585, 0.6006352974750342, 0.6807938343162148  …  -0.906337973219629, -0.8590198078516809, -0.8019686048138744, -0.7358307770331312, -0.6613556925255426, -0.5793871837573139, -0.4908539866857425, -0.39675921780972, -0.2981690084589743, -0.19620042510009775]
 [-0.15082405540983235, -0.06995238702801329, 0.007571289970564923, 0.09109733031221802, 0.18264136226457892, 0.2789388511374285, 0.3757518397366884, 0.46987726398827556, 0.5592383361909751, 0.6424348198831817  …  -0.9225497192282772, -0.8800971744123758, -0.8276727769413948, -0.7658705161864063, -0.6953906364566954, -0.6170317029519301, -0.5316815536993611, -0.44030723999420845, -0.3439440693220756, -0.24368387491144677]
 [-0.19506280973116166, -0.1091516811682829, -0.02932025439806001, 0.05238423456058722, 0.14049687283471254, 0.2342961351225284, 0.33042462177298787, 0.4254453572320646, 0.5167955957884462, 0.6027196311450813  …  -0.9363778532477574, -0.898870841542239, -0.8511792640314815, -0.7938434853452718, -0.7275131427663369, -0.6529397855803579, -0.5709683597222529, -0.4825276342011144, -0.38861967777621725, -0.29030850511757145]
 [-0.24000087685544447, -0.15040974717974073, -0.0675107244153527, 0.014105302140080303, 0.09975092210464653, 0.19088477714543062, 0.2855996469413282, 0.38076889856495205, 0.4735474420264886, 0.5618116978027428  …  -0.9478005834138212, -0.9153062767881148, -0.8724411843085672, -0.8196909846009951, -0.7576533585059949, -0.6870312171162, -0.6086247375127336, -0.5233222964391426, -0.4320904046379159, -0.3359627558975765]
 [-0.28499548641412303, -0.19306848542937066, -0.107153421196027, -0.024439700056455183, 0.059892192392528104, 0.14867050568774398, 0.2415519209789444, 0.3361874838599748, 0.4297728796645624, 0.5199112783232672  …  -0.9568022227181977, -0.9293749524157447, -0.8914174925430326, -0.843359916287451, -0.7857467351908072, -0.7192307296159881, -0.6445655524846017, -0.5625971900882375, -0.4742543767257438, -0.380538071772513]
 [-0.32959396973260013, -0.23653941069701864, -0.14810989353777096, -0.06366402363205848, 0.02040072219643715, 0.10742699511421148, 0.19839397858821978, 0.2919838099937205, 0.38576506390988546, 0.4772530326013572  …  -0.9633731993902743, -0.9410543893969674, -0.9080730571455374, -0.8648028943090658, -0.8117341691223683, -0.7494681714469985, -0.6787103999168476, -0.600262568356723, -0.5150135220443646, -0.4239291667382822]
 [-0.37348672904475805, -0.28029410571460744, -0.19007591015495548, -0.10375609665537283, -0.019169082435306434, 0.066835968594659, 0.15609871625047522, 0.24835092253193028, 0.34179483731360727, 0.4340883168993156  …  -0.9675100425850026, -0.9503281739077059, -0.9223787072881647, -0.8839783213512176, -0.8355621078135452, -0.7776786417206679, -0.7109837658724775, -0.636233159862253, -0.5542737779249187, -0.46603425260712966]
 [-0.4164302549546262, -0.3239177620076702, -0.2327025810905798, -0.1446944931730752, -0.05910332355945962, 0.026563442849627056, 0.11451518054357172, 0.2053666551444189, 0.2980799876998681, 0.390667574765698  …  -0.9692153659541224, -0.9571859762039943, -0.9343112868713654, -0.9008504773266359, -0.8571826719122091, -0.8038026443121993, -0.7413152115688935, -0.6704283812639819, -0.5919453295097618, -0.5067553006428421]
 [-0.4582467054271522, -0.36706402467655974, -0.275613097741568, -0.18633219429309494, -0.09956287083578431, -0.013688421992499738, 0.07343420686201867, 0.16301645286034597, 0.25476976514400407, 0.3472154778137774  …  -0.9684978232661664, -0.9616235334388481, -0.9438536647324002, -0.9153895568597458, -0.8765537195367621, -0.827786178318332, -0.7696394889329812, -0.7027724766056839, -0.6279427713048108, -0.5459982234918539]
 ⋮
 [0.4626925757171845, 0.5197084598679994, 0.5712692134162024, 0.6168181822798823, 0.6557570122794477, 0.6876770836710681, 0.7120936112055404, 0.7287874503415496, 0.7374084585409517, 0.7379581253728655  …  -0.2516496284340157, -0.17950349669530338, -0.10553336202045746, -0.03047221894533244, 0.04493290450785265, 0.11993694943535411, 0.19378101310349005, 0.26572207380770535, 0.33500253217754505, 0.4009042983578619]
 [0.43323611378949567, 0.49243203232941896, 0.5464619847195149, 0.5947919900783437, 0.6367414299742509, 0.6719429633371161, 0.6997882078263767, 0.7201301041886025, 0.7324443863971354, 0.736848945840973  …  -0.2839546653703741, -0.21308309009825352, -0.1400497650095259, -0.06558047822270283, 0.00957956952938366, 0.0846918955320644, 0.1589942087366732, 0.23175070105783024, 0.30219145855632507, 0.3696073260204317]
 [0.402960729393416, 0.46415762277161926, 0.5205432836689168, 0.5714745996292214, 0.6163637082377562, 0.6546664394538947, 0.6859294664462038, 0.7097364668886798, 0.7258100954795961, 0.7338795142215306  …  -0.3154698404686832, -0.24602840758922911, -0.17408853504098395, -0.100373876188072, -0.02561922892084522, 0.049433011607817996, 0.12403341291979374, 0.19743216443569345, 0.2688786631323459, 0.3376317112423159]
 [0.3719246485200948, 0.4349805667093558, 0.49355332341510333, 0.5469693474200821, 0.5946356886155895, 0.6359621875061007, 0.6704889221502508, 0.6977374622172758, 0.7174246552074339, 0.7292030044952984  …  -0.3461264803603616, -0.27826644005799156, -0.20757715885140077, -0.13477546799717816, -0.06058962982123367, 0.014240837043848058, 0.0889714103498765, 0.16285099880519022, 0.2351323912394057, 0.30506678813024307]
 [0.34020451894955983, 0.4049734941515185, 0.4655645532962203, 0.521341532561757, 0.5716221641016415, 0.6158839844287055, 0.6535201142966574, 0.6841712846956429, 0.7073272034102265, 0.722838742542415  …  -0.3758564770541817, -0.30972648144288145, -0.24044272048507254, -0.1687108836168709, -0.0952557762125724, -0.02080816115698858, 0.053885782233547154, 0.12808446754071243, 0.2010307038617931, 0.27198885556511576]
 [0.3078869076807189, 0.37419710670496265, 0.43666770671090566, 0.49463565268738896, 0.5474193324310155, 0.5944506449000322, 0.6351272962126203, 0.6690200352924757, 0.6956342992474922, 0.7147206569015206  …  -0.4045922698795907, -0.34034002772151495, -0.27261243033777866, -0.20210754144178367, -0.1295411151726065, -0.05563977275079725, 0.018856285917847132, 0.09320609043457158, 0.16665660355894282, 0.2384675302455796]
 [0.27504474605097606, 0.3427318176558183, 0.40692930723637954, 0.466931011872749, 0.522082076962302, 0.5717388309299939, 0.6153474089255594, 0.652356659710928, 0.6823582656934998, 0.7049177763688533  …  -0.4322695778456389, -0.3700388897444762, -0.304015899880477, -0.2348923357530011, -0.16337171932436925, -0.09017765665539308, -0.016041460144607996, 0.05829446748902336, 0.13208561944175704, 0.20458286482771043]
 [0.24175073794727003, 0.3106593808518278, 0.3764156519279145, 0.438309892875599, 0.4956645881943721, 0.5478308070040124, 0.5942152193825742, 0.6342637256626067, 0.6675059975895505, 0.6935147251685762  …  -0.4588280264340868, -0.398755353115982, -0.3345849250952787, -0.26699269264869585, -0.19667486500835263, -0.12434586151215882, -0.050732521496541845, 0.023428078279629515, 0.09739288134558519, 0.17041530069832192]
 [0.20808360692835845, 0.278053748948399, 0.34520476697177355, 0.40884072564798424, 0.4682422180216393, 0.5227848528199837, 0.5718012765816797, 0.6147838354544468, 0.6511415246987899, 0.6805325648306076  …  -0.4842064423143692, -0.4264258514993202, -0.3642512273807662, -0.2983385787292899, -0.22937757629759, -0.15807052301827698, -0.08514088417491864, -0.011317253646062056, 0.06265630959920761, 0.13604088786385748]
 [0.17412401071548386, 0.24498660642566886, 0.3133788433537509, 0.37858809730723114, 0.4398981396037934, 0.49665250992753623, 0.5481895479190471, 0.5939478826432547, 0.6333503789953584, 0.6659727774989552  …  -0.5083446103193974, -0.4529898433902406, -0.3929474251210088, -0.3288617991969665, -0.26140742672410716, -0.19127903168654795, -0.11919060189604765, -0.04586685781737683, 0.027954484464797202, 0.10153420832610925]
 [0.13994636098692478, 0.21153849366891359, 0.2810087219858571, 0.34763343592108725, 0.41069520410733057, 0.4695152832451113, 0.5234288656700916, 0.5718370343823277, 0.6141599600471367, 0.6499167464715976  …  -0.5311902709721679, -0.47838543219437507, -0.42061064771959183, -0.3584937531864719, -0.2926948696185687, -0.2238981335153392, -0.15280832585832235, -0.08014423672120385, -0.006637618782509085, 0.06697412885458202]
 [0.11717981009678732, 0.18916610244313495, 0.2592614358228818, 0.32673777467567405, 0.39087475176395786, 0.45098372436041867, 0.5063935384134988, 0.5564883474460175, 0.600678647774868, 0.6384585180661744  …  -0.5456136470047357, -0.4945673816268969, -0.43836994710398103, -0.37763626498392044, -0.3130167304324964, -0.2451871887752768, -0.1748461025661458, -0.10270736422191966, -0.029498793410425217, 0.04404446832969359]

In [7]:

plot(x, sol.u[end], ylims=(-1.1, 1.1), label="solução no tempo final $(sol.t[end])")

Out[7]:

In [8]:

anim = @animate for (t,u) in zip(sol.t, sol.u)
    plot(x, u, ylims=(-1.1, 1.1), label="t=$(round(t,digits=2))",
        title="Uma solução da equação de adveção com a=$a", titlefont=10)
end
gif(anim, "../../../assets/attachments/img/anim_adveccao_dir.gif", fps = 20)
nothing

┌ Info: Saved animation to 
│   fn = /Users/rmsrosa/Documents/git-repositories/modelagem_matematica/_assets/attachments/img/anim_adveccao_dir.gif
└ @ Plots /Users/rmsrosa/.julia/packages/Plots/MzlNY/src/animation.jl:130

advection

Deslocamento para a esquerda¶

In [12]:

"""
Operador diferenças finitas para frente, com condições de contorno periódicas.
"""
function δ⁺(u, h, ::Val{:per})
    du = zero(u)
    for j = 1:length(u)-1
        du[j] = (u[j+1] - u[j])/h
    end
    du[end] = (u[1] - u[end])/h
    return du
end

Out[12]:

δ⁺

In [13]:

function dudt_advection⁺!(dudt, u, p, t)
    a, h = p
    du = δ⁺(u, h, Val(:per))
    dudt .= - a * du
    return nothing
end

Out[13]:

dudt_advection⁺! (generic function with 1 method)

In [14]:

a = -0.5 # velocidade de propagação
p = [a, h] # parâmetros
prob = ODEProblem(dudt_advection⁺!, u₀, tspan, p, saveat = τ)
sol = solve(prob, Tsit5())
sol.retcode

Out[14]:

:Success

In [15]:

anim = @animate for (t,u) in zip(sol.t, sol.u)
    plot(x, u, ylims=(-1.1, 1.1), label="t=$(round(t,digits=2))",
        title="Uma solução da equação de adveção com a=$a", titlefont=10)
end
gif(anim, "../../../assets/attachments/img/anim_adveccao_esq.gif", fps = 20)
nothing

┌ Info: Saved animation to 
│   fn = /Users/rrosa/Documents/git_repositories/modelagem_matematica/_assets/attachments/img/anim_adveccao_esq.gif
└ @ Plots /Users/rrosa/.julia/packages/Plots/MzlNY/src/animation.jl:130

Out[15]:

advection left

Exercícios¶

Implemente o operador diferenças finitas para trás com condições de contorno de Dirichlet homogênea $u(0)=0$ e faça a simulação para algum parâmetro $a>0$.
Implemente o operador diferenças finitas para frente com condições de contorno de Dirichlet homogênea $u(0)=0$, faça a simulação no caso $a>0$ e observe a instabilidade do método nesse caso.
Implemente o operador diferenças finitas para frente com condições de contorno de Neumann homogênea $u'(L)=0$ e faça a simulação para algum parâmetro $a<0$.
Implemente o operador diferenças finitas para trás com condições de contorno de Neumann não-homogênea $u'(0)=b$, com um novo parâmetro $b>0$, e faça a simulação para algum parâmetro $a>0$.
Implemente o operador diferenças centradas, com condições de contorno periódicas, faça a simulação para algum parâmetro $a>0$ e observe a instabilidade do método.
Com condições de contorno periódicas, em um intervalo $I$, mostre que a "massa" $\int_I u(t,x)\;\mathrm{d}x$ é constante em $t$.
Da mesma forma, mostre que $\int_I u^n(t,x)\;\mathrm{d}x$ também é constante em $t$, para todo $n\in \mathbb{N}$.
Calcule $m[k] = \sum_i u(k\tau,x_i)h$, $k=1, \ldots$, no exemplo de condições periódicas de contorno e observe a evolução dessa quantidade através de um gráfico.