Resolución numérica de ecuacións de Fokker-Planck nunha dimensión espacial

T aballo Fin de G ao
Resolución numé ica de ecuacións de
Fokke -Planck nunha dimensión
espacial
Ja ie López Ped a es
2019/2020
UNIVERSIDADE DE SANTIAGO DE COMPOSTELA
GRAO DE MATEMÁTICAS
T aballo Fin de G ao
Resolución numé ica de ecuacións de
Fokke -Planck nunha dimensión
espacial
Ja ie López Ped a es
Xullo 2020
UNIVERSIDADE DE SANTIAGO DE COMPOSTELA
T aballo p opos o
Á ea de Coñecemen o: Ma emá ica Aplicada
Tí ulo: Resolución numé ica de ecuacións de Fokke -Planck
nunha dimensión espacial
B e e desc ición do con ido
Resol e anse algúns modelos de anspo e de pa ículas ca -
gadas, como elec óns ou ións pesados. Os mé odos numé i-
cos es a án baseados en mé odos de di e enzas ni as.
Ti o : Ósca López Pouso
iii

Índice xe al
Resumo
ii
In odución
ix
1. Esquema numé ico 1
2. In e p e ación ísica 9
3. Implemen ación do esquema 13
4. Resul ados numé icos 21
Anexos 31
Anexo I: Fó mulas de de i ación numé ica 33
Anexo II: Códigos Ma lab 37
Bibliog a ía 43
Resumo
O aballo comeza in oducindo un coñecido p oblema no mundo da ísica, o cal o mu-
lamos ma ema icamen e. Pos e io men e con inuamos desc ibindo un esquema numé ico
que pe mi e ob e unha solución pa a o p oblema baseado nas di e enzas ni as.
Ademais dunha o mulación ma emá ica apó ase unha isión e unha explicación ísi-
ca do p oblema. Finalmen e, desc ibimos o algo i mo i e a i o que esol e o p oblema e
implemen ámolo en Ma lab pa a ob e esul ados numé icos.
Abs ac
The a icle s a s in oducing a well-known physics p oblem, ha we go o o mula-
e ma hema ically. Then, we con inue desc ibing a nume ic scheme ha allows ob ain a
solu ion o ou p oblem based on ni e die ences.
In addi ion o a ma hema ical o mula ion we show a li le sigh and a physical expla-
na ion o he p oblem. A he end, we nish wi h he i e a i e algo i hm's desc ip ion ha
sol es he p oblem and we implemen i on Ma lab o ob ain nume ical esul s.
ii
4
CAPÍTULO 1. ESQUEMA NUMÉRICO
Así acabamos de disc e iza unha p imei a pa e do noso p oblema, ago a p ecisamos
disc e iza o e mo da di usi idade. Dependendo do pun o da malla onde nos a opemos
es e e mo pode se nulo (cando
µ=±1
); en ón conside amos os seguin es casos:
(1) Se nos a opamos en
i= 1
, ou sexa
µ1=−1
emos que
D(µ1) = 1−µ2
1= 1−(−1)2=
1−1=0
, i.e.,
D1= 0
. Nes e caso podemos ealiza a seguin e disc e ización:
∂
∂µ D(µ)∂ψ
∂µ (µ,z)=(µ1,zm)
O(h2)
≈4D2∂ψ
∂µ (µ2, zm)−D3∂ψ
∂µ (µ3, zm)
2h
(1.7)
e cando nos a opamos con
∈ {2,3}
, emp egamos a ó mula es ánda cen ada en
dous pun os pa a disc e iza as de i adas de p imei a o de que apa ecen na ecuación
an e io :
∂ψ
∂µ (µ , zm)
O(h2)
≈ψm
+1 −ψm
−1
2h.
(1.8)
(2) Aquí analizamos os nodos in e io es da malla, ou sexa,
i∈ {2, ... , I −1}
. Aquí bas a
emp ega a ó mula es ánda de segunda o de cen ada:
∂
∂µ D(µ)∂ψ
∂µ (µ,z)=(µi,zm)
O(h2)
≈
Di−1
2ψm
i−1−Di−1
2+Di+1
2ψm
i+Di+1
2ψm
i+1
h2.
(1.9)
(3) Queda es uda a disc e ización no nodo nal, ou sexa en
µI= 1
. A opamos de no o
unha di usión nula, logo p ocedemos analogamen e que no p imei o caso:
∂
∂µ D(µ)∂ψ
∂µ (µ,z)=(µI,zm)
O(h2)
≈DI−2∂ψ
∂µ (µI−2, zm)−4DI−1∂ψ
∂µ (µI−1, zm)
2h
(1.10)
e de no o pa a
∈ {I−2, I −1}
p océdese igual que en (1.8).
Ago a que xa emos disc e izadas as de i adas do noso p oblema, podemos p ocede a
desc ibi o esquema numé ico comple o que nos le a a a opa a solución
ψ
.
Te emos que p es a especial a ención ao núme o de ecuacións esul an es, pois de-
pendendo da pa idade de
I
a opa emos que o núme o de ecuacións non coincide co de
incógni as. En is a ás o mulas de disc e ización sinaladas an es p ecisamos esixi alome-
nos que
I≥4
e
N > 2
. P ocedamos así a desc ibi o esquema numé ico de no o sepa ando
en casos como xemos pa a disc e iza :

