Translating Differential Equation Models Into Kernel ...



Translating Differential Equation Models Into Kernel Methods for Data Analysis – Phase IV Emphasis on Simulations

Presented to:

Drs. William Macready and Ashok Srivastava, NASA

CAMCOS team:

Shikha Naik, Afshin Tiraie,

Virginia Banh, Bao Fang,

Efrem Rensi, Lawrence Varela

Supervising faculty:

Drs. Igor Malyshev and Maria Cayco

Department of Mathematics

San Jose State University

May 2005

CAMCOS

Center for Applied Mathematics and Computer Science

SJSU

Table of Contents

Summary

Introduction

I. Analysis

1. Fundamental solutions for the Fokker-Planck type equations with variable coefficients

2. Generalized Kernel Method for parabolic PDEs.

3. Eigenfunction expansion for the 1-D Fokker-Planck equation.

4. Project background information and resources.

II. Theoretical aspects of the Statistical Analysis performed in the project

1. Mean Square Error (MSE)

2. PDFs in use

3. Maximum Likelihood Inference.

III. Simulation concepts, techniques and outcomes.

IV. The “data generators” and statistical experiments

1. Regular double well Fokker-Planck equation case

2. General double well Fokker-Planck equation case

3. Log-Likelihood and the best fit parametric PDFs

V. Kernel-based Solutions of the Fokker-Planck family of equations

VI. Eigenfunction expansion – numerical experiments and analysis

VII. The “Cloud Evolution” demo

Project extensions

References

Summary

This is the fourth and final phase of the NASA/CAMCOS project on Kernel Methods for Differential Equation Models, initiated by Dr. William Macready at NASA Ames in the fall of 2003.

This semester we investigated a set of four one-dimensional Fokker-Plank equations with single and double well potentials and variable diffusion. Since the kernel method we used required knowledge of the Green’s functions, we investigated simulation techniques that would generate all necessary supporting functions. The kernel method itself has been streamlined, and now it can be applied to an arbitrary one-dimensional parabolic problem with, generally speaking, a non-symmetric and non-positive definite Green’s function or kernel. The outline of the method and the proof of its convergence are provided in the report.

Part of the project consisted of studying reconstruction techniques that can be applied to (hypothetically) incomplete and noisy initial data which is later used in the generation of the kernel-based approximations. Numerous simulations have been conducted under a variety of noise and data loss conditions, using several interpolating options. Statistical data was collected and analyzed, and recommendations were made.

A large number of MATLAB programs have been written (or modified from the previous projects) to simulate the above mentioned procedures, to investigate the properties of the approximate kernel solution(s) to the Fokker-Planck equations, and to perform other tasks. The mathematical treatment, numerical experiments and software developments are presented in this report.

In addition, a known eigenfunction expansion of the solution of the single well Fokker-Planck equation with constant diffusion has been independently derived and some numerical investigations have been conducted.

Introduction

The goal of this project is to continue the investigation of the methods of translating differential equation models into kernel functions. More specifically, the project will explore simulation techniques for parabolic PDEs, modify and streamline the Kernel Method in terms of its application to the general Fokker-Planck equations, and investigate effects of noise and data loss on the construction of the kernel functions that yield good approximation to the solutions of the given equations. The long-term goals of the project include the development of techniques suitable for the models described by

more general multidimensional evolution partial differential equations and systems of such equations.

Kernel methods represent a recent advance in machine learning algorithms, with applications including the analysis of very large data sets. Central to the kernel method is the kernel function, which is essentially a mathematical model that allows for the analysis of regression and for density distribution estimates. It is therefore critical to choose a kernel function that reflects the properties of the physical domain and the corresponding differential equation models.

The current project serves as a fourth step toward the spectrum of physically motivated models involving multi-dimensional domains and equations with variable coefficients, noisy data, etc. These complex models are central to NASA’s long-term goals for this project.

This semester we conducted an investigation of the kernel method for the general Fokker-Planck equation

[pic] (0.1)

with zero boundary conditions. Here the solution [pic] is a probability density function of the position (or velocity) of a particle, [pic] is a constant or variable diffusion coefficient, and [pic] is a drift coefficient function, where [pic] is a single or double well potential of the field of external forces. The solution of the problem is presumed to be known at a finite (sufficiently large) empirical set of points (sample) which is given in the form:

[pic] (0.2)

Our goal is to find an approximation to the solution of (1) in the form of a linear combination of pde-based kernels (Green’s functions measured at the sample):

[pic], (0.3)

that fits the data (2), thus utilizing the techniques of the generic kernel method.

We also investigated effects of noisy and incomplete initial data on the construction of the kernel functions.

We performed the following steps of investigation:

• A new approach has been applied to achieve theoretical justification of the Kernel Method for all cases of the one-dimensional Fokker-Planck equation with constant and variable diffusion coefficient. We believe that this approach can be applied to arbitrary parabolic PDE (sec. I-3).

