Métodos computacionais / 2019

Davi C. Rodrigues
Aula 11

Item 1 - revisão dos exercícios para hoje.

As notas desta aula 11 não seguem nenhum livro em especial, foram inspiradas no help do mathematica, comentários no stackexchange e outros.

Detalhes adicionais sobre o NDSolve

Dando continuidade ao exercício enviado...

Sobre o help do Mathematica

O help do Mathematica é especialmente desenvolvido para resolução numérica de equações diferencais. Além do help do próprio NDSolve, ver esses dois tutoriais

Curiosidade: os links acima foram gerados a partir de...

Hyperlink["Tutorial introdutório","paclet:tutorial/NumericalSolutionOfDifferentialEquations", BaseStyle->{"Text"}]
Hyperlink["Tutorial avançado","paclet:tutorial/NDSolveOverview",BaseStyle->{"Text"}]

Um exemplo 1D de problema de condições de contorno (inspirado no help do Mathematica).

Primeiramente, tentemos resolver a eq. abaixo analiticamente.

xsol[t_]= x[t]/.DSolve[{x''[t] + Sin[x[t]] == 0}, x[t],t]

Duas soluções foram encontradas, mas sem garantia de que não há outras.

O Solve é em geral um método menos robusto que o Reduce de resolução de sistema de equações algébricas. O DSolve pode chamar o Solve para resolver algum problema algébrico, mas neste caso aparece uma emnsagem indicando que seria indicado tentar o Reduce. Segue um breve exemplo do uso de Reduce vs Solve.

Exemplo Solve vs Reduce

Solve[{(x^2) + b x + π == y, Sin[y] ==  x},x]

O Solve tenta resolver a equação para x para quaisquer valores das outras variáveis e sem usar condições. O último comportamente pode ser alterado usando MaxExtraConditions -> 1 (o padrão é zero).

Solve[{(x^2) + b x + π == y, Sin[y] ==  x},x, MaxExtraConditions->1]

Reduce usa condições extras de como padrão, e em geral é mais pesado para execução que o Solve,

Reduce[{(x^2) + b x + π == y, Sin[y] ==  x}, x]

A fim de evitar que o DSolve tente usar funções inversas, e nos dê a resposta que foi possível encontrar com a integração, pode-se usar

SetOptions[Solve,InverseFunctions-> False];
DSolve[{x''[t] + Sin[x[t]] == 0}, x[t],t]
SetOptions[Solve,InverseFunctions-> Automatic];

Com o resultado acima, pode-se finalmente aplicar Reduce, ao invés de Solve e....

Reduce[((4 ((EllipticF[((x[t])/2),(4/(2+C[1]))])^2))/(2+C[1]))==(((t+C[2]))^2), x[t]]

Como pode ver, neste caso mesmo o Reduce não encontra solução. Ressalta-se que em alguns casos um NDSolve leva a um Solve sem resposta, mas um Reduce pode chegar até a resposta desejada, logo o método acima é útil em várias situações.

Uma opção é usar métodos numéricos como FindRoot para conseguir inverter essa função, isolando x[t]. Mas não vamos seguir nesse ramo. Busquemos solução numérica diretamente da solução da eq. diferencial,

Usando o NDSolve diretamente

NDSolve vs NDSolveValue: só diferem na forma de expressar a solução final.

NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0}, x[t],t]
NDSolveValue[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0}, x[t],t]

sol1 = x[t] /. First@NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0}, x[t],t];
Plot[sol1, {t,0,10}]

A solução encontrada é só uma. Não poderia haver mais soluções? -- Lembrar que trata-se de problema de condição de contorno, nõa de condições iniciais.

Certamente há pelo menos mais uma solução, a trivial, em que x[t] = 0. Como encontrar outras soluções? Elas diferem no valor da derivada num ponto. Simplesmente aidicionar uma condição extra na derivada leva a problemas...

NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0, x'[0] == 0},x,{t,0,10}]

Outra opção é forçar o problema proposto transformando-o num problema de condições iniciais