5
- Pa a
(i, n)∈ {1}×{1, ... , N −1}
,
−µ1
k+αn
1
2+σn
1D2
2h2ψn
1+−σn
1D3
8h2ψn
2+−σn
1D2
2h2ψn
3+
+σn
1D3
8h2ψn
4+µ1
k+αn+1
1
2+σn+1
1D2
2h2ψn+1
1+
+−σn+1
1D3
8h2ψn+1
2+−σn+1
1D2
2h2ψn+1
3+
+σn+1
1D3
8h2ψn+1
4=Wn
1+Wn+1
1
2.
(1.11)
- Pa a
(i, n)∈ {2, ... , I −1}×{1, ... , N −1}
,
−
σn
iDi−1
2
2h2!ψn
i−1+
+
−µi
k+αn
i
2+
σn
iDi−1
2+Di+1
2
2h2
ψn
i+
+ −
σn
iDi+1
2
2h2!ψn
i+1 + −
σn+1
iDi−1
2
2h2!ψn+1
i−1+
+

µi
k+αn+1
i
2+
σn+1
iDi−1
2+Di+1
2
2h2
ψn+1
i+
+ −
σn+1
iDi+1
2
2h2!ψn+1
i+1 =Wn
i+Wn+1
i
2.
(1.12)
- Pa a
(i, n)∈ {I}×{1, ... , N −1}
,
σn
IDI−2
8h2ψn
I−3+−σn
IDI−1
2h2ψn
I−2+−σn
IDI−2
8h2ψn
I−1+
+−µI
k+αn
I
2+σn
IDI−1
2h2ψn
I+ σn+1
IDI−2
8h2!ψn+1
I−3+
+ −σn+1
IDI−1
2h2!ψn+1
I−2+ −σn+1
IDI−2
8h2!ψn+1
I−1+
+ µI
k+αn+1
I
2+σn+1
IDI−1
2h2!ψn+1
I=Wn
I+Wn+1
I
2.
(1.13)
Temos así comple amen e desc i o o esquema numé ico que nos pe mi e esol e o
p oblema denido polas ecuacións (1), (2) e (3). Pe o a opámonos cun p oblema pois nun
de e minado caso o núme o de ecuacións non coincide co de incógni as.
6
CAPÍTULO 1. ESQUEMA NUMÉRICO
Pa a es uda di a p oblemá ica bas a conside a un
I
impa . O núme o de incógni as
que emos é
I×N
, pe o pola con a o esquema que acabamos de desc ibi unicamen e posúe
I×(N−1)
ecuacións. Ago a ben emos que comp oba se as ecuacións (condicións de uxo
en an e) (2) e (3) p opo cionan ecuacións sucien es pa a sol en a es e desaxus e. Así
pois depende á da pa idade de
I
o núme o de ecuacións p opo cionadas. Cando
I
é pa
po po cionan
I
ecuacións, pe o cando é
I
impa unicamen e p opo cionan
I−1
ecuacións.
Así emos que cando
I
é impa al a unha ecuación e non podemos esol e o esquema
p opos o an e io men e.
Se
I
é pa esol emos o p oblema emp egando o esquema al como oi desc i o. A
con inuación desc ibimos como comple a o esquema pa a
I
impa .
Sexa
I
impa , logo é ob io que
µI+1
2= 0
simplemen e po cons ución. En p imei a
ins ancia podemos pensa en emp ega unha das condicións numé icas inicial ou nal:
ψ1
I+1
2
= (0)
ou
ψN
I+1
2
=g(0),
(1.14)
pe o en ealidade a in e ese adica en impoñe as dúas condicións debido a que es amos
in e esados en ob e solucións con inuas.
En ón o que acemos é elimina as
N−1
ecuacións co esponden es a
i=I+1
2
pa a
impoñe as condicións numé icas an e io es e así conside a un no o conxun o de
N−2
ecuacións pa a
µ= 0
e os nodos do in e alo
(Z
ini
, Z
n
)
. En ón pa a
µ= 0
a ecuación (1)
eesc íbese como:
α(0, z)ψ(0, z)−σ(0, z)∂2ψ
∂µ2(0, z) = W(0, z)
pa a
z∈[Z
ini
, Z
n
].
(1.15)
Vexamos que e ec i amen e se e ica o an e io . Pa a subs i uí en
µ= 0
p ecisamos
an es de i a o e mo que le a a di usión emp egando a eg a de de i ación do p odu o e
pos e io men e e alua en
µ= 0
. De i emos en ón:
∂
∂µ (1 −µ2)∂ψ
∂µ =−2µ∂ψ
∂µ + (1 −µ2)∂
∂µ ∂ψ
∂µ =−2µ∂ψ
∂µ + (1 −µ2)∂2ψ
∂µ2.
(1.16)
Ago a se subs i uímos pa a
µ= 0
emos:
∂
∂µ (1 −µ2)∂ψ
∂µ (µ,z)=(0,z)
=−2·0∂ψ
∂µ (0, z) + (1 −02)∂2ψ
∂µ2(0, z) = ∂2ψ
∂µ2(0, z),
(1.17)
que e ec i amen e é o que emp egamos en (1.15).
Se pensamos no in e alo abe o
(Z
ini
, Z
n
)
a ecuación (1.15) suxi e:
αn
i∗ψn
i∗−σn
i∗
ψn
i∗−1−2ψn
i∗+ψn
i∗+1
h2=Wn
i∗
pa a
n∈ {2, ... , N −1},
(1.18)
7
sendo
i∗=I+1
2
. O an e io pode eesc ibi se do seguin e xei o:
−σn
i∗
h2ψn
i∗−1+αn
i∗+2σn
i∗
h2ψn
i∗+−σn
i∗
h2ψn
i∗+1 =Wn
i∗
(1.19)
pa a
n∈ {2, ... , N −1}
.
En ón o esquema nal pa a o caso
I
impa consis e no desc i o es ándolle as
N−1
ecuacións co esponden es a
i=I+1
2
, engandindo as condicións de con o no numé icas
(1.14) e as
N−2
ecuacións (1.19) que acabamos de ob e .
Como comen a emos máis adian e implemen a emos o esquema pa a o caso
I
impa
pois ob éñense mello es esul ados.
8
CAPÍTULO 1. ESQUEMA NUMÉRICO
Capí ulo 2
In e p e ación ísica
A ecuación de Fokke -Planck (EFP) en unha i e esan e in e p e ación ísica. Imos
mos a a elación do noso p oblema coa ecuación do anspo e de Bolzmann (ETB).
A impo ancia do p oblema de Fokke -Planck adica pois en que a solución des e,
ψ
,
é unha ap oximación da solución do p oblema de Bolzmann en ce as condicións. O in e-
esan e é es uda di as condicións. Así pois a EFP é unha ap oximación da de Bolzmann
cando alamos do anspo e de pa ículas que so en pequenas des iacións na súa axec-
o ia ao impac a con ou as pa ículas e pequenas pe das de ca ga. Nes e caso alamos de
dispe sión ou sca e ing. Exemplos comúns de pa ículas que so en di os compo amen-
os son pa ículas ca gadas como ións pesados ou elec óns da codia do á omo que posúen
ca ga nega i a.
Ámbalas dúas ecuacións an e io es desc iben a densidade de uxo angula de pa ículas,
ψ
. Ademais es ablecen o equilib io en e as pe das e as ganancias de
ψ
ao longo das
di eccións de p opagación. Di o dou o xei o es udan o g adien e de
ψ
ao longo de cada
di ección do espazo idimensional
R3
.
Conside emos
Ω⊂R3
un dominio espacial e omemos
ω
como di ección de p opagación.
Se esc ibimos
ω
en coo deadas es é icas emos:
ω=ω(ϕ, θ) = (sen ϕcos θ, sen ϕsen θ, cos ϕ)∈
S
2,
sendo
S
2
a es e a uni a ia do espazo con
ϕ∈[0, π]
o ángulo pola e
θ∈[0,2π)
o ángulo
acimu al.
Se o uxo depende da posición
x= (x1, x2, x3)∈Ω
unicamen e a a és de
x3=z
,
emos que
ω· ∇ψ=ω3∂ψ
∂z
. Se na ecuación (1) conside amos
µ=ω3
ob emos exac amen e
o an e io , sendo
ω
a di ección de p opagación. No noso caso es amos aballando sob e
a ecuación unidimensional, en ón o noso dominio se á da o ma
Ω = R2×(Z
ini
, Z
n
)
, o
9