• Knowledge of the Green’s functions is essential for the Kernel Method application to PDE models. But, since explicit formulae for the PDEs with variable coefficients are seldom available, we pursued an additional line of investigation that consisted of the simulation of the Dirac’s delta used as initial condition for the parabolic PDE of interest (sections I-4, III).

• Several MATLAB programs have been written to conduct a series of numerical experiments for the method validation and property investigation purposes. These programs, in different combinations depending on the task to be performed, have been organized in sets for the user’s convenience (section V, CD).

• We considered the possibility of noisy and “defective” data and designed equation specific “data generators” to implement those effects and perform statistical experiments and analysis to identify the best interpolation techniques. For the Gaussian and pde-based kernel approximation cases we used the Kernel Method mentioned above (section IV).

• We investigated some numerical aspects of the eigenfunction expansions of the solution of the standard one-dimensional Fokker-Planck equation in an infinite domain (section VI).

I. Analysis

1. Fundamental solutions for the Fokker-Planck type equations with variable coefficients

1.1. Fokker-Planck equation case

a) Direct transformation

The Fokker-Planck equation

[pic] (1.1.1)

is a particular case of the equation ([1], p. 93, #2.2):

[pic] (1.1.2)

for the choice of [pic], where the following substitution of variables

[pic] (1.1.3)

[pic] (1.1.4)

reduces (1) to the heat (diffusion) equation with constant coefficients

[pic] (1.1.5)

where

[pic]. (1.1.6)

The verification is quite straightforward but useful, since it shows exact relations between operators that effect some calculations ahead:

[pic] (1.1.7)

Let us denote the operator on the left side of (7) as [pic], then the fundamental solution of the Fokker-Planck equation at the initial moment [pic] corresponding to the instantaneous delta-source at a fixed point [pic], satisfies:

[pic] (1.1.8)

If compared to (7), we obtain

[pic], (1.1.9)

where [pic] and [pic] need to be expressed in terms of [pic] and [pic]:

[pic] (1.1.10)

Thus,

[pic] (1.1.11)

where [pic] is fixed and is not a subject of the substitution performed above.

The fundamental solution of the heat operator (left side of (11)) is

[pic] (1.1.12)

where [pic] is the Heaviside function. Formally, the solution of the non-homogeneous equation (11) can be written as a convolution of [pic] with the term [pic] on the right ([2], [pic] for [pic] and, respectively, [pic] for [pic]):

[pic]

[pic] (1.1.13)

Applying the substitution of variables

[pic] (1.1.14)

to the integrals in (13) and using the property of the delta-function ([pic], twice), we can simplify (13) as follows:

[pic]

[pic]

[pic] (1.1.15)

Thus the solution of (11) has the form

[pic] (1.1.16)

Remark. Another technique of obtaining the same formula for [pic] will be given in the case of a modified Fokker-Planck equation to expand the application base.

b) Inverse transformation.

Going back to the original problem and its solution we recall (3), (4) and (6) to obtain

[pic]. (1.1.17)

Thus we need to see how the return to [pic] will affect the function [pic]. Since

[pic] (1.1.18)

we find [pic]

[pic] (1.1.19)

where [pic] is the fundamental solution of (1) ([3], p.100) at [pic].

Remark. It suffices to replace [pic] with [pic] in the formula (19) to reflect the shift of the initial condition ([pic] in our case) from [pic]moment to [pic]. A straightforward substitution of the time variable in (1) ([pic], etc.) does not affect the equation and the derivation of (19). Therefore, the final formula for the fundamental solution [pic]of (1) takes on the form:

[pic] (1.1.20)

1.2. Modified Fokker-Planck equation case

The modified FP equation has the form (see [4]):

[pic] (1.1.21)

It corresponds to another particular case ([1], p.85) of

[pic] (1.1.22)

with the obvious choice of coefficients.

A substitution of variables similar to the one in part I

[pic] (1.1.23)

brings up again the equation (5) where

[pic]. (1.1.24)

The arbitrary constant [pic] can be chosen to be [pic] which makes [pic] for [pic] (compare to (3)).

In an effort to find a fundamental solution for (21) we shall apply a technique different from the one used in part I. Consider the initial condition for (21) in the form

[pic] (1.1.25)

The substitution of (23) into (24) with [pic] implies:

[pic] (1.1.26)

since [pic] at [pic]. Thus we obtain the following initial-value problem:

[pic] (1.1.27)

And, since the distribution [pic], we finally have

[pic] (1.1.28)

The fundamental solution of (28) is equal to [pic] times the constant [pic], where K is defined in (12). Therefore, following the reverse substitution of variables (18)-(19), and using (24), the fundamental solution of (21) can be written in the form:

[pic]

[pic]

[pic] (1.1.29)

Incidentally,

[pic], (1.1.30)

where [pic] is the fundamental solution of (1) ([3], p.100) at [pic]. Transition to the case of [pic] is identical to part I (see (20)), thus

[pic] (1.1.31)

