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

Traballo Fin de Grao Resolución numérica de ecuacións de Fokker-Planck nunha dimensión espacial Javier López Pedrares 2019/2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Resolución numérica de ecuacións de Fokker-Planck nunha dimensión espacial Javier López Pedrares Xullo 2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Traballo proposto Área de Coñecemento: Matemática Aplicada Título: Resolución numérica de ecuacións de Fokker-Planck nunha dimensión espacial Breve descrición do contido Resolveranse algúns modelos de transporte de partículas cargadas, como electróns ou ións pesados. Os métodos numéricos estarán baseados en métodos de diferenzas nitas. Titor: Óscar López Pouso iii Índice xeral Resumo vii Introdución ix 1. Esquema numérico 1 2. Interpretación física 9 3. Implementación do esquema 13 4. Resultados numéricos 21 Anexos 31 Anexo I: Fórmulas de derivación numérica 33 Anexo II: Códigos Matlab 37 Bibliografía 43 v Resumo O traballo comeza introducindo un coñecido problema no mundo da física, o cal formulamos matematicamente. Posteriormente continuamos describindo un esquema numérico que permite obter unha solución para o problema baseado nas diferenzas nitas. Ademais dunha formulación matemática apórtase unha visión e unha explicación física do problema. Finalmente, describimos o algoritmo iterativo que resolve o problema e implementámolo en Matlab para obter resultados numéricos. Abstract The article starts introducing a well-known physics problem, that we go to formulate mathematically. Then, we continue describing a numeric scheme that allows obtain a solution of our problem based on nite dierences. In addition to a mathematical formulation we show a little sight and a physical explanation of the problem. At the end, we nish with the iterative algorithm's description that solves the problem and we implement it on Matlab to obtain numerical results. vii 4 CAPÍTULO 1. ESQUEMA NUMÉRICO Así acabamos de discretizar unha primeira parte do noso problema, agora precisamos discretizar o termo da difusividade. Dependendo do punto da malla onde nos atopemos este termo pode ser nulo (cando µ=±1 ); entón consideramos os seguintes casos: (1) Se nos atopamos en i= 1 , ou sexa µ1=−1 temos que D(µ1) = 1−µ2 1= 1−(−1)2= 1−1=0 , i.e., D1= 0 . Neste caso podemos realizar a seguinte discretización: ∂ ∂µ D(µ)∂ψ ∂µ (µ,z)=(µ1,zm) O(h2) ≈4D2∂ψ ∂µ (µ2, zm)−D3∂ψ ∂µ (µ3, zm) 2h (1.7) e cando nos atopamos con r∈ {2,3} , empregamos a fórmula estándar centrada en dous puntos para discretizar as derivadas de primeira orde que aparecen na ecuación anterior: ∂ψ ∂µ (µr, zm) O(h2) ≈ψm r+1 −ψm r−1 2h. (1.8) (2) Aquí analizamos os nodos interiores da malla, ou sexa, i∈ {2, ... , I −1} . Aquí basta empregar a fórmula estándar de segunda orde centrada: ∂ ∂µ 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 estudar a discretización no nodo nal, ou sexa en µI= 1 . Atopamos de novo unha difusión nula, logo procedemos analogamente que no primeiro caso: ∂ ∂µ D(µ)∂ψ ∂µ (µ,z)=(µI,zm) O(h2) ≈DI−2∂ψ ∂µ (µI−2, zm)−4DI−1∂ψ ∂µ (µI−1, zm) 2h (1.10) e de novo para r∈ {I−2, I −1} procédese igual que en (1.8). Agora que xa temos discretizadas as derivadas do noso problema, podemos proceder a describir o esquema numérico completo que nos leva a atopar a solución ψ . Teremos que prestar especial atención ao número de ecuacións resultantes, pois dependendo da paridade de I atoparemos que o número de ecuacións non coincide co de incógnitas. En vista ás formulas de discretización sinaladas antes precisamos esixir alomenos que I≥4 e N > 2 . Procedamos así a describir o esquema numérico de novo separando en casos como xemos para discretizar: 5 - Para (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) - Para (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) - Para (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í completamente descrito o esquema numérico que nos permite resolver o problema denido polas ecuacións (1), (2) e (3). Pero atopámonos cun problema pois nun determinado caso o número de ecuacións non coincide co de incógnitas. 6 CAPÍTULO 1. ESQUEMA NUMÉRICO Para estudar dita problemática basta considerar un I impar. O número de incógnitas que temos é I×N , pero pola contra o esquema que acabamos de describir unicamente posúe I×(N−1) ecuacións. Agora ben temos que comprobar se as ecuacións (condicións de uxo entrante) (2) e (3) proporcionan ecuacións sucientes para solventar este desaxuste. Así pois dependerá da paridade de I o número de ecuacións proporcionadas. Cando I é par porporcionan I ecuacións, pero cando é I impar unicamente proporcionan I−1 ecuacións. Así temos que cando I é impar falta unha ecuación e non podemos resolver o esquema proposto anteriormente. Se I é par resolvemos o problema empregando o esquema tal como foi descrito. A continuación describimos como completar o esquema para I impar. Sexa I impar, logo é obvio que µI+1 2= 0 simplemente por construción. En primeira instancia podemos pensar en empregar unha das condicións numéricas inicial ou nal: ψ1 I+1 2 =f(0) ou ψN I+1 2 =g(0), (1.14) pero en realidade a interese radica en impoñer as dúas condicións debido a que estamos interesados en obter solucións continuas. Entón o que facemos é eliminar as N−1 ecuacións correspondentes a i=I+1 2 para impoñer as condicións numéricas anteriores e así considerar un novo conxunto de N−2 ecuacións para µ= 0 e os nodos do intervalo (Z ini , Z n ) . Entón para µ= 0 a ecuación (1) reescríbese como: α(0, z)ψ(0, z)−σ(0, z)∂2ψ ∂µ2(0, z) = W(0, z) para z∈[Z ini , Z n ]. (1.15) Vexamos que efectivamente se verica o anterior. Para substituír en µ= 0 precisamos antes derivar o termo que leva a difusión empregando a regra de derivación do produto e posteriormente evaluar en µ= 0 . Derivemos entón: ∂ ∂µ (1 −µ2)∂ψ ∂µ =−2µ∂ψ ∂µ + (1 −µ2)∂ ∂µ ∂ψ ∂µ =−2µ∂ψ ∂µ + (1 −µ2)∂2ψ ∂µ2. (1.16) Agora se substituímos para µ= 0 temos: ∂ ∂µ (1 −µ2)∂ψ ∂µ (µ,z)=(0,z) =−2·0∂ψ ∂µ (0, z) + (1 −02)∂2ψ ∂µ2(0, z) = ∂2ψ ∂µ2(0, z), (1.17) que efectivamente é o que empregamos en (1.15). Se pensamos no intervalo aberto (Z ini , Z n ) a ecuación (1.15) suxire: αn i∗ψn i∗−σn i∗ ψn i∗−1−2ψn i∗+ψn i∗+1 h2=Wn i∗ para n∈ {2, ... , N −1}, (1.18) 7 sendo i∗=I+1 2 . O anterior pode reescribirse do seguinte xeito: −σn i∗ h2ψn i∗−1+αn i∗+2σn i∗ h2ψn i∗+−σn i∗ h2ψn i∗+1 =Wn i∗ (1.19) para n∈ {2, ... , N −1} . Entón o esquema nal para o caso I impar consiste no descrito restándolle as N−1 ecuacións correspondentes a i=I+1 2 , engandindo as condicións de contorno numéricas (1.14) e as N−2 ecuacións (1.19) que acabamos de obter. Como comentaremos máis adiante implementaremos o esquema para o caso I impar pois obtéñense mellores resultados. 8 CAPÍTULO 1. ESQUEMA NUMÉRICO Capítulo 2 Interpretación física A ecuación de Fokker-Planck (EFP) ten unha iteresante interpretación física. Imos mostrar a relación do noso problema coa ecuación do transporte de Bolzmann (ETB). A importancia do problema de Fokker-Planck radica pois en que a solución deste, ψ , é unha aproximación da solución do problema de Bolzmann en certas condicións. O interesante é estudar ditas condicións. Así pois a EFP é unha aproximación da de Bolzmann cando falamos do transporte de partículas que sofren pequenas desviacións na súa traxectoria ao impactar con outras partículas e pequenas perdas de carga. Neste caso falamos de dispersión ou scattering. Exemplos comúns de partículas que sofren ditos comportamentos son partículas cargadas como ións pesados ou electróns da codia do átomo que posúen carga negativa. Ámbalas dúas ecuacións anteriores describen a densidade de uxo angular de partículas, ψ . Ademais establecen o equilibrio entre as perdas e as ganancias de ψ ao longo das direccións de propagación. Dito doutro xeito estudan o gradiente de ψ ao longo de cada dirección do espazo tridimensional R3 . Consideremos Ω⊂R3 un dominio espacial e tomemos ω como dirección de propagación. Se escribimos ω en coordeadas esféricas temos: ω=ω(ϕ, θ) = (sen ϕcos θ, sen ϕsen θ, cos ϕ)∈ S 2, sendo S 2 a esfera unitaria do espazo con ϕ∈[0, π] o ángulo polar e θ∈[0,2π) o ángulo acimutal. Se o uxo depende da posición x= (x1, x2, x3)∈Ω unicamente a través de x3=z , temos que ω· ∇ψ=ω3∂ψ ∂z . Se na ecuación (1) consideramos µ=ω3 obtemos exactamente o anterior, sendo ω a dirección de propagación. No noso caso estamos traballando sobre a ecuación unidimensional, entón o noso dominio será da forma Ω = R2×(Z ini , Z n ) , o 9 10 CAPÍTULO 2. INTERPRETACIÓN FÍSICA Figura 2.1: Sección do slab unidimensional. cal adoita chamarse slab unidimensional. Unha pequena ilustración do anterior é a gura 2.1. Agora ben, partimos de tres variables espaciais e imos reducilas a unicamente unha. Isto pódese facer grazas a que o uxo, ψ , non varía nin con x1 nin con x2 e ademais o dominio anterior son copias do dominio unidimensional (Z ini , Z n ) . Se continuamos co descrito temos ω3=µ= cos ϕ , con isto observamos trivialmente a razón de porque o parámetro µ varía entre 1 e −1 pois µ é o coseno dun ángulo. Ao pensar na EFP unidimensional deixamos de depender do ángulo acimutal, entón unicamente temos o ángulo polar ou equivalentemente µ e así o operador ∂ ∂µ (1 −µ2)∂· ∂µ (2.1) coñecido comunmente como o operador continuo de dispersión, o cal non é máis que o laplaciano sobre a esfera. Podemos dicir así que a ecuación (1) é a EFP monoenerxética supoñendo que a enerxía é constante no dominio unidimensional con simetría de xeometría plana e as ecuacións (2) e (3) impoñen a densidade de uxo de partículas que entra no dominio. Por isto último ditas ecuacións teñen o nome de condicións de uxo entrante. Como mencionamos ao nal da Introdución o noso problema ten unha característica peculiar dentro das ecuacións en derivadas parciais: a ausencia de condicións de contorno para |µ|= 1 , ou sexa o operador (2.1) non proporciona ditas condicións. Isto débese a que o equilibrio dado pola EFP é suciente para obter o uxo angular unha vez coñecemos o uxo a través da fronteira física. Entón impoñendo unhas condicións de contorno o que 11 faríamos sería sobredeterminar o problema. Finalmente poñamos un signicado físico a cada un dos termos que aparecen na ecuación (1): - Por ψ(µ, z) denotamos á densidade de uxo angular de partículas en (µ, z) . Dito doutro xeito é o número de partículas que se moven desde z ao longo da dirección ω(µ) por unidade de área normal a ω(µ) , por unidade de ángulo sólido 1 , por unidade de enerxía e por unidade de tempo. - Cando falamos de movernos dende z na dirección de ω(µ) referímonos claramente a movernos en (x1, x2, z) para calquera (x1, x2)∈R2 xos e ao longo da dirección ω(µ) = p1−µ2cos θ, p1−µ2sen θ, µ para calquera ángulo θ∈[0,2π) xo. - O primeiro termo da ecuación, µ∂ψ ∂z , é a derivada direccional de ψ(µ, z) ao longo da dirección ω(µ) . - Denotamos por α ao coeciente de absorción, así pois αψ describe as perdas por absorción. - O termo σ∂ ∂µ h(1 −µ2)∂ψ ∂µ i é o termo de difusión na variable angular. - A función de dúas variables W representa unha fonte interna de partículas onde teña sigmo positivo, e un sumidoiro interno onde o signo sexa negativo. 1 Por ángulo sólido entendemos un ángulo espacial que abarca un obxeto cando é visto dende un punto. 12 CAPÍTULO 2. INTERPRETACIÓN FÍSICA Capítulo 3 Implementación do esquema Como mencionamos no Capítulo 1 obtivemos dous esquemas distintos segundo a paridade de I . Tras a realización de diversos experimentos numéricos chegouse á conclusión de que era máis recomendable a implantación do esquema para I impar. Tras os diversos experimentos realizados observouse que para o caso I impar non aparecían inestabilidades arredor de µ= 0 . Porén, no caso I par cando nos achegamos a Q0 poden aparecer solucións que son inestables, presentando así oscilacións espurias. En termos de implementación do esquema poderíamos considerar dous tipos de esquema: - Método directo: un esquema no que resolvemos a ecuación en todo o dominio, Q . Consideramos unha gran matriz cadrada de dimensión I×N e tentamos en resolver o sistema linear completo asociado ao problema. - Método iterativo: algoritmo no que se resolve o problema partindo dunha semente inicial en Q0 , resolvemos por ascenso e descenso, respectivamente en Q+ e Q− , os problemas de valor inicial e nal denidos polas funcións f e g . Realizamos este proceso ata que os valores obtidos en dúas iteracións sucesivas para a solución en Q0 disten unha cantidade real positiva épsilon moi pequena. Ou tamén se non se acada converxencia dado un número de iteracións máximas. Neste traballo centrarémonos na implementación do algoritmo iterativo para o caso I impar. Consideremos entón de novo os conxuntos do dominio Q descritos ao comezo do Capítulo 1. Sexan entón Q+= (0,1] ×[Z ini , Z n ] , o segmento Q0={0} × [Z ini , Z n ] e Q−= [−1,0) ×[Z ini , Z n ] (ver gura 3.1). 13 20 CAPÍTULO 3. IMPLEMENTACIÓN DO ESQUEMA Capítulo 4 Resultados numéricos Nesta sección mostraremos unha serie de resultados numéricos obtidos a partires do algoritmo iterativo que implementamos en Matlab. Para todos os experimentos consideraremos Z ini = 0 , Z n = 1 , itmax = 2000 (número de iteracións máximas) e tomaramos a seguinte tolerancia, ε= 10−8 . Ademais, agás que así se indique, en todos os tests realizados emprégase como semente para inicializar o algoritmo a dada pola ecuación (3.8). En canto notación, por E abs (Q) entendemos o seguinte: E abs (Q) = m´ax Q|ψmalla −ψ|, sendo ψmalla a solución aproximada e o máximo tómase sobre o conxunto de todos os nodos da malla. As fórmulas que foron empregadas para a obtención do esquema numérico son exactas cando a solución exacta é polinómica de grao menor ou igual que 1 na variable µ e de grao menor ou igual que 2 para a variable z . - TEST 1. Supoñamos que a solución exacta é coñecida, entón sexa: ψ(µ, z) = µ2z2. Logo realizando uns simples cálculos obtemos as funcións f , g e W substituíndo nas 21 22 CAPÍTULO 4. RESULTADOS NUMÉRICOS ecuacións (1)-(3): f(µ) = ψ(µ, Z ini ) = ψ(µ, 0) = 0, g(µ) = ψ(µ, Z n ) = ψ(µ, 1) = µ2, W(µ, z) = µ∂ψ ∂z +αψ −σ∂ ∂µ (1 −µ2)∂ψ ∂µ = = 2µz3+α(µ, z)µ2z2−σ(µ, z)z2(2 −6µ2), considerando: α(µ, z) = |sen(12µz)| e σ(µ, z) = 1 + sen(12µz) cos(12µz). Como podemos observar no cadro 4.1 aparece a terminoloxía (µ, z) , con isto referímonos ao punto da malla (µi, zn) no que se atopa o erro máximo en valor absoluto. Tomaremos dita notación para os posteriores test e cadros. Ademais tamén se intenta reexar a orde de converxencia do esquema para iso realízase o seguinte: Sexa E(h, k) o erro máximo cometido na aproximación para unha malla de tamaño h - k e sexa c unha constante, logo se a orde é p tanto para h coma para k debe vericarse asintoticamente ( h, k →0 ) a seguinte igualdade: E(ch, ck) = cpE(h, k). Tomando logaritmos podemos despexar a orde, tense así: p= ln E(ch,ck) E(h,k) ln c. (I, N)E abs (Q) (µ, z) Iteracións Orde (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 Cadro 4.1: Resultados numéricos para o test 1 considerando ω= 2 . 23 Por exemplo se consideramos h 2 e k 2 , temos a seguinte constante c=1 2 . Para este exemplo observamos que ψ é polinómica de grao 2 na variable z . Pero non é polinómica de grao menor ou igual que 1 na variable µ , entón para obter mellores resultados o que facemos é renar a malla en dita variable, ou sexa en h . Dito comportamento reíctese tamén no cadro 4.1. - TEST 2. Supoñamos de novo que coñecemos a solución exacta ao problema, sexa así: ψ(µ, z) = µz3. Realizando os cálculos pertinentes obtemos: f(µ) = ψ(µ, Z ini ) = ψ(µ, 0) = 0, g(µ) = ψ(µ, Z n ) = ψ(µ, 1) = µ, W(µ, z) = µ∂ψ ∂z +αψ −σ∂ ∂µ (1 −µ2)∂ψ ∂µ = = 3µ2z2+µz3(α(µ, z)−2σ(µ, z)) . Como podemos observar nos cadros 4.2 e 4.3 aparece o termo non converxencia, con isto referímonos a que non se acada a converxencia do método chegado o número de iteracións máximas. Porén, o algoritmo está a comportarse de maneira converxente pois o residuo calculado tras a iteración nal achégase á tolerancia. Por exemplo para o cadro 4.2 o último residuo calculado, onde aparece a non converxencia, é res = 1.18673 ×10−7 , que como observamos achégase á tolerancia épsilon. (I, N)E abs (Q) (µ, z) Iteracións Orde (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 converxencia 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 Cadro 4.2: Resultados numéricos para o test 2 considerando ω= 2 e tomando de novo α(µ, z) = |sen(12µz)| e σ(µ, z) = 1 + sen(12µz) cos(12µz) . 24 CAPÍTULO 4. RESULTADOS NUMÉRICOS Figura 4.1: Solución aproximada obtida e erro cometido para o test 2 no caso (I, N) = (1001,901) . Neste segundo test observamos que a solución exacta ψ é polinómica de grao 1 na variable µ , porén non é de grao menor ou igual que 2 na variable z , logo aquí é de interese renar a malla en z , ou sexa en k , como facemos no cadro 4.2. - TEST 3. Supoñamos que a solución exacta é coñecida, entón sexa: ψ(µ, z) = ln(2 + µ2+z3). Logo realizando uns sinxelos cálculos obtemos: f(µ) = ψ(µ, 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 novo no cadro 4.3 aparece a non converxencia debido a que non se satisfai o test de parada chegado o número máximo de iteracións. Neste caso o último residuo calculado é res = 7.17688 ×10−7 25 (I, N)E abs (Q) (µ, z) Iteració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 converxencia Cadro 4.3: Resultados numéricos para o test 3 considerando ω= 2 , α(µ, z) = |sen(12µz)| e σ(µ, z) = 1 + sen(12µz) cos(12µz) . . Figura 4.2: Solución aproximada obtida e erro cometido para o test 3 no caso (I, N) = (11,10) . - TEST 4. Problema de Kim-Tranquilli: Na referencia [1] Kim e Tranquilli describen empregando a ecuación de Fokker-Planck o movemento de propagación da luz no interior dun tecido biolóxico. Presentamos na gura 4.3 imaxes que nos permiten comparar os seus resultados. Para este problema consideramos de novo ω= 2 e as seguintes funcións: α(µ, z) = 0,02 , σ(µ, z)=0,01 , f(µ)=1 , g(µ)=2 e W(µ, z)=0. 26 CAPÍTULO 4. RESULTADOS NUMÉRICOS Figura 4.3: Solución aproximada obtida e solución aproximada en z=Z n = 1 para o test de Kim-Tranquilli no caso (I, N) = (21,20) . - TEST 5. Agora vamos comprobar a ecacia do algoritmo cambiando o parámetro de relaxación ω para observar como afecta dito parámetro ao número de iteracións necesarias para acadar a converxencia. Consideramos de novo a función do TEST 3: ψ(µ, z) = ln(2 + µ2+z3). Consideramos entón tamén as seguintes funcións: f(µ) = 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 cadro 4.4 aparece repetidas veces o termo inestable. Este termo o que nos indica é que o esquema numérico explota no sentido de que se volve inestable. Na guras 4.4 27 Iteracións con (I, N) = (11,10) E(11,10) abs (Q) Iteració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 Diverxencia 91 5.77668 ×10−4 ω= 3,5 Diverxencia 78 5.77681 ×10−4 ω= 4 Diverxencia Diverxencia Cadro 4.4: Resultados numéricos para o test 5. e 4.5 representamos un exemplo desta inestabilidade ocasionada ao ir aumentando o parámetro de relaxación ω . Ademais desta reexión e deste último test tamén é interesante comprobar o que ocorre cos resultados numéricos cando modicamos a semente. Así pois tras múltiples experimentos numéricos a elección feita neste artigo é a que mellores resultados de converxencia proporciona. Consideraremos logo un novo test no que se mostra unha elección de semente distinta. Figura 4.4: Solución aproximada obtida e solución aproximada en z=Z n = 1 para o test 5 no caso (I, N) = (11,10) e ω= 3 . 28 CAPÍTULO 4. RESULTADOS NUMÉRICOS Figura 4.5: Solución exacta e erro cometido aproximando para o test 5 no caso (I, N) = (11,10) e ω= 3 . - TEST 6. Neste último test vamos mostrar o que ocorre cando escollemos unha semente distinta á proporcionada pola ecuación (3.8). Consideremos entón a seguinte semente: (ψn i∗)[0] =arctan(π) 1 + π2 para n∈ {2, ... , N −1}, (4.2) onde i∗=I+1 2 e de novo tomamos a función ψ do TEST 3: ψ(µ, z) = ln(2 + µ2+z3). De novo considéranse as funcións f , g e W dadas polas ecuacións (4.1). (I, N)E abs (Q) (µ, z) Iteració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 Cadro 4.5: Resultados numéricos para o test 6 considerando ω= 2 e a nova semente. 29 Como podemos ver no cadro 4.5 ao cambiar a semente e escoller a dada pola ecuación (4.2) observamos como o número de iteracións necesarias para acadar a converxencia aumentou. Con istos resultados mostramos así que a semente dada pola ecuación (3.8) proporciona mellores resultados. Sexan nalmente as seguintes sementes: (ψn i∗)[0] =π 1 + arctan(eπ) para n∈ {2, ... , N −1}, (4.3) (ψn i∗)[0] = 106 para n∈ {2, ... , N −1} (4.4) e (ψn i∗)[0] =−5×106 para n∈ {2, ... , N −1}. (4.5) Como podemos observar no cadro 4.6 o algoritmo segue levándonos á solución a pesares de considerar unha semente moi lonxe da óptima como poden ser as dadas pola ecuación (4.4) ou pola (4.5). Así pois observamos o comportamento de aumento de iteracións para acadar a solución a medida que consideramos sementes máis dispares. Acabamos de ver así que non é necesario tomar unha semente próxima aos valores da solución exacta no segmento Q0 para que o algoritmo converxa, feito que é compatible coa converxencia global. (I, N) Iteracións (3.8) Iteracións (4.3) Iteracións (4.4) Iteració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 Cadro 4.6: Resultados numéricos para o test 6 considerando ω= 2 e as distintas sementes. Ademais para este test consideramos itmax = 4000 . 36 ANEXO I: FÓRMULAS DE DERIVACIÓN NUMÉRICA Porén as fórmulas que se empregan para discretizar o termo h(1 −µ2)∂ψ ∂µ i cando nos atopamos nos casos i= 1 e i=I son, respectivamente, a fórmula (3) e a (4). Ditas fórmulas aproximan a derivada de primeira orde empregando nodos cara adiante e cara atrás, e nalmente empregamos a fórmula (3) para discretizar a derivada de primeira orde que nos aparece. Se nos xamos na derivación do esquema para estes dous últimos nodos non empregamos as fórmulas completas, pois tanto D1 coma DI son nulos e polo tanto o termo que acompaña a f(α) tamén o é. Anexo II: Códigos Matlab No anexo presente tratamos de amosar liñas de código Matlab que nos permitiron obter resultados numéricos para o noso problema empregando o algoritmo iterativo que deseñamos con anterioridade. Creamos un programa principal, o cal chama a distintas funcións que permiten resolver as distintas etapas do método. A continuación móstranse distintas imaxes que reexan dita estrutura: Figura 6: Neste programa principal chamamos á function que realiza os cálculos e as distintas functions que empregamos para a representación dos resultados obtidos. 37 38 ANEXO II: CÓDIGOS MATLAB Figura 7: Parte principal da function que realiza o método iterativo. Como podemos ver na gura 7 o bucle do método iterativo chama a distintas funcións. Cada unha destas funcións serve para realizar as etapas do método iterativo descritas no Capítulo 3. A function Qplus resolve o problema no conxunto Q+ , dito doutro xeito resolve o problema de valor inicial dado pola función f . Por outra banda, a function Qminus resolve o problema en Q− , onde temos un problema de valor nal dado pola función g . Finalmente, chamamos á function actualizar que serve para actualizar os valores da solución empregando o método de relaxación descrito. Na gura 8 observamos o bucle en n que resolve un sistema por ascenso para obter as aproximacións da solución en Q+ . No bucle denimos unha matriz, na cal a meirande parte das entradas son nulas, e un vector segundo membro, a partires da discretización da ecuación. Finalmente, para obter as aproximacións resolvemos a ecuación matricial dada pola matriz e o vector anteriormente mencionados. 39 Figura 8: Estrato da function Qplus no cal se amosa a construción do sistema linear a resolver. Analogamente, na function Qminus denimos unha matriz e un vector segundo membro para resolver o sistema e obter a aproximación da solución en Q− , pero neste caso o sistema que se resolve é por descenso. Na gura 9 observamos esta pequena diferencia, onde o sistema é resolto por descenso. 40 ANEXO II: CÓDIGOS MATLAB Figura 9: Estrato da function Qminus no cal se amosa o bucle que resolve o PVF denido pola función g . Despois destes cálculos chamamos á function da gura 10 que actualiza os valores da solución para i∗=I+1 2 . Figura 10: Esta función corresponde coa actualización empregando o parámetro de relaxación ω . Tras todo este proceso obtemos unha gran matriz U que almacena a solución aproximada para o mallado creado. Aprobeitaremos dita matriz para obter erros cometidos na aproximación, para así nalmente representar ditos erros e a solución obtida. 41 Ademais das functions de cálculo comezamos con functions para introducir os datos do problema como poden ser os termo D(µ) ou as funcións α , σ e W . Figura 11: Exemplos de fontes/sumidoiros que empregamos nos nosos test. Figura 12: Amosamos un exemplo do termo de difusividade externa empregada nos nosos experimentos numéricos. 42 ANEXO II: CÓDIGOS MATLAB Bibliografía [1] A. D. Kim e P. Tranquilli . Numerical solution of the Fokker-Planck equation with variable coecients, Journal of Quantitative Spectroscopy & Radiative Transfer 109 , no. 5 (2008) 727-740. [2] Ó. López Pouso e N. Jumaniyazov . Numerical experiments with the Fokker-Planck equation in 1D slab geometry, Journal of Computational and Theoretical Transport 45 , no. 3 (2016) 184-201. [3] Ó. López Pouso e N. Jumaniyazov . Direct versus iterative methods for forwardbackward equations. Numerical comparisons on a particular transport kinetic model , en proceso de publicación. [4] D. Stein e I. B. Bernstein . Boundary value problem involving a simple Fokker- Planck equation, Physics of Fluids 19 (1976) 811-814. [5] V. Vanaja e R. B. Kellogg . Iterative methods for a forward-backward heat equation, SIAM Journal on Numerical Analysis 27 , no. 3 (1990) 622-635. 43