10
CAPÍTULO 2. INTERPRETACIÓN FÍSICA
Figu a 2.1: Sección do slab unidimensional.
cal adoi a chama se slab unidimensional. Unha pequena ilus ación do an e io é a gu a
2.1.
Ago a ben, pa imos de es a iables espaciais e imos educilas a unicamen e unha.
Is o pódese ace g azas a que o uxo,
ψ
, non a ía nin con
x1
nin con
x2
e ademais o
dominio an e io son copias do dominio unidimensional
(Z
ini
, Z
n
)
.
Se con inuamos co desc i o emos
ω3=µ= cos ϕ
, con is o obse amos i ialmen e a
azón de po que o pa áme o
µ
a ía en e
1
e
−1
pois
µ
é o coseno dun ángulo.
Ao pensa na EFP unidimensional deixamos de depende do ángulo acimu al, en ón
unicamen e emos o ángulo pola ou equi alen emen e
µ
e así o ope ado
∂
∂µ (1 −µ2)∂·
∂µ
(2.1)
coñecido comunmen e como o ope ado con inuo de dispe sión, o cal non é máis que o
laplaciano sob e a es e a.
Podemos dici así que a ecuación (1) é a EFP monoene xé ica supoñendo que a ene xía
é cons an e no dominio unidimensional con sime ía de xeome ía plana e as ecuacións (2)
e (3) impoñen a densidade de uxo de pa ículas que en a no dominio. Po is o úl imo
di as ecuacións eñen o nome de condicións de uxo en an e.
Como mencionamos ao nal da In odución o noso p oblema en unha ca ac e ís ica
peculia den o das ecuacións en de i adas pa ciais: a ausencia de condicións de con o no
pa a
|µ|= 1
, ou sexa o ope ado (2.1) non p opo ciona di as condicións. Is o débese a que
o equilib io dado pola EFP é sucien e pa a ob e o uxo angula unha ez coñecemos o
uxo a a és da on ei a ísica. En ón impoñendo unhas condicións de con o no o que
11
a íamos se ía sob ede e mina o p oblema.
Finalmen e poñamos un signicado ísico a cada un dos e mos que apa ecen na ecua-
ción (1):
- Po
ψ(µ, z)
deno amos á densidade de uxo angula de pa ículas en
(µ, z)
. Di o
dou o xei o é o núme o de pa ículas que se mo en desde
z
ao longo da di ección
ω(µ)
po unidade de á ea no mal a
ω(µ)
, po unidade de ángulo sólido
1
, po unidade
de ene xía e po unidade de empo.
- Cando alamos de mo e nos dende
z
na di ección de
ω(µ)
e e ímonos cla amen e a
mo e nos en
(x1, x2, z)
pa a calque a
(x1, x2)∈R2
xos e ao longo da di ección
ω(µ) = p1−µ2cos θ, p1−µ2sen θ, µ
pa a calque a ángulo
θ∈[0,2π)
xo.
- O p imei o e mo da ecuación,
µ∂ψ
∂z
, é a de i ada di eccional de
ψ(µ, z)
ao longo da
di ección
ω(µ)
.
- Deno amos po
α
ao coecien e de abso ción, así pois
αψ
desc ibe as pe das po
abso ción.
- O e mo
σ∂
∂µ h(1 −µ2)∂ψ
∂µ i
é o e mo de di usión na a iable angula .
- A unción de dúas a iables
W
ep esen a unha on e in e na de pa ículas onde eña
sigmo posi i o, e un sumidoi o in e no onde o signo sexa nega i o.
1
Po ángulo sólido en endemos un ángulo espacial que aba ca un obxe o cando é is o dende un pun o.
12
CAPÍTULO 2. INTERPRETACIÓN FÍSICA
Capí ulo 3
Implemen ación do esquema
Como mencionamos no Capí ulo 1 ob i emos dous esquemas dis in os segundo a pa i-
dade de
I
. T as a ealización de di e sos expe imen os numé icos chegouse á conclusión de
que e a máis ecomendable a implan ación do esquema pa a
I
impa .
T as os di e sos expe imen os ealizados obse ouse que pa a o caso
I
impa non apa-
ecían ines abilidades a edo de
µ= 0
. Po én, no caso
I
pa cando nos achegamos a
Q0
poden apa ece solucións que son ines ables, p esen ando así oscilacións espu ias.
En e mos de implemen ación do esquema pode íamos conside a dous ipos de esque-
ma:
- Mé odo di ec o: un esquema no que esol emos a ecuación en odo o dominio,
Q
.
Conside amos unha g an ma iz cad ada de dimensión
I×N
e en amos en esol e
o sis ema linea comple o asociado ao p oblema.
- Mé odo i e a i o: algo i mo no que se esol e o p oblema pa indo dunha semen e
inicial en
Q0
, esol emos po ascenso e descenso, espec i amen e en
Q+
e
Q−
, os
p oblemas de alo inicial e nal denidos polas uncións
e
g
. Realizamos es e
p oceso a a que os alo es ob idos en dúas i e acións sucesi as pa a a solución en
Q0
dis en unha can idade eal posi i a épsilon moi pequena. Ou amén se non se acada
con e xencia dado un núme o de i e acións máximas.
Nes e aballo cen a émonos na implemen ación do algo i mo i e a i o pa a o caso
I
impa .
Conside emos en ón de no o os conxun os do dominio
Q
desc i os ao comezo do Ca-
pí ulo 1. Sexan en ón
Q+= (0,1] ×[Z
ini
, Z
n
]
, o segmen o
Q0={0} × [Z
ini
, Z
n
]
e
Q−= [−1,0) ×[Z
ini
, Z
n
]
( e gu a 3.1).
13
20
CAPÍTULO 3. IMPLEMENTACIÓN DO ESQUEMA