Remark. It can be proved that [pic] satisfies the modified FP equation in [pic] and the adjoint one in [pic]. It behaves like a delta-function when [pic]. It also decays to [pic] with [pic].

Remark. The connection between [pic] for the modFP and [pic] for the FP can be established by other means ([4]), but the point here was to find [pic] independently of the fact whether those equations were connected or not.

1.3. General Fokker-Planck equation (special case)

Let us consider equation ([1], p. 99, #6):

[pic] , (1.1.32)

which can be classified as a “special case” of a general FP equation in one space variable.

The substitution of [pic] defined by the one-to-one function

[pic] (1.1.33)

and

[pic] (1.1.34)

changes (32) into an equation with constant coefficients

[pic] (1.1.35)

whose fundamental solution (see also (12)) is given by

[pic] (1.1.36)

The initial condition for (32) in the form [pic] implies the corresponding initial condition for (35) as follows:

[pic] (1.1.37)

According to [2], the solution of the problem (35), (37), extended by 0 for [pic], satisfies the non-homogeneous equation

[pic] (1.1.38)

and it can be written in the form

[pic] (1.1.39)

Using the substitution

[pic] (1.1.40)

and the property of the delta-function

[pic] (1.1.41)

transforms (39) into

[pic]

[pic] (1.1.42)

By performing reverse substitution ([pic], (33)), we obtain the fundamental solution for (32) in the form:

[pic], (1.1.43)

where

[pic] (1.1.44)

Remarks. 1) It can be verified that (43) satisfies (32) in [pic] and its adjoint in [pic]. It behaves like delta-function when [pic], and since [pic], it decays to 0 with [pic]. 2) All MATLAB simulations, including the Kernel Method for this case, can be found in the “SpecFPdelta” directory on the enclosed CD.

1.4. Heat equation with variable diffusion coefficient

Let us consider the equation

[pic] (1.1.45)

with the delta initial condition

[pic] (1.1.46)

The substitution of the space variable in the form [pic] [pic] allows reduction of the problem to the one with constant coefficients:

[pic] (1.1.47)

The fundamental solution of the operator in (47) can be found using Fourier transforms in the form:

[pic] (1.1.48)

It can be easily verified to satisfy the properties:

[pic]. (1.1.49)

Using (48), the solution of the problem (47) can be written in the form:

[pic]. (1.1.50)

By reversing the substitution of variables, the fundamental solution of the problem (45)-(46) is:

[pic]. (1.1.51)

Remark. MATLAB simulations for this case, can be found in the “SP_Diff” directory on the enclosed CD.

2. Generalized Kernel Method for parabolic PDEs.

2.1. Given a “data function” [pic]. (Think of [pic]as a continuous function – a “spline” or other interpolation of the data [pic]). We presume it to be a solution to a given PDE at some (unknown) time value [pic]. Then, using the integral representation of the solution via Green’s function of the problem in the form

[pic] (1.2.1)

we, ideally, expect

[pic], (1.2.2)

which is an integral equation of the first kind for the unknown function [pic]. Such problem is a so-called ill posed one. We shall use some regularization technique later in the process to guarantee stability of its (numerically obtained) solution.

Let’s introduce a generic solution of the given PDE in the form (1)

[pic] (1.2.3)

with a yet unknown continuous coefficient function [pic]. We shall try to identify [pic] from the requirement that [pic] in (3) makes a good approximation to the given data function [pic] in (2) for some fixed value of [pic].

Thus, we reduce the solution of (2) to minimizing the following (least squares) functional in the class of functions defined by (3):

[pic], (1.2.4)

Its minimum, using a Frechet differential, can be found from the equation:

[pic], (1.2.5)

where [pic] is from the same class as [pic], that is

[pic] (1.2.6)

Equation (5) can be now put in the form (we shall temporarily drop [pic] from the expressions for [pic] and [pic] for ease of notation):

[pic]

[pic] , (which follows by Fubini’s Theorem)

[pic]

[pic] (1.2.7)

for any [pic]. Since [pic] is dense in [pic] ,

[pic] (1.2.8)

for all [pic] (due to the continuity of the integrand), or in discrete form (a Riemann sum for the uniform partitioning with [pic]fixed):

[pic]. (1.2.9)

Any information about [pic] derived from (9) will be an approximation to that in (8).

Since [pic] is arbitrary, we can replace it with any [pic] from the same partitioning (“sample”) to obtain a system of equations:

[pic] (1.2.10)

Then, introducing matrix [pic] and vectors [pic]

[pic] (1.2.11)

[pic] (1.2.12)

we reduce (10) to the matrix form:

[pic] or [pic] (1.2.13)

Now it is time to recall (3) for [pic]. It allows us to express every component of the vector [pic] in the form

[pic], (1.2.14)

that is,

[pic]. (1.2.15)

And from (13), we get

[pic] (1.2.16)

Let constant [pic] be “absorbed” by the unknown vector of coefficients [pic] thus producing (15)-(16) in the form:

[pic] (1.2.17)

