Wedig, Walter

English

Repository KITopen

On the optimization of non-negative density solutions of stationary Fokker-Planck equations

KITopen

On the Optimization of Non-Negative DensitySolutions of Stationary Fokker-Planck EquationsWalter V. Wedig ∗Universiẗat Karlsruhe, D-76128 Karlsruhe, GermanyClassical solution methods of Fokker-Planck equations (see e.g. [1],[2],[3]) lead to negative densities at least in the tales wherethe density solutions strongly approach zero. To avoid these inconsistencies the paper proposes to apply central differencesschemes in combination with large-scaled quadratic optimization programming [4].1 Introduction to the problemStochastic dynamical systems are investigated in the time domain by direct path-wise integration along the process trajectoriesor indirect in the space domain by solving associated Fokker-Planck equations to determine relative frequencies or probabilitiesof the processes. To test details of integration, consider the example of a linear scalar system described by the Itô equationdXt = −ωgXtdt + σdWt, E(dW 2t ) = dt. (1)Herein,ωg is a limiting frequency of the low-pass processXt andσ is the intensity parameter of its excitation by the incrementsdWt of the Wiener processWt. It is suitable to introduce dimensionless time and state byXt = X̄t√E(X2t ), with E(X2t ) = σ2/(2ωg), (2)dτ = ωgdt, dWτ =√ωg dWt. (3)Theintroduction of both dimensionless quantities into the system equation (1) leads to the parameter free equationdX̄τ = −X̄τdτ +√2dWτ , E(dW 2τ ) = dτ. (4)−4 −3 −2 −1 0 1 2 3 400.050.10.150.20.250.30.350.4x →↑ p(x)−1.5 −1 −0.5 0 0.5 1 1.500.10.20.30.40.50.60.70.8ϕ →↑ p(ϕ) λ = .5, 1, 2Fig. 1 a) Normalized Gaussian density. b) Density transformed byϕ = arctan λx.It possesses the Fokker-Planck equation (5) which is easily integrated under the stationarity condition∂p/∂τ = 0.∂p(x, τ)∂τ= − ∂∂x[xp(x, τ)]− 12∂2∂x2[2p(x, τ)], (5)p(x) =1√2πexp(−12x2), for |x| ≤ ∞. (6)Herein,p(x) is the well-known stationary density solution of the Fokker-Planck equation (5). The left picture in figure 1shows the numerical evaluation ofp(x) in the range|x| ≤ 4.∗ Correspondingauthor: e-mail:wedig@itm.uni-karlsruhe.de, Phone: 0049 721 608 2658, Fax: 0049 721 608 6070PAMM · Proc. Appl. Math. Mech. 4, 23–26 (2004) / DOI 10.1002/pamm.200410006  © 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim © 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim First published in: EVA-STAR (Elektronisches Volltextarchiv – Scientific Articles Repository) http://digbib.ubka.uni-karlsruhe.de/volltexte/10000040552 Transformation to finite solution rangesUnder the stationarity condition∂p/∂τ = 0 the Fokker-Planck equation (5) leads to an ordinary differential equation whichhas to be solved under the two boundary conditionsp(x = 0) = C andp(x = ±∞) = 0. Both are necessary to fulfillthe normalization condition that the integration ofp(x) over the entire range|x| ≤ ∞ is one. Because of the homogeneityof the Fokker-Planck equation the constant C represents a simple factorization and can be chosen to one without any loss ofgenerality. Vanishing density boundaries onx = ±∞ are necessary to satisfy the radiation conditions that the solution processremains inside the range|x| ≤ ∞. To investigate these singular boundaries numerically, the transformationφ = arctan λxis introduced with0 < λ < ∞. Provided thatλ is a free parameter with finite values, thearctan-transformation reduces thesingular density range|x| ≤ ∞ to the finite one|ϕ| ≤ π/2 leading to the transformed density distributionp(ϕ) =1λ√2π cos2 ϕexp(− tan2 ϕ2λ2), for |ϕ| ≤ π/2, (7)It is shown in the right side of figure 1 for the three parameter valuesλ = .5, 1 andλ = 2. Associated to the transformeddensity, above, there is a transformed Itô equation of the formdΦτ = −(1 + 2λ2 cos2 Φτ ) sin Φτ cosΦτdτ +√2λ cos2 ΦτdWτ , (8)which is derived applying It̂o’s calculus [5] to the angle processΦτ = arctan λX̄τ . The transformed stochastic differentialequation (8) leads to the Fokker-Planck equation∂∂ϕ[(1λ2+ 2 cos2 ϕ) sin ϕ cosϕp(ϕ)] +∂2∂ϕ2[cos4 ϕp(ϕ)] = 0, (9)when the stationary density solution is of interest only.−1.5 −1 −0.5 0 0.5 1 1.500.10.20.30.40.50.60.70.8ϕ →↑ p(ϕ)N = 10−1.5 −1 −0.5 0 0.5 1 1.500.10.20.30.40.50.60.70.8ϕ →↑ p(ϕ) λ = .5, 1, 2N = 23Fig. 2 Discrete and exact solutions of the Fokker-Planck equation (10) for the density parametersλ = .5, 1, 2.3 Application of central differences schemesThe stationary Fokker-Planck equation (9) of the transformed system can once be integrated leading to the first order form(1λ2− 2 cos2 ϕ) sin ϕ cos ϕp(ϕ) + cos4 ϕ ∂∂ϕ[p(ϕ)] = C. (10)Herein,C is a constant of integration which is vanishing when symmetry respectively radiation conditions have to be satisfied.Subsequently, one cosine function in (10) can be cancelled out. Applying the central difference scheme bydp/dϕ ≈ (pn+1 −pn−1)/(2∆ϕ) with the step size∆ϕ = π/(2N), the discrete form of the Fokker-Planck equation (10) is derived toanpn + c3npn−1 − c3npn+1 = 0, cn = cos n∆ϕ, (11)an = (2c2n −1λ2)2sn∆ϕ, sn = sin n∆ϕ. (12)It is valid for all n = 1, 2, . . . , N − 1. Provided that the density solution is symmetric, the following boundary conditions areto be satisfied:p0 = 1, (Normalization), pN = 0. (Radiation Condition) (13)Minisymposium MA2 24  © 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Obviously, this discrete formulation leads toN − 1 inhomogeneous equations to be evaluated for givenλ−parameter and stepsize∆ϕ. Figure 2 shows evaluations for the parametersλ = .5, 1, 2 and forN = 10 steps on the left picture andN = 23 onthe right-hand side. Correspondingly, the discrete solutions are stable showing good coincidences with the exact solutions forlow parameter valuesλ and small step sizes∆ϕ. For high parametersλ, the step numberN has to be sufficiently increased toobtain a smooth density close to the exact solution.4 Application of central differences schemesPerforming all differentiations with respect toϕ, the second order Fokker-Planck equation (9) takes the formc4d2pdϕ2− 6sc3 dpdϕ− 2c2(c2 − 3s2)p + 1λ2[scdpdϕ+ (c2 − s2)p] = 0, |ϕ| ≤ π/2. (14)Herein,c ands are abbreviations forc = cos ϕ ands = sin ϕ, respectively. Applying the central differences scheme with thestep size∆ϕ = π/(2N) leads to the discrete version of (14)αnpn + βnpn+1 + γnpn−1 = 0, n = 1, 2, . . . , N − 1. (15)This recursion formula possesses the coefficientsα , βn andγn calculated toαn = 2c4n + 2c2n(c2n − 3s3n)∆ϕ2 − (c2n − s2n)∆ϕ2/λ2, sn = sin n∆ϕ, (16)βn(γn) = −c4n ± 3c3nsn∆ϕ∓ cnsn∆ϕ/(2λ2), cn = cos n∆ϕ. (17)Figure 3 shows evaluations of (15) together with the exact solutions (7) forN = 10 andλ = .5 in the left picture, respectivelyfor N = 40 andλ = 2 in the right one. The left sides in both pictures are calculated under the same boundary conditionsp0 = 1 andpN = 0, as previously applied. The right sides in both pictures show solutions which satisfyp0 = 1 andp1 = p−1.Obviously, these boundary conditions are less optimal. The associated solutions diverge strongly on the right-hand side at theboundaryn = N andϕ = π/2, respectively.−1.5 −1 −0.5 0 0.5 1 1.500.10.20.30.40.50.60.70.8ϕ →↑ p(ϕ)N = 10λ = 0.5−1.5 −1 −0.5 0 0.5 1 1.5−0.2−0.100.10.20.30.40.50.60.7ϕ →↑ p(ϕ) λ = 2N = 40Fig. 3 Discrete and exact solutions of the second order Fokker-Planck equation (14) forλ = 0.5 and2.0.5 Large-scale quadratic optimization programming.For further improvements higher order schemes may be applied or more efficiently large-scaled quadratic optimization to avoidnegative results in the density tales. For these purposes the scalar equations (15) are written into the vector formAp = b.α1 β1 0 . . . 0γ2 α2 β2 . . . 0......0 0 . . . γN−1 αN−1p1p2...pN−1 =−γ1p00...−βN−1pN . (18)Solving this system under the boundary conditionsp0 = 1 andpN = 0, the firstNeq-density values, obtained, are positive.Only in the tale range forNeq < n < N , there are negative density values beginning withpNeq < 0. For the special exampleN = 24 andλ = .5, 1, 2 shown on the left side of figure 4, these numbers areNeq = 17, 20 and22. Naturally, negative densityPAMM · Proc. Appl. Math. Mech. 4 (2004)  25  © 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim values are physically meaningless and can only be explained by the approximating procedure in the tale range where densitysolutions strongly approach zero. To avoid these inconsistencies, it is proposed to apply large-scale quadratic programming,given e.g. in Matlab 6.1. Accordingly, the best approximation of the equation systemAp = b under the constraint conditionsp ≥ 0, (N-1 lower bounds), Aeqp = beq, (Neq equality constraints) (19)is calculated by minimizing the quadratic equation error(Ap− b)T (Ap− b) = Min! , ⇒ pT AT Ap− 2bT Ap = Min! (20)Typical numerical results of this minimizing procedure are shown on the right-hand side of figure 4. Particularly forλ = 2 andN = 24 the optimization results are still good. In contrast to the classical results on the left side of figure 4 there are no negativedensity values on the right-hand side. Of course, the non-negative density values, obtained by the quadratic optimization,converge towards the exact solutions shown in figure 4 for increasingN and decreasing step size∆ϕ, respectively.−1.5 −1 −0.5 0 0.5 1 1.500.10.20.30.40.50.60.70.8 ↑ p(ϕ)λ = .5, 1, 2neq=17, 20, 22N = 24ϕ →Fig. 4 Discrete, exact and optimized solutions of the Fokker-Planck equation (14).6 Summary and concluding remarks.To solve Fokker-Planck equations by central differences schemes with optimally selected boundary conditions, the infinitesolution range is reduced to a finite one utilizing suitable mappings. The paper investigates large-scale quadratic programmingto avoid negative probability densities in the tales. Further improvements can be achieved by higher order differences schemesor by adaptive finite elements.References[1] S. F. Wojtkiewicz, L. A. Bergman, E. A. Johnson, B. F. Spencer Jr., M. C. Grigoriu, Stochastic Response to Additive Gaussian andPoission White Noises, in: Stochastic Structural Dynamics, Spencer & Johnson (eds), Balkema, Rotterdam, 53-60, (1999).[2] U. v. Wagner & W. V. Wedig. On the Calculation of Stationary Solutions of Multi-Dimensional Fokker-Planck Equations by OrthogonalFunctions, in: Nonlinear Dynamics21, 3, 289-306, (2000).[3] A. Naess & V. Moe, New Techniques for Path Integral Solution of the Random Vibration of Non-Linear Oscillators, in: ProceedingsIcossar ’97, Kyoto, Japan, Rotterdam, Balkema (1997).[4] T. F. Coleman and Y. Li, A Reflective Newton Method for Minimizing a Quadratic Function Subject to Bounds on some of theVariables, Siam Journal on Optimization,6, 4, 1040-1058, (1996).[5] L. Arnold, Random Dynamical Systems, (Berlin, Springer, 1998).Minisymposium MA2 26  © 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 

https://publikationen.bibliothek.kit.edu/1000004055/774183

On the optimization of non-negative density solutions of stationary Fokker-Planck equations

Abstract

Similar works

Full text

Available Versions

Repository KITopen

KITopen