Capí ulo 4
Resul ados numé icos
Nes a sección mos a emos unha se ie de esul ados numé icos ob idos a pa i es do
algo i mo i e a i o que implemen amos en Ma lab.
Pa a odos os expe imen os conside a emos
Z
ini
= 0
,
Z
n
= 1
,
i max = 2000
(núme o
de i e acións máximas) e oma amos a seguin e ole ancia,
ε= 10−8
.
Ademais, agás que así se indique, en odos os es s ealizados emp égase como semen e
pa a inicializa o algo i mo a dada pola ecuación (3.8).
En can o no ación, po
E
abs
(Q)
en endemos o seguin e:
E
abs
(Q) = m´ax
Q|ψmalla −ψ|,
sendo
ψmalla
a solución ap oximada e o máximo ómase sob e o conxun o de odos os nodos
da malla.
As ó mulas que o on emp egadas pa a a ob ención do esquema numé ico son exac as
cando a solución exac a é polinómica de g ao meno ou igual que 1 na a iable
µ
e de g ao
meno ou igual que 2 pa a a a iable
z
.
- TEST 1.
Supoñamos que a solución exac a é coñecida, en ón sexa:
ψ(µ, z) = µ2z2.
Logo ealizando uns simples cálculos ob emos as uncións
,
g
e
W
subs i uíndo nas
21
22
CAPÍTULO 4. RESULTADOS NUMÉRICOS
ecuacións (1)-(3):
(µ) = ψ(µ, Z
ini
) = ψ(µ, 0) = 0,
g(µ) = ψ(µ, Z
n
) = ψ(µ, 1) = µ2,
W(µ, z) = µ∂ψ
∂z +αψ −σ∂
∂µ (1 −µ2)∂ψ
∂µ =
= 2µz3+α(µ, z)µ2z2−σ(µ, z)z2(2 −6µ2),
conside ando:
α(µ, z) = |sen(12µz)|
e
σ(µ, z) = 1 + sen(12µz) cos(12µz).
Como podemos obse a no cad o 4.1 apa ece a e minoloxía
(µ, z)
, con is o e e í-
monos ao pun o da malla
(µi, zn)
no que se a opa o e o máximo en alo absolu o.
Toma emos di a no ación pa a os pos e io es es e cad os.
Ademais amén se in en a eexa a o de de con e xencia do esquema pa a iso ealí-
zase o seguin e:
Sexa
E(h, k)
o e o máximo come ido na ap oximación pa a unha malla de amaño
h
-
k
e sexa
c
unha cons an e, logo se a o de é
p
an o pa a
h
coma pa a
k
debe
e ica se asin o icamen e (
h, k →0
) a seguin e igualdade:
E(ch, ck) = cpE(h, k).
Tomando loga i mos podemos despexa a o de, ense así:
p=
ln E(ch,ck)
E(h,k)
ln c.
(I, N)E
abs
(Q) (µ, z)
I e acións O de
(11,10)
1.26187
×10−2
(-1,0.888...) 32
(33,10)
7.54505
×10−4
(-1,0.888...) 88 2.2421760954
(101,10)
7.82894
×10−5
(0.8,1) 183 1.988398736
(321,10)
8.08662
×10−6
(0.85625,1) 256 1.951768725
(1001,10)
4.98727
×10−7
(0.914,1) 1 2.444991618
(1001,901)
8.50361
×10−7
(0.814,1) 504
Cad o 4.1: Resul ados numé icos pa a o es 1 conside ando
ω= 2
.
23
Po exemplo se conside amos
h
2
e
k
2
, emos a seguin e cons an e
c=1
2
.
Pa a es e exemplo obse amos que
ψ
é polinómica de g ao 2 na a iable
z
. Pe o
non é polinómica de g ao meno ou igual que 1 na a iable
µ
, en ón pa a ob e
mello es esul ados o que acemos é ena a malla en di a a iable, ou sexa en
h
.
Di o compo amen o eíc ese amén no cad o 4.1.
- TEST 2.
Supoñamos de no o que coñecemos a solución exac a ao p oblema, sexa así:
ψ(µ, z) = µz3.
Realizando os cálculos pe inen es ob emos:
(µ) = ψ(µ, Z
ini
) = ψ(µ, 0) = 0,
g(µ) = ψ(µ, Z
n
) = ψ(µ, 1) = µ,
W(µ, z) = µ∂ψ
∂z +αψ −σ∂
∂µ (1 −µ2)∂ψ
∂µ =
= 3µ2z2+µz3(α(µ, z)−2σ(µ, z)) .
Como podemos obse a nos cad os 4.2 e 4.3 apa ece o e mo non con e xencia, con
is o e e ímonos a que non se acada a con e xencia do mé odo chegado o núme o de
i e acións máximas. Po én, o algo i mo es á a compo a se de manei a con e xen e
pois o esiduo calculado as a i e ación nal achégase á ole ancia.
Po exemplo pa a o cad o 4.2 o úl imo esiduo calculado, onde apa ece a non con e -
xencia, é
es =
1.18673
×10−7
, que como obse amos achégase á ole ancia épsilon.
(I, N)E
abs
(Q) (µ, z)
I e acións O de
(11,10)
2.81389
×10−3
(-1,0) 33
(11,29)
2.98535
×10−4
(-1,0) 405 1.976630514
(11,91)
2.90146
×10−5
(-1,0) Non con e xencia 1.996469446
(11,281)
3.00019
×10−6
(-1,0) 195 1.999263129
(11,901)
2.89858
×10−7
(-1,0) 15 2.001566823
(1001,901)
2.21538
×10−7
(-1,0) 102
Cad o 4.2: Resul ados numé icos pa a o es 2 conside ando
ω= 2
e omando de no o
α(µ, z) = |sen(12µz)|
e
σ(µ, z) = 1 + sen(12µz) cos(12µz)
.
24
CAPÍTULO 4. RESULTADOS NUMÉRICOS
Figu a 4.1: Solución ap oximada ob ida e e o come ido pa a o es 2 no caso
(I, N) =
(1001,901)
.
Nes e segundo es obse amos que a solución exac a
ψ
é polinómica de g ao 1 na
a iable
µ
, po én non é de g ao meno ou igual que 2 na a iable
z
, logo aquí é de
in e ese ena a malla en
z
, ou sexa en
k
, como acemos no cad o 4.2.
- TEST 3.
Supoñamos que a solución exac a é coñecida, en ón sexa:
ψ(µ, z) = ln(2 + µ2+z3).
Logo ealizando uns sinxelos cálculos ob emos:
(µ) = ψ(µ, Z
ini
) = ψ(µ, 0) = ln(2 + µ2),
g(µ) = ψ(µ, Z
n
) = ψ(µ, 1) = ln(3 + µ2),
W(µ, z) = µ∂ψ
∂z +αψ −σ∂
∂µ (1 −µ2)∂ψ
∂µ =3µz2
2 + µ2+z3+α(µ, z) ln(2 + µ2+z3)−
−σ(µ, z)−4µ2
2 + µ2+z3−2(1 −µ2)µ2−z3−2
(2 + µ2+z3)2.
De no o no cad o 4.3 apa ece a non con e xencia debido a que non se sa is ai o
es de pa ada chegado o núme o máximo de i e acións. Nes e caso o úl imo esiduo
calculado é
es =
7.17688
×10−7
25
(I, N)E
abs
(Q) (µ, z)
I e acións
(11,10)
7.05208
×10−3
(-1,0) 42
(33,29)
5.77633
×10−4
(-0.0625,0.1787514) 136
(101,91)
5.83623
×10−5
(-0.04,0.1666...) 405
(321,281)
4.73001
×10−6
(-0.09375,0.1214286) 1191
(1001,901)
2.33608
×10−4
(0,222...) Non con e xencia
Cad o 4.3: Resul ados numé icos pa a o es 3 conside ando
ω= 2
,
α(µ, z) = |sen(12µz)|
e
σ(µ, z) = 1 + sen(12µz) cos(12µz)
.
.
Figu a 4.2: Solución ap oximada ob ida e e o come ido pa a o es 3 no caso
(I, N) =
(11,10)
.
- TEST 4.
P oblema de Kim-T anquilli:
Na e e encia [1] Kim e T anquilli desc iben emp egando a ecuación de Fokke -Planck
o mo emen o de p opagación da luz no in e io dun ecido biolóxico. P esen amos na
gu a 4.3 imaxes que nos pe mi en compa a os seus esul ados.
Pa a es e p oblema conside amos de no o
ω= 2
e as seguin es uncións:
α(µ, z) =
0,02
,
σ(µ, z)=0,01
,
(µ)=1
,
g(µ)=2
e
W(µ, z)=0.