Finding [pic] requires the inversion of the matrix [pic] which is symmetric and positive definite, therefore invertible. Unfortunately it is also close to singular which makes the problem (17) ill posed. To obtain an approximate solution of the matrix equation (17) we shall consider an “epsilon-regularized” version of it first:

[pic] (1.2.18)

It can be proved that [pic] and [pic] respectively, thus delivering an approximate solution to (17) (see sec. 3. below)

2.2. A Different look at (17) – the “Kernel Method connection”

We start with (3) and obtain a discrete (vector) form of it (via discretization of the integral and computing [pic] at prescribed sample set of points [pic] ):

[pic] (1.2.19)

Since the solution is presumed to be given at some time value as a vector [pic], that is [pic], we need to find the vector of coefficients [pic] satisfying the equation

[pic] (1.2.20)

Let us apply a nonzero matrix [pic] to both sides of (20) to produce

[pic], (1.2.21)

that is

[pic] , (1.2.22)

where [pic] is a symmetric and positive definite matrix, [pic] is the vector to be found and [pic] is a transformed data vector. Any solution of (20) is also a solution of (22). But (22) has a unique solution due to the properties of [pic] . Therefore the same is true for (20).

Let us presume for a moment that we found an exact solution of (22). Thus the unique vector [pic] can be used in construction of [pic] in the form

[pic] (which is simply equal to [pic]) (1.2.23)

The advantage of (22) lies in its connection to the generic kernel method. Since (22) is basically an ill posed operator equation, its (numerical) solution requires regularization.

Let us denote the image of the approximate solution [pic] under [pic] as [pic] and require that it delivers a minimum to the regularized empirical functional as in [4]:

[pic] (1.2.24)

where [pic] is the [pic]-based Reproducing Kernel Hilbert Space (RKHS).

Via a Frechet differential technique, retracing exactly the motions of the generic Kernel Method, we can prove that unique minimizer for (24) can be found in the form

[pic] (1.2.25)

with the vector of coefficients being a solution of the equation

[pic] (1.2.26)

which is nothing else but equation (18), thus providing connection between the approach in 2.1. and the Kernel Method.

2.3. Convergence issues.

a) Theorem. ([5], p. 219 “Inverse Operator”) Let [pic] be a bounded linear operator in Banach space [pic]. If [pic] then [pic] is invertible, [pic] is bounded, and [pic]. (1.2.27)

As a result, the following can be established:

Corollary. Let [pic] and [pic] be bounded operators in a Banach space [pic] with [pic] invertible. Then [pic] is invertible if [pic] , and

[pic]. (1.2.28)

b) Consider equation (22) and its regularized (epsilon-perturbed) version (26) together:

[pic] (1.2.29)

We need to prove that [pic] and [pic] when [pic].

Let

[pic]. (1.2.30)

[pic] is invertible since [pic] is symmetric and positive definite. Then, from the Corollary it follows that we need to require [pic], or

[pic] (1.2.31)

Given that condition (31) is met, from (28) we find:

[pic], (1.2.32)

that is, the inverse of [pic] is also a bounded operator.

c) Next, denote

[pic] (1.2.33)

that is, we have a sequence of bounded invertible operators [pic] convergent to a bounded invertible operator [pic]. For their inverses (which we know exist) we find that

[pic]

[pic]

[pic]

Thus we proved that the sequence of inverses [pic] is convergent to [pic] .

d) Since [pic], (1.2.34)

this immediately leads to other inequalities:

[pic] (1.2.35)

and

[pic] (1.2.36)

which imply the expected convergences and approximation error estimates.

e) It may be of interest to see how the convergence [pic] and other properties can be established for the case of bounded linear operators in [pic], that is, all the operators are matrices. In this case, the matrices [pic] and [pic] can be diagonalized: [pic] where [pic] is a diagonal matrix of positive eigenvalues [pic] of [pic] (symmetry and positive definiteness of [pic]). Therefore equations in (29) can be put in the form

[pic] (1.2.37)

and subsequently reduced to

[pic] (1.2.38)

The diagonal matrix [pic] is obviously invertible and

[pic](1.2.39)

The rest follows automatically since [pic] is invertible:

[pic]. (1.2.40)

f) The “C_epsilon_allcases” directory on the attached CD contains MATLAB programs that demonstrate convergence of the Kernel Method by means of computing estimates (35)-(36) and some other related norms for all cases covered by this report (single-well FP both formula based and simulated, double-well FP, and 2 cases of general simulated FP). We assembled, in the illustration, some figures and data for 2 cases: formula-based “swFP” and simulated “dw_generalFP”.

formula-based “swFP” case

[pic]

The figure above shows the behavior of the right side of the estimate (36) in time and epsilon for [pic] and [pic].

The following table contains several time-slices of the data presented in a surface form above:

|epsilon |t_ind = 2 |t_ind = 12 |t_ind = 22 |t_ind = 32 |t_ind = 42 |t_ind = 51 |

| | | | | | | |