Clear[sol]
sol = NDSolveValue[{x''[t] + Sin[x[t]] == 0, x[0]== 0, x'[0] == 0},x,{t,0,10}]
Plot[sol[t],{t,0,10}, Axes->False, Frame-> True]

Ou seja, a solução trivial foi encontrada. Encontrar outras soluções é fácil, é só trocar a deriva de x em t=0, contudo essas soluções em geral não satisfazem o problema da condição de contorno dada originalmente,

Clear[sol]
sol = NDSolveValue[{x''[t] + Sin[x[t]] == 0, x[0]== 0, x'[0] == 1},x,{t,0,10}]
Plot[sol[t],{t,0,10}, Axes->False, Frame-> True]
sol[10]

A solução acima não satisfaz a condição de contorno x[0]== x[10]==0.

Sobre o método de tiro/chute: “Shooting Method”

A descrição geral desse método no Mathematica pode ser vista (HyperlinkTemplate["aqui", paclet:tutorial/NDSolveBVP, paclet:tutorial/NDSolveBVP, Link, HyperlinkActive, ]) . O ponto chave é usar as condições de contorno desejadas e deixar o Mathematica buscar bela derivada em x[0] tal que x[10]==0. Para isso, tal como num FindRoot, indica-se um ponto de partida na busca pela solução.

sol=NDSolveValue[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0},x,t, Method->{"Shooting", "StartingInitialConditions"->{ x'[0] == 2}}]
Plot[sol[t],{t,0,10}, Axes->False, Frame-> True]

sols = First[NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0},x,t, Method->{"Shooting", "StartingInitialConditions"->{ x'[0] == #}}]]& /@Range[-5,5,0.1];
Plot[Evaluate[x[t] /. sols], {t,0,10}, PlotRange->All]

Questão: encontramos 7 soluções diferentes. O método não garante que essas sejam todas as possíveis, pode haver outras com essas condições de contorno?

Resolução numérica das equações de Friedmann

As equações de Friedmann para tempo conforme, com constante cosmológica e poeira podem ser escritas da seguinte forma

Clear[H, Λ, a, ϵ, ϵ0]
Feq1[t_]  = ((H[t])^2) - (1/3)Λ ((a[t])^2)== (8/3) π G ((a[t])^2) ϵ[t];
Feq2[t_] =  - 2 H'[t]-((H[t])^2) + Λ ((a[t])^2)== 8 π G ((a[t])^2) p[t];
H[t_] = a'[t]/a[t];
ϵ[t_] = ϵ0/((a[t])^3);
p[t_] = 0;
G=1; (*Escolha de sistema de unidades*)
(*Acima está sendo usado G =1*)

DSolve[Feq1[t], a[t],t]

Se tentamos impor a condição de contorno a[0] == 0, com Λ=0,

Block[{Λ=0},DSolve[{Feq1[t], a[0]==0}, a[t],t]]

Note que em a[t]=0 há uma aparente singularidade da EDO que pode ser removida multiplicando tudo por ((a[t])^2) (obs: no tempo físico a’[0] é infinito, mas no tempo conforme não; em ambos os casos H[0] é infinito). Λ≠0 leva a uma expressão que o mathematica não consegue resolver nessa forma, assim por enquanto não vamos considerar esse caso.

Pode-se notar que as equações são invariantes por um fator multiplicativo constante de a[t], desde que o tempo t seja reescalado apropriadamente (a saber t -> C t => a -> a/C). E como t só aparece nas derivadas, as equações de Friedmann são invariantes por translação temporal, logo t=0 ou t=1 não têm nehum valor intrínseco. Logo, sem perda de generalidade, podemos usar a translação temporal para fazer com que t=0 seja o instante em que a=0 (supondo que haja um instante que satisfaça a=0). Podemos ainda usar a invariância por escalonamento de a[t] por constante para fazermos a[t_hoje] = 1. Note que não temos liberdade de escolher quanto vale t_hoje, dado que já usamos a liberdade de translação temporal para fixar o significado de t=0.

Com as convenções acima, parece natural usar a[0]=0 como condição inicial/de contorno. E dado que a primeira eq. de Friedmann é de primeira ordem, nem é preciso especificar a’[t]. Entretanto, a própria equação indica que a’[0] é infinito, logo não dá para começar uma integração numérica a partir desse ponto. Usando a segunda eq. de Friedmann também não é de ajuda para esse problema, pois o valor de a’[0] teria de entrar como condição inicial.

Para resolver numericamente o caso com Λ=0, temos de especificar ϵ0. Podemos escolher ϵ0=1 (note que ele só contribui com um fator multiplicativo para a[t], e nao fixamos ate o momento a constante multiplicativa de a[t], logo essa escolha pode ser feita sem perda de generalidade).

Clear[sol]
sol[t_]= a[t]/.Block[{ϵ0=1, Λ=0},NDSolve[{Feq1[t], a[0]==0}, a[t],{t,0,5}]]
Plot[sol[t],{t,0,5}]

Ambas as soluções correspondem à solução trivial, a[t] = 0! Como encontrar as não triviais?

Duas possibilidades: 1. abandona-se a idéia de usar a[0]=0 como uma condições de contorno (às vezes é útil, mas não vamos seguir essa direção); 2. Resolver a primeira eq. de Friedmann usando uma pequena e atingível variação na condição inicial, por exemplo usa-se a[0] = (10^(-10)). A opção 2 possibilita também usar a segunda eq. de Friedmann, mas para isso, a primeira tem de ser usada para encontrar a condição inicial para a’[0].

Clear[sol]
sol[t_]= a[t]/.Block[{ϵ0=1, Λ=0},NDSolve[{Feq1[t], a[0]==(10^(-10))}, a[t],{t,0,5}]]
Plot[sol[t],{t,0,5}]

Foi encontrada duas “soluções”, a primeira na verdade não tem sentido, é complexa.

sol[0.05][[1]]

Options@NDSolve

Clear[sol, sol1,sol2]
{sol1[t_], sol2[t_]}= a[t]/.Block[{ϵ0=1, Λ=0},NDSolve[{Feq1[t], a[0]==(10^(-10))}, a[t],{t,-5,5},  WorkingPrecision->20, AccuracyGoal->∞, PrecisionGoal->5]]
Show[{Plot[sol1[t],{t,-5,-0.001}], Plot[sol2[t],{t,0,5}, PlotStyle->Orange]}, PlotRange-> All]

Uma é a solução decrescente, a outra a crescente; ambas contidas na solução anlítica.

Às vezes é relevante verificar quais os pontos da EDO foram efetivamente calculados. O uso de EavaluaitonMonitor isoladamente é útil para contar o total de pontos usados. A combinação de EvaluationMonitor com Reap & Sow mostra todos os pontos.

out  = Block[{ϵ0=1, Λ=0, c=0},NDSolve[{Feq1[t], a[0]==(10^(-10))}, a[t],{t,-5,5},  WorkingPrecision->20, AccuracyGoal->∞, PrecisionGoal->5, EvaluationMonitor:> c++]; c]

Ou seja, 318 pontos foram calculados.  Abaixo, usa-se Reap & Sow,

out  = Reap@Block[{ϵ0=1, Λ=0},NDSolve[{Feq1[t], a[0]==(10^(-10))}, a[t],{t,-5,5},  WorkingPrecision->20, AccuracyGoal->∞, PrecisionGoal->5, EvaluationMonitor:> Sow[t]]]

Pode-se fazer um ListPlot pra saber como foi a distribuição dos pontos usados.

ListPlot[Sort[out[[2,1]]], PlotRange->All]

Solução com Λ

Clear[sol, sol1,sol2]
{sol1[t_], sol2[t_]}= a[t]/.Block[{ϵ0=1, Λ=(10^(-3))},NDSolve[{Feq1[t], a[0]==(10^(-10))}, a[t],{t,-5,5},  WorkingPrecision->20, AccuracyGoal->∞, PrecisionGoal->5]]
Show[{Plot[sol1[t],{t,-5,0}], Plot[sol2[t],{t,0,5}, PlotStyle->Orange]}, PlotRange-> All]

Clear[sol, sol1,sol2]
{sol1[t_], sol2[t_]}= a[t]/.Block[{ϵ0=1, Λ=(10^(-2))},NDSolve[{Feq1[t], a[0]==(10^(-10))}, a[t],{t,-5,5},  WorkingPrecision->20, AccuracyGoal->∞, PrecisionGoal->5]]
Show[{Plot[sol1[t],{t,-5,0}], Plot[sol2[t],{t,0,5}, PlotStyle->Orange]}, PlotRange-> All]

Questão: Considere Λ=(10^(-2)). Vc é capaz de integrar este caso no mesmo intervalo? Dica: entenda primeiro o caso ϵ0=0 e Λ≠0. Lembrar que estas contas estão sendo feitas no tempo conforme.

Fisicamente, o ramo da solução de interesse é a crescente, a outra é a inversão temporal dessa (que sempre será encontrada, dado que as eqs. são invariantes por t -> -t).

Numericalmente, é agora fácil de encontrar soluções para universos com outros fluidos adicionados, em particular, de interesse físico, um de radiação (gás de fótons, dado por (Subscript[ϵ, rad]) = (Subscript[ϵ, rad])/(a^4) e (Subscript[p, rad]) = (Subscript[ϵ, rad])/3).

Exemplo  de problema de contorno em EDP

Laplacian[u[x,y,z],{x,y,z}]

Primeiramente, um breve comentário sobre solução analítica da eq. de Laplace em duas dimensões

DSolve[Laplacian[u[x,y],{x,y}]==0,u[x,y],{x,y}]

Note que acima C1 e C2 não são constantes, mas sim funções arbitrárias dos argumentos ix+y e ix-y. Usando ix como tempo, o que temos são duas ondas em um espaço unidimensional, cada uma se propagando em sentidos opostos. É fácil verificar que a solução acima realmente se aplica para funções arbitrárias.

Eq. de Poisson com condições tipo Dirichlet,

ufun=NDSolveValue[{-Laplacian[u[x,y],{x,y}]==10, u[x,0]==0,u[x,1]==1},u,{x,0,1},{y,0,1}];
Plot3D[ufun[x,y],{x,0,1},{y,0,1}]

Eq. de Helmholtz também com condições de Dirichlet

ufun=NDSolveValue[{-Laplacian[u[x,y],{x,y}]==10 u[x,y], u[x,0]==0,u[x,1]==1},u,{x,0,1},{y,0,1}];
Plot3D[ufun[x,y],{x,0,1},{y,0,1}]

Às vezes é mais conveninete usar DirichletCondition. No primeiro caso abaixo, não há um mérito aparente, mas veja o segundo

ufun=NDSolveValue[{-Laplacian[u[x,y],{x,y}]==10, DirichletCondition[u[x,y]==0, y==0],DirichletCondition[u[x,y]==1 , y==1]},u,{x,0,1},{y,0,1}];
Plot3D[ufun[x,y],{x,0,1},{y,0,1}]

ufun=NDSolveValue[{-Laplacian[u[x,y],{x,y}]==10, DirichletCondition[u[x,y]==0, True]},u,{x,0,1},{y,0,1}];
Plot3D[ufun[x,y],{x,0,1},{y,0,1}]

Acima, usa-se u[x,y]==0 para todos os pontos do limite de integração.

Questão: As soluções acima encontradas são as únicas para as dadas condições de contorno? Como encontrar as demais?

Equação diferencial com condições de Dirichlet em região não trivial

Ω=ImplicitRegion[((x-5)^2+(y-5)^2>3^2),{{x,0,5},{y,0,10}}];
RegionPlot[Ω]

op=-Laplacian[u[x,y],{x,y}]-20;
(Subscript[Γ, D])={DirichletCondition[u[x,y]==0, x==0&&8≤y≤10],
DirichletCondition[u[x,y]==100,(x-5)^2+(y-5)^2==3^2]};
ufun=NDSolveValue[{op==0,(Subscript[Γ, D])},u,{x,y}∈Ω, WorkingPrecision->100, PrecisionGoal->10, AccuracyGoal->∞]

Plot3D[ufun[x,y],{x,0,5},{y,0,10}, WorkingPrecision->30, MaxRecursion->6]

As “falhas” acima claramente são efeitos numéricos do Plot e não da solução da eq. diferencial, pois trocar parâmetreos de precisão do plot alteram essas falhas. Eliminá-las requer ou aumento substancial de maxrecursion (logo mais tempo), ou indicação de região que  deve ser plotada, como abaixo.

Plot3D[ufun[x,y],{x,y}∈Ω, WorkingPrecision->30]

ContourPlot[ufun[x,y],{x,y}∈Ω]

Questão: Reconsidere um dos exemplos anteriores de EDP, mas use condições de contorno tipo Neumann.

Mais sobre localização de expressões

Block & variations

Block[{a,b},
a= 2;
b= π;
fblock[x_]:= a x + b
];

Module[{a,b},
a= 2;
b= π;
fmodule[x_]:= a x + b
];

With[{ a=2, b=π},
fwith[x_] := x a + b
]

fblock[c]
fmodule[c]
fwith[c]

Block[{a,b},
a[x_]:= 55 x;
b= π;
gblock[x_]:= a[x] x + b
];

Module[{a,b},
a[x_]:= 55 x;
b= π;
gmodule[x_]:= a[x] x + b
];

With[{a[x_]= 55 x, b=π},
gwith[x_] := x a[x] + b
]

gblock[c]
gmodule[c]
gwith[c]

f[x_]:=x
Block[{f},
Print@DownValues[f];
f[x_,y_]:=x y;
Print@DownValues[f];];
Print@DownValues[f];

Nota sobre índices mudos

Atenção. O caso abaixo, e variações (como Sum, Do, Table...), requer atenção. O bloqueio da variável em Integrate não é como o gerado por um Module, mas sim como um Block.

f[x_] := Integrate[x, {a, 0,1}]
f[g]
f[a]

Qualquer que seja o argumento de f, ele será integrado em “a”. Esse comportamento pode ou não ser desejado, assim como nem sempre Block gera a localização mais indicada e às vezes Module é o mais indicado.

Nota-se que o “a” em Integrate está localizado no sentido de um Block (não interfere com nenhum “a” definido antes e nem deixa traço da operação de Integrate na variável “a”). Contudo, qualquer função de “a”  no argumento de f terá comportamente diferente com respeito a o uso de outra variável.   A definição abaixo via Module evita esse comportamento (que nem sempre é desejado).

Clear[f];
fm[x_] := Module[{a},Integrate[x, {a, 0,1}]];
fm[g]
fm[a]

Uma outra função semelhante a Block, mas que está “escondida” no contexto Internal`. Mais sobre contextos logo abaixo.

f[x_]:=x
Internal`InheritedBlock[{f},Print@DownValues[f];
f[x_,y_]:=x y;
Print@DownValues[f];];
Print@DownValues[f];

Não é  raro encontrar exemplos no stackexchange usando esse Internal`InheritedBlock. Ele é especialmente útil quando se deseja fazer uma modificação numa função que já existe, deixanda a possibilidade de facilmente voltar para a definição original. Por exemplo

originalPath  = $Path;
Internal`InheritedBlock[{$Path},Print[$Path == originalPath];
AppendTo[$Path, "nova pasta aqui"];
Print[$Path == originalPath];];
Print[$Path == originalPath]

O False acima indica que a definição de $Path foi alterada, e o True no final indica que voltou à definição original.

O Block usual leva a uma inconsistência...

originalPath  = $Path;
Block[{$Path},Print[$Path == originalPath];
AppendTo[$Path, "nova pasta aqui"];
Print[$Path == originalPath];];
Print[$Path == originalPath]

Uma aplicação comum é com SetOptions, enquanto se bloqueia uma dada função. Se a função estivesse bloqueada no sentido do Block (ou Module), a definição original da função não seria acessível.

Context

O contexto primário do Mathematica é o Global`

Context[]

Localizações como Module e Block são independentes do contexto,

Module[{x}, Context[]]
Block[{x}, Context[]]

Símbolos definidos numa sessão que acabou de ser aberta no Mathematica serão alocados em Global`,

x = π;
Context[x]

Variáveis de um contexto são chamadas usando o nome do contexto seguido do nome da variável naquele contexto. Assim, como o atual contexto é o Global`, tanto faz usar

x
Global`x

Mas podemos definir outro x em outro contexto...

x = π
UpsideDown`x = 1/π

x
UpsideDown`x
x UpsideDown`x

Context[UpsideDown`x]
Context[x]

$Context

Note que várias funções do Mathematica não estão definidas em Global`, mas em System`.

Context[π]
Context[Solve]

O motivo de podermos acessar a essas expressões sem precisar digitar System`Pi, como abaixo,

System`Pi
asdfasdf`Pi

se deve ao contexto System encontrar-se na $ContextPath (assim como Global`),

$ContextPath

Adicionar UpsideDown` à $ContextPath vai trazer problemas, pois nesse caso haverá duas variáveis x acessíveis ao usar x. Felizmetne o Mathematica avisa  que isso ocorre (“shadowed  variables”)

AppendTo[$ContextPath, "UpsideDown`"]

Notar que imediatamente ao rodar a linha acima todas as ocorrências de “x” ficaram vermelhas.

x

O x de UpsideDown pode ser acessado, mas devido ao “shadowing” só o pode ser feito através da escrita explícita de todo o contexto,

UpsideDown`x

Como agora UpsideDown está na $ContextPath, podemos definir símbolos em UpsideDown que todos esses serão acessíveis sem precisar especificar o contexto,

UpsideDown`y[a_]= (a^(-2));

y[a]

Context[y]

Uma forma de remover o conflito na variável x é remover UpsideDown da path, mas se isso não for desejável pode-se remover a variável UpsideDown`x ou Global`x. Em ambos os casos, um simples Clear não terá efeito, pois basta usar uma variável uma vez para ela ser atribuída a um contexto, e limpar com Clear apenas apaga seu valor, não seu contexto.

x
Clear[x]

x

Context[x]

Remove[Global`x]

x
Context[x]

Associando um contexto a um ambiente num código

Context[]

Begin["test`"]

Context[]

gg = 67;
gg
Context[gg]

End[]

gg

test`gg

Se todas as expressões definidas num certo trecho do código tiverem de ser associadas a um contexto específico, é conveniente usar o Begin[...] e End[] para demarcar essa região do código. É bom frisar que Begin e End não fazem nada que não possa ser feito escrevendo explicitamente o contexto de cada variável. Seja como for, seu uso adequado pode tanto economizar caracteres digitados como, em especial, pode deixar o código bem mais fácil de ser lido.

Context & Packages

Em essência, um pacote no Mathematica é  um conjunto de definições associadas a dado contexto. Usualmente destina-se esse nome a um sentido um pouco mais específico, espera-se normalmente que um pacote: i. esteja em um arquivo, ou conjunto de arquivos, separado de outros programas e pacotes; ii. No sentido restrito do termo no Mathematica, devem conter em seu início a instrução BeginPackage[...].

O que o BeginPackage faz é muito parecido com o Begin. A principal diferença é a forma de uso, dado que BeginPackage explicitamente demarca o início de um pacote. A segunda mais importante diferença é que o contexto criado com BeginPackage é adicionado à $ContextPath logo após o EndPackage. Assim, qualquer expressão definada (ou apenas mencionada) logo após o BeginPackage estará automaticamente disponível em qualquer contexto, sem precisar digitar explicitamente o contexto da expressão. Por fim, qualquer pacote chamado dentro de BeginPackage (certametne qualquer pacote chamado com Needs, tem de verificar para os com Get) é adicionado normalmente à $ContextPath, como usual, porém é removido da $ContextPath após o EndPackage. Há uma sintaxe própria do BeginPackage para chamar outros pacotes e não remover os contextos dos pacotes chamados depois de EndPackage, a saber:

BeginPackage["PacoteMassa`",{"Pacote1`", "Pacote2`"}] (*Não executado aqui*)

No caso acima, Pacote1 e Pacote2 serão chamados, serão adicionados à $ContextPath e lá permanecerão, mesmo após o EndPackage.

Detalhes e outras sutilezes de BeginPackage podem ser vistas aqui... https://mathematica.stackexchange.com/questions/13958/whats-the-difference-between-these-two-ways-of-specifying-dependent-packages

Vale lembrar que a estrutra padrão de criação de um pacote no Mathematica é (como já vimos na aula passada) (exemplo extraído de (HyperlinkTemplate["Setting Up Wolfram Language Packages", paclet:tutorial/SettingUpWolframLanguagePackages, paclet:tutorial/SettingUpWolframLanguagePackages, Link, HyperlinkActive, ]))

BeginPackage["Collatz`"]

Collatz::usage =
        "Collatz[n] gives a list of the iterates in the 3n+1 problem,
        starting from n. The conjecture is that this sequence always
        terminates."

Begin["`Private`"]

Collatz[1] := {1}

Collatz[n_Integer]  := Prepend[Collatz[3 n + 1], n] /; OddQ[n] && n > 0

Collatz[n_Integer] := Prepend[Collatz[n/2], n] /; EvenQ[n] && n > 0

End[ ]

EndPackage[ ]

HDInterval & observações sobre estatística no Mathematica

ETInterval:: usage = "ETInterval[dof, sigmas] outputs the equal tails interval for a ChiSquare distribution (not to be confused with the Extra Terrestrial interval). It is such that: i) when a ChiSquare distribution with dof degrees of freedom is integrated along these limits the resulting probability is Probσ[sigmas]. ii) The lower and upper bounds are such that the the exterior regions have equal probabilities.";

HDInterval:: usage = "HDInterval[dof, sigmas] outputs the highest density interval for a ChiSquare distribution (not to be confused with the High Definition interval). It is such that: i) when a ChiSquare distribution with dof degrees of freedom is integrated along the interval limits the resulting probability is Probσ[sigmas]. ii) The region is the one with highest density. Equivalently, if dof > 2, the mode is between the lower and upper bounds and these bounds have the same PDF value. If dof ≤ 2, the lower bound is zero.";


(*ETInterval definition*)
Module[{nsigma},
Probσ[nsigma_]= Integrate[PDF[NormalDistribution[0,1],x], {x, -nsigma,nsigma}];
ETInterval[dof_,nsigma_]:={Quantile[ChiSquareDistribution[dof],((1-Probσ[nsigma])/2)],
Quantile[ChiSquareDistribution[dof], ((1+Probσ[nsigma])/2) ] };
];


(*HDInterval definition*)
Module[{dof, nsigma, x, pdf, cdf, Xmax, FirstX, LevelProbability, lowerBoundCriticalX, upperBoundCriticalX, CriticalX},

Clear[LevelProbability, CriticalX];

pdf[dof_,x_]=PDF[ChiSquareDistribution[dof],x];
cdf[dof_,x_]=CDF[ChiSquareDistribution[dof],x];
Xmax[dof_]= If[dof≥ 2 , dof-2, 0]; (*The value of x that maximizes the PDF*)
 
FirstX[dof_,pdfvalue_]:= If[dof>2 ,-(-2+dof) ProductLog[-(((((2^(dof/2)) pdfvalue Gamma[(dof/2)]))^(2/(-2+dof)))/(-2+dof)) ],0]; (*Finds the value of the leftmost x for a given PDF value.*)

LevelProbability[dof_, x_?NumberQ] := cdf[dof,x] - cdf[dof,FirstX[dof,pdf[dof,x]]]  /; x>Xmax[dof];
LevelProbability[dof_, x_?NumberQ] := 0 /; x≤ Xmax[dof];
(*LevelProbability gives the PDF integral from FirstX to x (both have the same PDF values)*)

lowerBoundCriticalX[dof_,nsigma_]=Quantile[ChiSquareDistribution[dof], Probσ[nsigma] ];
upperBoundCriticalX[dof_,nsigma_]=Quantile[ChiSquareDistribution[dof], ((1+Probσ[nsigma])/2) ];

CriticalX[dof_, nsigma_]:= CriticalX[dof, nsigma] = x/.FindRoot[LevelProbability[dof, x] == Probσ[nsigma], {x, lowerBoundCriticalX[dof,nsigma], upperBoundCriticalX[dof,nsigma] }, WorkingPrecision -> 30, PrecisionGoal-> 5, AccuracyGoal->∞];
(*CriticalX find the value of the rightmost x such that LevelProbability has the probability inferred from nsgima.*)

HDInterval[dof_, nsigma_] := {FirstX[dof,pdf[dof,CriticalX[dof,nsigma]]],  CriticalX[dof,nsigma]};
];

Testing HDInterval and ETInterval

Module[{ν, σ, rangeHD, rangeET, pdf},
ν=4;
σ = 2;
rangeHD = HDInterval[ν,σ];
rangeET = N@ETInterval[ν,σ];
pdf[ν_,x_] := PDF[ChiSquareDistribution[ν],x];
Show[Plot[pdf[ν,x], {x,0,1.1 rangeET[[2]]}],Graphics[{Line[{{rangeHD[[1]],pdf[ν, rangeHD[[2]]]},{rangeHD[[2]],pdf[ν, rangeHD[[2]]]}}], Blue,Line[{{rangeET[[1]],0},{rangeET[[1]],1} }], Line[{{rangeET[[2]],0},{rangeET[[2]],1} }]}]]
]