26
CAPÍTULO 4. RESULTADOS NUMÉRICOS
Figu a 4.3: Solución ap oximada ob ida e solución ap oximada en
z=Z
n
= 1
pa a o es
de Kim-T anquilli no caso
(I, N) = (21,20)
.
- TEST 5.
Ago a amos comp oba a ecacia do algo i mo cambiando o pa áme o de elaxación
ω
pa a obse a como a ec a di o pa áme o ao núme o de i e acións necesa ias pa a
acada a con e xencia.
Conside amos de no o a unción do TEST 3:
ψ(µ, z) = ln(2 + µ2+z3).
Conside amos en ón amén as seguin es uncións:
(µ) = ln(2 + µ2),
g(µ) = ln(3 + µ2),
W(µ, z) = 3µz2
2 + µ2+z3+|sen(12µz)|ln(2 + µ2+z3)−
−(1 + sen(12µz) cos(12µz)) −4µ2
2 + µ2+z3−2(1 −µ2)µ2−z3−2
(2 + µ2+z3)2.
(4.1)
No cad o 4.4 apa ece epe idas eces o e mo ines able. Es e e mo o que nos indica é
que o esquema numé ico explo a no sen ido de que se ol e ines able. Na gu as 4.4
27
I e acións con
(I, N) = (11,10) E(11,10)
abs
(Q)
I e acións con
(I, N) = (33,29) E(33,29)
abs
(Q)
ω= 1
87 5.05206
×10−3
264 5.77510
×10−4
ω= 1,5
57 5.05207
×10−3
180 5.77600
×10−4
ω= 2
42 5.05208
×10−3
136 5.77633
×10−4
ω= 2,5
113 5.05210
×10−3
109 5.77652
×10−4
ω= 3
Di e xencia 91 5.77668
×10−4
ω= 3,5
Di e xencia 78 5.77681
×10−4
ω= 4
Di e xencia Di e xencia
Cad o 4.4: Resul ados numé icos pa a o es 5.
e 4.5 ep esen amos un exemplo des a ines abilidade ocasionada ao i aumen ando o
pa áme o de elaxación
ω
.
Ademais des a eexión e des e úl imo es amén é in e esan e comp oba o que
oco e cos esul ados numé icos cando modicamos a semen e. Así pois as múl iples
expe imen os numé icos a elección ei a nes e a igo é a que mello es esul ados de
con e xencia p opo ciona. Conside a emos logo un no o es no que se mos a unha
elección de semen e dis in a.
Figu a 4.4: Solución ap oximada ob ida e solución ap oximada en
z=Z
n
= 1
pa a o es
5 no caso
(I, N) = (11,10)
e
ω= 3
.
28
CAPÍTULO 4. RESULTADOS NUMÉRICOS
Figu a 4.5: Solución exac a e e o come ido ap oximando pa a o es 5 no caso
(I, N) =
(11,10)
e
ω= 3
.
- TEST 6.
Nes e úl imo es amos mos a o que oco e cando escollemos unha semen e dis in a
á p opo cionada pola ecuación (3.8).
Conside emos en ón a seguin e semen e:
(ψn
i∗)[0] =a c an(π)
1 + π2
pa a
n∈ {2, ... , N −1},
(4.2)
onde
i∗=I+1
2
e de no o omamos a unción
ψ
do TEST 3:
ψ(µ, z) = ln(2 + µ2+z3).
De no o considé anse as uncións
,
g
e
W
dadas polas ecuacións (4.1).
(I, N)E
abs
(Q) (µ, z)
I e acións
(11,10)
7.05208
×10−3
(-1,0) 42
(33,29)
5.77633
×10−4
(-0.0625,0.1787514) 155
(101,91)
5.83623
×10−5
(-0.04,0.1666...) 520
(321,281)
4.73001
×10−6
(-0.09375,0.1214286) 1697
Cad o 4.5: Resul ados numé icos pa a o es 6 conside ando
ω= 2
e a no a semen e.
29
Como podemos e no cad o 4.5 ao cambia a semen e e escolle a dada pola ecuación
(4.2) obse amos como o núme o de i e acións necesa ias pa a acada a con e xencia
aumen ou. Con is os esul ados mos amos así que a semen e dada pola ecuación
(3.8) p opo ciona mello es esul ados.
Sexan nalmen e as seguin es semen es:
(ψn
i∗)[0] =π
1 + a c an(eπ)
pa a
n∈ {2, ... , N −1},
(4.3)
(ψn
i∗)[0] = 106
pa a
n∈ {2, ... , N −1}
(4.4)
e
(ψn
i∗)[0] =−5×106
pa a
n∈ {2, ... , N −1}.
(4.5)
Como podemos obse a no cad o 4.6 o algo i mo segue le ándonos á solución a pe-
sa es de conside a unha semen e moi lonxe da óp ima como poden se as dadas pola
ecuación (4.4) ou pola (4.5). Así pois obse amos o compo amen o de aumen o de
i e acións pa a acada a solución a medida que conside amos semen es máis dispa es.
Acabamos de e así que non é necesa io oma unha semen e p óxima aos alo es da
solución exac a no segmen o
Q0
pa a que o algo i mo con e xa, ei o que é compa ible
coa con e xencia global.
(I, N)
I e acións (3.8) I e acións (4.3) I e acións (4.4) I e acións (4.5)
(11,10)
42 51 88 92
(33,29)
136 184 316 331
(101,91)
405 600 1030 1081
(321,281)
1191 1945 3339 3503
Cad o 4.6: Resul ados numé icos pa a o es 6 conside ando
ω= 2
e as dis in as semen es.
Ademais pa a es e es conside amos
i max = 4000
.
36
ANEXO I: FÓRMULAS DE DERIVACIÓN NUMÉRICA
Po én as ó mulas que se emp egan pa a disc e iza o e mo
h(1 −µ2)∂ψ
∂µ i
cando nos
a opamos nos casos
i= 1
e
i=I
son, espec i amen e, a ó mula (3) e a (4). Di as
ó mulas ap oximan a de i ada de p imei a o de emp egando nodos ca a adian e e ca a
a ás, e nalmen e emp egamos a ó mula (3) pa a disc e iza a de i ada de p imei a o de
que nos apa ece.
Se nos xamos na de i ación do esquema pa a es es dous úl imos nodos non emp ega-
mos as ó mulas comple as, pois an o
D1
coma
DI
son nulos e polo an o o e mo que
acompaña a
(α)
amén o é.