|10^ -1 |2.12E+16 |4.62E+18 |5.76E+18 |1.43E+19 |8.94E+18 |1.81E+19 |

|10^ -2 |2.09E+15 |4.58E+17 |5.73E+17 |1.43E+18 |8.92E+17 |1.80E+18 |

|10^ -3 |2.09E+14 |4.57E+16 |5.73E+16 |1.43E+17 |8.92E+16 |1.80E+17 |

|10^ -4 |2.09E+13 |4.57E+15 |5.73E+15 |1.43E+16 |8.92E+15 |1.80E+16 |

|… |2.09E+12 |4.57E+14 |5.73E+14 |1.43E+15 |8.92E+14 |1.80E+15 |

| |2.09E+11 |4.57E+13 |5.73E+13 |1.43E+14 |8.92E+13 |1.80E+14 |

| |2.09E+10 |4.57E+12 |5.73E+12 |1.43E+13 |8.92E+12 |1.80E+13 |

|… |2.09E+09 |4.57E+11 |5.73E+11 |1.43E+12 |8.92E+11 |1.80E+12 |

| |2.09E+08 |4.57E+10 |5.73E+10 |1.43E+11 |8.92E+10 |1.80E+11 |

| |2.09E+07 |4.57E+09 |5.73E+09 |1.43E+10 |8.92E+09 |1.80E+10 |

| |2.09E+06 |4.57E+08 |5.73E+08 |1.43E+09 |8.92E+08 |1.80E+09 |

| |2.09E+05 |4.57E+07 |5.73E+07 |1.43E+08 |8.92E+07 |1.80E+08 |

| |20885 |4.57E+06 |5.73E+06 |1.43E+07 |8.92E+06 |1.80E+07 |

| |2088.5 |4.57E+05 |5.73E+05 |1.43E+06 |8.92E+05 |1.80E+06 |

|… |208.85 |45730 |57302 |1.43E+05 |89159 |1.80E+05 |

| |20.885 |4573 |5730.2 |14277 |8915.9 |18006 |

| |2.0885 |457.3 |573.02 |1427.7 |891.59 |1800.6 |

| |0.20885 |45.73 |57.302 |142.77 |89.159 |180.06 |

| |0.020885 |4.573 |5.7302 |14.277 |8.9159 |18.006 |

| |0.0020885 |0.4573 |0.57302 |1.4277 |0.89159 |1.8006 |

| |0.00020885 |0.04573 |0.057302 |0.14277 |0.089159 |0.18006 |

|… |2.09E-05 |0.004573 |0.0057302 |0.014277 |0.0089159 |0.018006 |

|10^ -23 |2.09E-06 |0.0004573 |0.00057302 |0.0014277 |0.00089159 |0.0018006 |

|10^ -24 |2.09E-07 |4.57E-05 |5.73E-05 |0.00014277 |8.92E-05 |0.00018006 |

|10^ -25 |2.09E-08 |4.57E-06 |5.73E-06 |1.43E-05 |8.92E-06 |1.80E-05 |

|10^ -26 |2.09E-09 |4.57E-07 |5.73E-07 |1.43E-06 |8.92E-07 |1.80E-06 |

|10^ -27 |2.09E-10 |4.57E-08 |5.73E-08 |1.43E-07 |8.92E-08 |1.80E-07 |

|10^ -28 |2.09E-11 |4.57E-09 |5.73E-09 |1.43E-08 |8.92E-09 |1.80E-08 |

|10^ -29 |2.09E-12 |4.57E-10 |5.73E-10 |1.43E-09 |8.92E-10 |1.80E-09 |

|10^ -30 |2.09E-13 |4.57E-11 |5.73E-11 |1.43E-10 |8.92E-11 |1.80E-10 |

With epsilon small enough, [pic] in (36) does not exceed [pic] across all times.

simulated “dw_generalFP” case

Similarly, in the double-well general FP case we find behavior of the right side of the estimate (36) in time and epsilon for [pic] and [pic].

[pic]

and

|epsilon |t_ind = 1 |t_ind = 11 |t_ind = 21 |t_ind = 31 |t_ind = 41 |t_ind = 51 |

| | | | | | | |

|10^ -1 |3.08E-01 |9.91E+16 |4.03E+17 |2.46E+18 |4.22E+18 |1.04E+18 |

|10^ -2 |3.04E-02 |9.81E+15 |4.00E+16 |2.44E+17 |4.19E+17 |1.03E+17 |

|10^ -3 |3.03E-03 |9.80E+14 |3.99E+15 |2.44E+16 |4.18E+16 |1.03E+16 |

|10^ -4 |3.03E-04 |9.80E+13 |3.99E+14 |2.44E+15 |4.18E+15 |1.03E+15 |

|… |3.03E-05 |9.80E+12 |3.99E+13 |2.44E+14 |4.18E+14 |1.03E+14 |

| |3.03E-06 |9.80E+11 |3.99E+12 |2.44E+13 |4.18E+13 |1.03E+13 |

