Métodos computacionais / 2019

Davi C. Rodrigues
Aula 10

e-mail: davi@cosmo-ufes.org
página do curso: davi.cosmo-ufes.org/mc2019

item 1: revisão da última aula

Tópico abaixo não segue o livro

Sobre  o uso e métodos de resolução de problemas com NIntegrate, NDSolve, FindMinimum e NMinimize

NIntegrate / NDSolve

Do[...., {i,1,5}, {j,4,7}]

Comecemos por um problema simples, que admite solução analítica.

ρ[r_] = ((ρs)/((r/rs)(((1 + (r/rs)))^2)));

 Integrate[((ρs)/((r/rs)(((1 + (r/rs)))^2)))  (r^2) Sin[θ], {r,0,R}, {θ,0,π}, {ϕ,0,2 π}]

4 π Integrate[((ρs)/((r/rs)(((1 + (r/rs)))^2)))  (r^2) , {r,0,R}]

Simplify[ConditionalExpression[4 π (rs^3) ρs (-1+(rs/(R+rs))-Log[rs]+Log[R+rs]),((Re[(rs/R)]≥0&&(rs/R)≠0)||(rs/R)∉(Reals)||Re[(rs/R)]<-1)&&((Im[R]≥0&&Im[rs]≥0)||(Im[R]≤0&&Im[rs]≤0)||((Im[rs])/(Im[R]))≤-1||((Im[rs] Re[R])/(Im[R]))≤Re[rs])], Assumptions -> {R>0, rs>0}]

M[R_] =4 π Integrate[((ρs)/((r/rs)(((1 + (r/rs)))^2)))  (r^2) , {r,0,R}, Assumptions -> {R>0, rs >0}]

M[0]
M[2]

Este é um problema simples que a solução analítica é também simples, mas usaremos este exemplo para ilustrar o emprego de métodos numéricos. E aqui teremos a vantagem de ser fácil saber se a solução numérica é boa ou não, pois já temos a solução analítica.

O primeiro passo a notar é que não se pode fazer de imediato a integração numérica... Relacionado a isso,  “Assumptions” não é compatível com NIntegrate.

NIntegrate[((ρs)/((r/rs)(((1 + (r/rs)))^2)))  (r^2) , {r,0,R} ]
NIntegrate[((ρs)/((r/rs)(((1 + (r/rs)))^2)))  (r^2) , {r,0,10}]

Nota-se que é fácil eliminar as constantes ρs e rs. 

- ρs é apenas uma constante multiplicativa, que pode ser facilmente inseria fora de NIntegrate.
- r pode ser reescalado, eliminando a dependência explícita em rs dentro da integral. 

Assim, temos

MN[Rs_] := 4 π ρs (rs^3) NIntegrate[(1/(r(((1 + r))^2))) (r^2) , {r,0,Rs}]

Usamos Rs só para servir de lembraça de que a entrada de MN não será examenete a mesma de M, pois mudamos os limites de integração, r -> r rs => R -> Rs = R rs

MN[0]

MN[2] /. rs-> 1
(M[2] /. rs -> 1)// N

MN[2] /. rs-> 2
(M[2] /. rs -> 2)// N

MNp[R_, rs_] := 4 π ρs (rs^3) NIntegrate[(1/(r(((1 + r))^2))) (r^2) , {r,0,(R/rs)}]

MNp[2,2]
M[2] /. rs -> 2 // N

Verificando o erro relativo, com rs=1, sem considerar a solução exata:

tab =Table[{r,MN[r] }, {r, 0, 50, 0.5}]/. {rs -> 1, ρs -> 1};

tab[[4]]

ListPlot[tab]

Clear[MNI]

MNI = Interpolation[tab, Method -> Spline]

MNI[2]

Plot[MNI[r], {r,0,50}, PlotStyle->Red]

Graphics[Circle[], ImageSize->10]

Show[ListPlot[tab, PlotMarkers-> Graphics[Circle[], ImageSize->5]], Plot[MNI[r], {r,0,50}, PlotStyle->Red]]