Anexo II: Códigos Ma lab
No anexo p esen e a amos de amosa liñas de código Ma lab que nos pe mi i on
ob e esul ados numé icos pa a o noso p oblema emp egando o algo i mo i e a i o que
deseñamos con an e io idade.
C eamos un p og ama p incipal, o cal chama a dis in as uncións que pe mi en esol e
as dis in as e apas do mé odo. A con inuación mós anse dis in as imaxes que eexan di a
es u u a:
Figu a 6: Nes e p og ama p incipal chamamos á unc ion que ealiza os cálculos e as dis-
in as unc ions que emp egamos pa a a ep esen ación dos esul ados ob idos.
37
38
ANEXO II: CÓDIGOS MATLAB
Figu a 7: Pa e p incipal da unc ion que ealiza o mé odo i e a i o.
Como podemos e na gu a 7 o bucle do mé odo i e a i o chama a dis in as uncións.
Cada unha des as uncións se e pa a ealiza as e apas do mé odo i e a i o desc i as no
Capí ulo 3.
A unc ion Qplus esol e o p oblema no conxun o
Q+
, di o dou o xei o esol e o
p oblema de alo inicial dado pola unción
.
Po ou a banda, a unc ion Qminus esol e o p oblema en
Q−
, onde emos un p oblema
de alo nal dado pola unción
g
.
Finalmen e, chamamos á unc ion ac ualiza que se e pa a ac ualiza os alo es da
solución emp egando o mé odo de elaxación desc i o.
Na gu a 8 obse amos o bucle en
n
que esol e un sis ema po ascenso pa a ob e
as ap oximacións da solución en
Q+
. No bucle denimos unha ma iz, na cal a mei ande
pa e das en adas son nulas, e un ec o segundo memb o, a pa i es da disc e ización da
ecuación. Finalmen e, pa a ob e as ap oximacións esol emos a ecuación ma icial dada
pola ma iz e o ec o an e io men e mencionados.
39
Figu a 8: Es a o da unc ion Qplus no cal se amosa a cons ución do sis ema linea a
esol e .
Analogamen e, na unc ion Qminus denimos unha ma iz e un ec o segundo memb o
pa a esol e o sis ema e ob e a ap oximación da solución en
Q−
, pe o nes e caso o sis ema
que se esol e é po descenso.
Na gu a 9 obse amos es a pequena di e encia, onde o sis ema é esol o po descenso.
40
ANEXO II: CÓDIGOS MATLAB
Figu a 9: Es a o da unc ion Qminus no cal se amosa o bucle que esol e o PVF denido
pola unción
g
.
Despois des es cálculos chamamos á unc ion da gu a 10 que ac ualiza os alo es da
solución pa a
i∗=I+1
2
.
Figu a 10: Es a unción co esponde coa ac ualización emp egando o pa áme o de elaxa-
ción
ω
.
T as odo es e p oceso ob emos unha g an ma iz
U
que almacena a solución ap oxi-
mada pa a o mallado c eado. Ap obei a emos di a ma iz pa a ob e e os come idos na
ap oximación, pa a así nalmen e ep esen a di os e os e a solución ob ida.
41
Ademais das unc ions de cálculo comezamos con unc ions pa a in oduci os da os do
p oblema como poden se os e mo
D(µ)
ou as uncións
α
,
σ
e
W
.
Figu a 11: Exemplos de on es/sumidoi os que emp egamos nos nosos es .
Figu a 12: Amosamos un exemplo do e mo de di usi idade ex e na emp egada nos nosos
expe imen os numé icos.