| |3.03E-07 |9.80E+10 |3.99E+11 |2.44E+12 |4.18E+12 |1.03E+12 |

|… |3.03E-08 |9.80E+09 |3.99E+10 |2.44E+11 |4.18E+11 |1.03E+11 |

| |3.03E-09 |9.80E+08 |3.99E+09 |2.44E+10 |4.18E+10 |1.03E+10 |

| |3.03E-10 |9.80E+07 |3.99E+08 |2.44E+09 |4.18E+09 |1.03E+09 |

| |3.03E-11 |9.80E+06 |3.99E+07 |2.44E+08 |4.18E+08 |1.03E+08 |

| |3.03E-12 |9.80E+05 |3.99E+06 |2.44E+07 |4.18E+07 |1.03E+07 |

| |3.03E-13 |9.80E+04 |3.99E+05 |2.44E+06 |4.18E+06 |1.03E+06 |

| |3.03E-14 |9.80E+03 |3.99E+04 |2.44E+05 |4.18E+05 |1.03E+05 |

|… |3.03E-15 |979.51 |3994.4 |2.44E+04 |41845 |1.03E+04 |

| |3.03E-16 |97.951 |399.44 |2440.6 |4184.5 |1030.8 |

| |3.03E-17 |9.7951 |39.944 |244.06 |418.45 |103.08 |

| |3.03E-18 |0.97951 |3.9944 |24.406 |41.845 |10.308 |

| |3.03E-19 |0.097951 |0.39944 |2.4406 |4.1845 |1.0308 |

| |3.03E-20 |0.0097951 |0.039944 |0.24406 |0.41845 |0.10308 |

| |3.03E-21 |0.00097951 |0.0039944 |0.024406 |0.041845 |0.010308 |

|… |3.03E-22 |9.80E-05 |0.00039944 |0.0024406 |0.0041845 |0.0010308 |

|10^ -23 |3.03E-23 |9.80E-06 |3.99E-05 |0.00024406 |0.00041845 |0.00010308 |

|10^ -24 |3.03E-24 |9.80E-07 |3.99E-06 |2.44E-05 |4.18E-05 |1.03E-05 |

|10^ -25 |3.03E-25 |9.80E-08 |3.99E-07 |2.44E-06 |4.18E-06 |1.03E-06 |

|10^ -26 |3.03E-26 |9.80E-09 |3.99E-08 |2.44E-07 |4.18E-07 |1.03E-07 |

|10^ -27 |3.03E-27 |9.80E-10 |3.99E-09 |2.44E-08 |4.18E-08 |1.03E-08 |

|10^ -28 |3.03E-28 |9.80E-11 |3.99E-10 |2.44E-09 |4.18E-09 |1.03E-09 |

|10^ -29 |3.03E-29 |9.80E-12 |3.99E-11 |2.44E-10 |4.18E-10 |1.03E-10 |

|10^ -30 |3.03E-30 |9.80E-13 |3.99E-12 |2.44E-11 |4.18E-11 |1.03E-11 |

The figure below shows epsilon-behavior of the column t_ind = 11.

[pic]

g) The estimates (35)-(36) use operator norms to show convergence in principle, but in reality the processes shows a much better “convergence property” than the estimates. The directories “dwFP_simulations”, “dwFPf0_S05”, “genFPdw200”, “genFPdw_f0”, “genFPsw200”, “genFPsw_f0”, “swFP_simulations”, “swFPf0_S05” each contain MATLAB programs that show convergence of corresponding KM to an “perfect surface” at a specified time-slice. For example, in a double-well FP case at the time-slice with t_ref = 4, we find that for [pic] the quality of fit improves rapidly, and is very reasonable even at[pic]:

[pic]

p = 1, mse = 4.5963e-005 norm(y-f) = 0.0067796

p = 2, mse = 6.3379e-007 norm(y-f) = 0.00079611

p = 3, mse = 6.7219e-009 norm(y-f) = 8.1987e-005

p = 4, mse = 6.9516e-011 norm(y-f) = 8.3376e-006

p = 5, mse = 7.3266e-013 norm(y-f) = 8.5595e-007

p = 6, mse = 1.5989e-014 norm(y-f) = 1.2645e-007

p = 7, mse = 1.2125e-015 norm(y-f) = 3.4821e-008

p = 8, mse = 2.5155e-017 norm(y-f) = 5.0155e-009

p = 9, mse = 2.7725e-019 norm(y-f) = 5.2654e-010

p = 10, mse = 2.8006e-021 norm(y-f) = 5.292e-011

p = 11, mse = 2.8022e-023 norm(y-f) = 5.2935e-012

All quoted cases demonstrate comparable quality of fit and fast convergence.

3. Eigenfunction expansion for the 1-D Fokker-Planck equation - derivation

Consider standard 1-D Fokker-Planck equation in the infinite domain with natural boundary conditions ([pic] at [pic] ):

[pic] (1.3.1)

A substitution of variables

[pic] (1.3.2)