Note que a versão interpolada pode ser derivada

MN'[2]

MNI'[2]

Plot[{MNI'[r], ((4 π)/(r(((1 + r))^2))) (r^2) }, {r,0,50}, PlotRange->All]

Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI'[r]}, {r,0,50}, PlotRange->All]
Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI'[r]}, {r,0,50}]

Qual seria a fonte do erro, a interpolação ou a integração numérica?
Para verificar se é erro da interpolação, podemos facilmente aumentar o número de pontos...

tab =ParallelTable[{r,MN[r] }, {r, 0, 50, 0.01}]/. {rs -> 1, ρs -> 1};
MNI = Interpolation[tab, Method -> Spline, InterpolationOrder->3];
Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI'[r]}, {r,0,50}, PlotRange->All]
Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI'[r]}, {r,0,50}]

MNI = Interpolation[tab, InterpolationOrder->5, Method -> Spline];
Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI'[r]}, {r,0.01,50}, PlotRange->All]
Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI'[r]}, {r,0.01,50}]
Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI'[r]}, {r,0.01,0.1}, PlotRange-> All]

Note que há uma sistemática, mas note também que são erros ainda assim muito pequenos.

Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI'[r]}, {r,0.01,0.1}, PlotRange-> All]

Agora comparemos com a solução exata.

Block[{rs=1, ρs=1},
Plot[{1 - M[r]/MN[r] }, {r,0.001,50}, WorkingPrecision->50]]
Block[{rs=1, ρs=1},
Plot[{1 - M[r]/MN[r] }, {r,0.001,50}, WorkingPrecision->50, PlotRange->All]]

E se desejássemos chegar até uma precisão não de ~13, mas de 20 casas?

Clear[MNp20]
MNp20[R_, rs_] := 4 π ρs (rs^3) NIntegrate[(1/(r(((1 + r))^2))) (r^2) , {r,0,(R/rs)}, WorkingPrecision->50, PrecisionGoal-> 20, AccuracyGoal->Infinity]
MNp20[r_] := MNp20[r,1];

MNp20[2,2]
MNp[2,2]

Block[{rs=1, ρs=1},
Plot[{1 - M[r]/MNp20[r,1] }, {r,0.001,50}, WorkingPrecision->50]]
Block[{rs=1, ρs=1},
Plot[{1 - M[r]/MNp20[r,1] }, {r,0.001,50}, WorkingPrecision->50, PlotRange->All]]

Block[{rs=1, ρs=1},
Plot[{1 - M[r]/MNp20[r,1] }, {r,0.001,50}, WorkingPrecision->80, PlotRange->All, MaxRecursion-> 15, PlotPoints-> 1000]]

Pedimos 20 de precisão, e ganhamos mais do que foi solicitado... Isso pode acontecer, mas em problemas complexos ocorre comumente o contrário!

Vale lembrar que o último teste acima comumente não pode ser feito, pois a solição exata não é conhecida.

MNp20[2]

tab20 =ParallelTable[{r,MNp20[r] }, {r, 0, 50, 0.01}]/. {rs -> 1, ρs -> 1};
MNI20 = Interpolation[tab20, Method -> Spline];
Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI20'[r]}, {r,0.01,50}, PlotRange->All]
Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI20'[r]}, {r,0.01,50}]
Plot[{1 - ( ((4 π)/(r(((1 + r))^2))) (r^2) )/MNI20'[r]}, {r,0.01,0.1}, PlotRange-> All]

Observacoes

- Em alguns casos, indicar para o Mathematica ponto(s) para serem excluídos da integração pode tanto melhorar precisão quanto performance.

- Eliminar “SymbolicProcessing” quando qualquer manipulação simbólica é desnecessária, acelera a computação.

MNp20F[R_, rs_] := 4 π ρs (rs^3) NIntegrate[(1/(r(((1 + r))^2))) (r^2) , {r,0,(R/rs)}, WorkingPrecision->50, PrecisionGoal-> 20, AccuracyGoal->Infinity, Method -> {Automatic, "SymbolicProcessing" -> False}]
MNp20F2[R_, rs_] := 4 π ρs (rs^3) NIntegrate[(1/(r(((1 + r))^2))) (r^2) , {r,0,(R/rs)}, WorkingPrecision->50, PrecisionGoal-> 20, AccuracyGoal->Infinity, Method -> {Automatic, "SymbolicProcessing" -> False}, Exclusions-> 0]

MNp20F2[34.2345, 2] // Timing
MNp20F[34.2345, 2] // Timing
MNp20[34.2345, 2] // Timing

Exemplo de caso mais complicado e sem solução exata

ρ[r_]= ((10^logRhoc)/((1 + (r/rc)) (1 + ((r^2)/(rc^2)))));

int[r_, θ_] =- (( (E^(-(( 1)/λ)(Sqrt[(r^2)+ (y^2) - 2 r y Cos[θ] ]))))/((Sqrt[(r^2)+ (y^2) - 2 r y Cos[θ] ]) ))  ρ[r]Sin[θ] (r^2);

The effective mass, apart from multiplicative factor on β , depends on r, rc and λ . But, apart from an overall λ multiplication, r and rc only appear dividing λ; thus we can consider only the ratios r/λ and rc/λ.

The intervls that the table below will explore are:

(10^(-4.3)) < r/λ < (10^3.3)
(10^(-4.3)) < rc/λ < (10^3.3)

Thus,  the following bounds satisfy those above:

0.01< r < 100
0.01 < rc < 200
0.1< λ < 200

Since 0.01/200 = (10^(-4.3)) and 200/0.1 = (10^3.3).

Other variation also compatible with this table, which will be used in version GY-9:

0.0503< r < 199.5
0.0503 < rc <199.5
0.1< λ < 1000

Plot[meffreal[r,1,1],{r,20,30}]

meffrealExact0[r_,rs_,l_] := 
With[{R = Rationalize[r,0], RS = Rationalize[rs,0], L = Rationalize[l,0]}, 
meffreal[R,RS,L]];

----

SetDirectory[NotebookDirectory[]]

aux[r_,rs_]=Exp[-r]( (E^(2 r)) (-1+r) (-π (Cos[rs]+Sin[rs] )- (E^( rs)) ExpIntegralEi[-r-rs]+
(Re[(E^( I rs))ExpIntegralEi[-r-I rs]]-Im[(E^( I rs))ExpIntegralEi[-r-I rs]]))+
(1+r) (- (E^( rs)) ExpIntegralEi[-rs]-π (Cos[rs]+Sin[rs] )+
Re[(E^( -I rs))ExpIntegralEi[r+I rs]]-Im[(E^( -I rs))ExpIntegralEi[r+I rs]]-2   Im[(E^( I rs))  ExpIntegralEi[-I rs]]+  (E^( -rs))  ExpIntegralEi[rs]- (E^( -rs))ExpIntegralEi[r+rs]));
meffreal[r_,rs_,l_]=Simplify[aux[r/l,rs/l]*l];
meffrealExact0[r_,rs_,l_] := 
With[{R = Rationalize[r,0], RS = Rationalize[rs,0], L = Rationalize[l,0]}, 
meffreal[R,RS,L]];

$MaxExtraPrecision=5000;
LaunchKernels[];
ParallelEvaluate[$MaxExtraPrecision=5000];

tab = ParallelTable[{{(10^(-4.3 + 0.025 i)), (10^(-4.3 + 0.025 j))}, N[meffrealExact0[(10^(-4.3 + 0.025 i)), (10^(-4.3 + 0.025 j)),1], 10]}, {i,0,304}, {j,0,304}];

Export["tableMeffBurkertYukawa2.dat", tab, "Real32"];

NMinimize

tutorial/ConstrainedOptimizationGlobalNumerical#469184158

Sobre o uso de pacotes

Várias aplicações foram apresentadas. Em especial com respeito ao xAct e ao MAGMA