42
ANEXO II: CÓDIGOS MATLAB
Bibliog a ía
[1]
A. D. Kim
e
P. T anquilli
.
Nume ical solu ion o he Fokke -Planck equa ion wi h
a iable coecien s, Jou nal o Quan i a i e Spec oscopy & Radia i e T ans e
109
,
no. 5 (2008) 727-740.
[2]
Ó. López Pouso
e
N. Jumaniyazo
.
Nume ical expe imen s wi h he Fokke -Planck
equa ion in 1D slab geome y, Jou nal o Compu a ional and Theo e ical T anspo
45
,
no. 3 (2016) 184-201.
[3]
Ó. López Pouso
e
N. Jumaniyazo
.
Di ec e sus i e a i e me hods o o wa d-
backwa d equa ions. Nume ical compa isons on a pa icula anspo kine ic model
,
en p oceso de publicación.
[4]
D. S ein
e
I. B. Be ns ein
.
Bounda y alue p oblem in ol ing a simple Fokke -
Planck equa ion, Physics o Fluids
19
(1976) 811-814.
[5]
V. Vanaja
e
R. B. Kellogg
.
I e a i e me hods o a o wa d-backwa d hea equa ion,
SIAM Jou nal on Nume ical Analysis
27
, no. 3 (1990) 622-635.
43