reduces (1) to

[pic] (1.3.3)

Another change of functions

[pic] (1.3.4)

leads to a modified equation:

[pic] (1.3.5)

where

[pic] (1.3.6)

Now, applying the separation of variables routine to (5)

[pic], (1.3.7)

we obtain the following eigenvalue problem:

[pic] (1.3.8)

If we presume zero boundary conditions at infinity, then using [6], we find a complete system of corresponding eigenvalues/eigenfunctions is given by:

[pic] [pic] (1.3.9)

[pic] (1.3.10)

where [pic] are the Chebyshev/Hermite polynomials.

[pic] is an orthonormal set of functions:

(1.3.11)

Therefore,

[pic], (1.3.12)

[pic] (1.3.13)

and

[pic] (1.3.14)

Thus, using (4), we found the eigenfunction expansion for the solution of (3), satisfying 0 boundary conditions at infinity:

[pic] (1.3.15)

Now, if we impose the initial condition function in (3)

[pic], (1.3.16)

the (Fourier) coefficients [pic] can be found as follows:

[pic] (1.3.17)

or

[pic] (1.3.17’)

To find an expression for the fundamental solution (Green’s function) of (3) with 0 b.c. we put

[pic]. (1.3.18)

And so,

[pic] (1.3.19)

Therefore, the eigenfunction expansion formula for the corresponding Green’s function is:

[pic] (1.3.20)

or

[pic] , (1.3.20’)

or even

[pic] (1.3.20’’)

depending on computational preferences.

The original Green’s function, W, from the original equation (1) can be calculated as follows. Since

[pic], (1.3.21)

using (15) we can find

[pic] (1.3.22)

and

[pic] (1.3.23)

It can be easily verified that

[pic] (1.3.24)

Then, similarly to (17), we find

[pic] (1.3.25)

If the initial condition for (1) is

[pic] (1.3.26)

then

[pic] (1.3.27)

and

[pic]

[pic] (1.3.28)

2. Remarks

1) A “shortcut” connection between [pic] and [pic] can be established as follows. We start with the integral representation of the solution [pic] in the form

[pic] (1.3.29)

and follow with the substitution (2) to obtain

[pic].

Now let, [pic]. Then [pic]. And we get,

[pic]. (1.3.30)

Thus,

[pic] (1.3.31)

as above (compare (20) and (28)).

2) Although the expansion of the type (28) is mentioned in quantum mechanics literature ([7]), we offer here another (more transparent) derivation of it.

3. Validation.

Starting with a known series ([8], p.329, case m = 0)

[pic] (1.3.32)

we replace [pic] with [pic] to find

[pic] , (1.3.33)

[pic]

[pic]

[pic], (1.3.34)

and, finally

[pic] =

[pic]

[pic]

[pic]

[pic]. (1.3.35)

That is, [pic] in (28) is identical to the expression of the fundamental solution (Green’s function) of the FP equation in [3].

4. Project background information and resources.

4.1. Equations

1) General Fokker-Planck equation

[pic] (1.4.1)

[pic] (1.4.2)

or

[pic] (1.4.3)

Single-well parabolic potential

[pic] (1.4.4)

Asymmetric double-well parabolic potential

[pic] (1.4.5)

where

[pic] (1.4.6)

[pic] is continuous at 0 if constants [pic] are chosen from the relation:

[pic] (1.4.7)

Then

[pic] (1.4.8)

([pic] is discontinuous at 0).

The standard FP corresponds to the choice of constant [pic]:

[pic] (1.4.9)

Further restrictions on coefficients (constant potential, therefore [pic]) reduces FP to a diffusion equation:

[pic] (1.4.10)

Remark. Depending on the choice of the potential [pic] and diffusion coefficient [pic] we shall have 4 cases to consider.

Special FP case requires special combination of coefficients:

[pic] (1.4.11)

Special diffusion case (optional):

[pic] (1.4.12)

4.2. All problems will require initial condition

[pic] (1.4.13)

and 0 boundary conditions at [pic] or at infinity. (1.4.14)

4.3. Known and derived fundamental solutions

diffusion equation (10):

[pic] (1.4.15)

standard FP equation (9) with a single-well potential (4):

[pic] (1.4.16)

3) special FP equation (11):

[pic] (1.4.17)

where

[pic] (1.4.18)

and

[pic] (1.4.19)

special diffusion case (optional):

[pic] (1.4.20)

Remark. In all cases, we shall try a simulation technique to generate both known fundamental solutions (for comparison/validation) and unknown ones (as the only source of information)

4.4. Integral representation via fundamental solution/Green’s function

The solution of all the problems outlined in 4.1-4.2 above can be expressed in the form:

[pic] (1.4.21)

(initial-value problem in infinite domain), or

[pic] (1.4.22)

(initial-value/boundary-value problem), where [pic] stands for either fundamental solution in case of (21), or Green’s function in case of (22).

These expressions make a “backbone” of the Kernel Method in application to PDEs.

Simulation of the delta-function and Green’s functions

Without entering the detail of the so-called generalized functions/distributions and their spaces, we define here delta sequence as a “cap-shaped function” with the property:

[pic], (1.4.23)

which means that [pic] in the distributional sense. Different examples of [pic] are available in the literature. We shall use the following ([2]):

[pic] (1.4.24)

where constant [pic] is chosen so that

[pic] (1.4.25)

Green’s functions for the problems in 4.1-4.2 can be found as solutions to corresponding PDEs with [pic] as initial condition, thus representing a response of the system at point [pic] to instantaneous unit delta source (disturbance) applied at the point[pic] . For the simulation purposes we shall use delta-sequence [pic] (instead of delta) as the initial condition and the pdepe Matlab solver.

II. Theoretical aspects of the Statistical Analysis performed in the project

1. Mean Square Error (MSE)

Let [pic] and [pic] be the empirical distribution and the hypothesized pdf, respectively. Then, the mean square error measures the mean of the summation of the differences at each point between [pic] and [pic].

MSE = [pic]

where n is the number of points and k is the number of parameters that where estimated. In our case, since n is very large, k and 1 can be ignored in the denominator.

We will perturb the original data by introducing noise or data reduction or both, then we will reconstruct it. We will perform this process 1000 times for each case. Then, we will analyze which type of distribution the MSEs

2. Probability Density Functions in use

Here are the distributions of the mean square error data after experiment of running 1000 times.

a. The Inverse Gaussian Distribution

A continuous random variable X is said to have an Inverse Gaussian Distribution if the pdf of X is:

[pic] where x > 0

= 0 [pic] otherwise

where the parameter [pic] is the mean.

[pic]

Figure 1-The Inverse Gaussian Function with Noise Introduced

b. Log logistic

The variable x has a log logistic distribution with location parameter µ and scale parameter [pic] > 0 if ln x has a logistic distribution with parameters µ and[pic]. The logistic distribution has the density function:

[pic] when x [pic] 0

= 0 otherwise

with location parameter µ and scale parameter [pic] > 0, for all real x.

[pic]

Figure 2-The Log-Logistic Function with Both Effects

c. Log-normal distribution

Let a random variable X be normally distributed with mean [pic] and variance[pic]. Then if we write X = ln Y, then Y is said to have a Log-Normal distribution. Then the pdf for Y is

[pic]

For a Log-Normal distribution Y, [pic] and [pic]

[pic]

Figure3-Log Normal Distribution

d. Birnbaum-Saunders distribution

The Birnbaum-Saunders distribution has the density function:

[pic] when x >0

= 0 otherwise

with scale parameter [pic] > 0 and shape parameter [pic] > 0, for x > 0. If x has a Birnbaum-Saunders distribution with parameters [pic] and[pic], then

[pic]

has a standard normal distribution.

[pic]

Figure4-The Birnbaum Saunders Distribution

3. Log Likelihood Inference

Let specify a probability density or probability mass function for the observations

[pic]

In this expression θ represents one or more unknown parameters that govern the distribution of Z.

Let Z has any distribution, then the likelihood function, L(θ; Z), is the product of the probability density function [pic]evaluated at the n data points.

[pic]

the probability of the observed data under the model [pic], and we can consider that [pic] as a function of θ with Z be fixed.

We denote the logarithm of [pic] by[pic], and this expression is called the log-likelihood, and each value [pic]is called a log-likelihood component.

III. Simulation concepts, techniques and outcomes.

1. Omega-epsilon Simulation

1. Discussion

The Green’s functions for the problems in 4.1-4.2 can be found as solutions to corresponding PDEs with [pic] as its initial condition, thus representing a response of the system at a point [pic] to instantaneous unit delta source (disturbance) applied at the point[pic]. As stated before, the set of equations for which we generated Green’s functions were the Fokker-Planck family of equations. The Fokker-Planck equations are defined with a drift [pic] and diffusion [pic] coefficients that influence the predicted distribution of particles in an external field.

In our experiments we simulated the Green’s functions with single and double well potentials f(x) (potential of the field of external forces resisting the motion of particles) whose profiles are illustrated below.

Single-well Double-Well

[pic][pic]

In order to initiate our simulation we needed to approximate the initial source function which is a Dirac’s delta. Since delta functions can not be used by MATLAB directly, we replaced it with a sequence of functions,[pic], which is a smooth substitute which has a spike at one point.

1.2 Experiment

We defined [pic] as a “cap-shaped” sequence of functions, by the formulas below:

It was especially important that we chose our constant [pic] in such a way that guaranteed our integration of [pic] to produce a value of 1. This requirement of maintaining an integration of 1 stems from our goal of maintaining probabilistic properties throughout our simulations. The following code describes the implementation we used to define our Omega-epsilon with the requirements stated above.

function u0 = swfp_deltaic(x)

% ic - initial cond. at t=0

%

global xp J

p = -5;

e = 10^p;

xJ = xp(J);

if abs(x-xJ) ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download