^{1,2,a)}, Bin Min

^{1,b)}and Zhiming Wang

^{1,c)}

### Abstract

The stochastic integral ensuring the Newton-Leibnitz chain rule is essential in stochastic energetics. Marcus canonical integral has this property and can be understood as the Wong-Zakai type smoothing limit when the driving process is non-Gaussian. However, this important concept seems not well-known for physicists. In this paper, we discuss Marcus integral for non-Gaussian processes and its computation in the context of stochastic energetics. We give a comprehensive introduction to Marcus integral and compare three equivalent definitions in the literature. We introduce the exact pathwise simulation algorithm and give the error analysis. We show how to compute the thermodynamic quantities based on the pathwise simulation algorithm. We highlight the information hidden in the Marcus mapping, which plays the key role in determining thermodynamic quantities. We further propose the tau-leaping algorithm, which advance the process with deterministic time steps when tau-leaping condition is satisfied. The numerical experiments and its efficiency analysis show that it is very promising.

The authors acknowledge the support from the National Science Foundation of China (NSFC) under Grant Nos. 11171009 and 91130005 and the National Science Foundation for Excellent Young Scholars (Grant No. 11222114).

I. INTRODUCTION

II. MARCUS INTEGRAL

A. Taylor series formulation by Di Paola and Falsone

B. Ordinary differential equation (ODE) formulation through Marcus mapping

C. Equivalence between the * integral and Marcus integral

D. General formulations

III. PATHWISE SIMULATION ALGORITHM

A. Algorithm and its convergence analysis

B. Numerical results

IV. COMPUTATION OF THERMODYNAMIC QUANTITIES

A. Computational strategy

B. The first law of thermodynamics in an overdamped Langevin equation

C. Heat measurement formula

V. TAU-LEAPING ALGORITHM

A. Algorithm construction

B. Tau-leaping condition

C. Efficiency analysis

VI. NUMERICAL RESULTS BY TAU-LEAPING METHOD

A. Random motion near two parallel walls

B. Langevin equation with double well potential

C. Langevin equation with periodic forcing

D. High dimensional case with multiplicative noise

VII. CONCLUSION

### Key Topics

- Poisson's equation
- 28.0
- Random noise
- 11.0
- Stochastic processes
- 6.0
- Brownian motion
- 4.0
- Chemical kinetics
- 4.0

## Figures

The smoothed function *L* _{δ}(*t*) of the single-jump realization *L*(*t*).

The smoothed function *L* _{δ}(*t*) of the single-jump realization *L*(*t*).

Shown in the figure is one specific realization of the solution. The blue solid line is the theoretical solution. The star symbol shows the numerical result by pathwise simulation algorithm and the red circle corresponds to the Ito integral. The second order Runge-Kutta methods are used for both drift and jump ODEs and the time stepsize δ*t* = 0.01. Some other parameters are λ = 20 and σ = 0.2.

Shown in the figure is one specific realization of the solution. The blue solid line is the theoretical solution. The star symbol shows the numerical result by pathwise simulation algorithm and the red circle corresponds to the Ito integral. The second order Runge-Kutta methods are used for both drift and jump ODEs and the time stepsize δ*t* = 0.01. Some other parameters are λ = 20 and σ = 0.2.

For ODEs (33) and (34) , we use the Runge-Kutta methods with *p* = 3 and *q* = 3. The time stepsize is chosen to be Δ*t* = 0.006, 0.005, 0.004, 0.003, and 0.002. Some parameters are λ = 20, σ = 0.2, *T* = 1. Three thousand samples are simulated. Shown in the left panel is the log-log plot of the error of mean versus time stepsizes, which gives numerical order 1.0157. The right panel shows the linear fitting of error against stepsize curve, which gives slope 4.1905.

For ODEs (33) and (34) , we use the Runge-Kutta methods with *p* = 3 and *q* = 3. The time stepsize is chosen to be Δ*t* = 0.006, 0.005, 0.004, 0.003, and 0.002. Some parameters are λ = 20, σ = 0.2, *T* = 1. Three thousand samples are simulated. Shown in the left panel is the log-log plot of the error of mean versus time stepsizes, which gives numerical order 1.0157. The right panel shows the linear fitting of error against stepsize curve, which gives slope 4.1905.

The red circle corresponds to the heat based on Marcus integral and the blue solid line corresponds to the total energy *U*. Some parameters are set as λ = 5, λ*I* ^{2} = 1, ε = 0.1 and the time stepsize δ*s* = 0.01 for solving both drift and jump ODEs.

The red circle corresponds to the heat based on Marcus integral and the blue solid line corresponds to the total energy *U*. Some parameters are set as λ = 5, λ*I* ^{2} = 1, ε = 0.1 and the time stepsize δ*s* = 0.01 for solving both drift and jump ODEs.

The red circle corresponds to the numerical heat *Q* and the blue solid line corresponds to the theoretical energy *E* which equals to *Q*. *L*(*t*) is a compound Poisson process with rate λ = 10 and *P*(*x* → *x* ± *I*) = λ/2, where *I* satisfies λ*I* ^{2} = 1. The stepsizes δ*s* = 0.01 for both drift and jump ODEs.

The red circle corresponds to the numerical heat *Q* and the blue solid line corresponds to the theoretical energy *E* which equals to *Q*. *L*(*t*) is a compound Poisson process with rate λ = 10 and *P*(*x* → *x* ± *I*) = λ/2, where *I* satisfies λ*I* ^{2} = 1. The stepsizes δ*s* = 0.01 for both drift and jump ODEs.

The distribution of *X* at time *T* = 10. The symbol * shows the distribution obtained by tau-leaping algorithm and ○ shows the distribution obtained by pathwise simulations. Some parameters are λ = 400, *A* = 500, and *X*(0) = *B*/2 = 500. The time step is δ*t* = 0.01 for pathwise simulations and Δ*t* = 0.06 for tau-leaping method. Five thousand samples are simulated.

The distribution of *X* at time *T* = 10. The symbol * shows the distribution obtained by tau-leaping algorithm and ○ shows the distribution obtained by pathwise simulations. Some parameters are λ = 400, *A* = 500, and *X*(0) = *B*/2 = 500. The time step is δ*t* = 0.01 for pathwise simulations and Δ*t* = 0.06 for tau-leaping method. Five thousand samples are simulated.

Comparison between the pathwise and tau-leaping simulations with adaptive stepsize for double well potential. (a) Plot of the double well potential and (b) the distribution of *X* at time *T* = 10. Shown with blue stars is the distribution obtained by tau-leaping algorithm and red circles is the distribution obtained by pathwise algorithm. Five thousand realizations are simulated.

Comparison between the pathwise and tau-leaping simulations with adaptive stepsize for double well potential. (a) Plot of the double well potential and (b) the distribution of *X* at time *T* = 10. Shown with blue stars is the distribution obtained by tau-leaping algorithm and red circles is the distribution obtained by pathwise algorithm. Five thousand realizations are simulated.

Comparison between the pathwise and tau-leaping simulations with adaptive switching for periodic forcing. (a) The distribution of *X* at time *T* = 10. Shown with blue stars is the distribution obtained by tau-leaping algorithm and red circles is the distribution obtained by pathwise algorithm. Five thousand realizations are simulated. (b) Time history of utilized stepsizes by adaptive tau-leaping algorithm. The red horizontal line corresponds to the threshold Δ*t* _{0}. The blue curve corresponds to the stepsizes used in tau-leaping steps while the green dots corresponds to the stepsizes used in pathwise simulation steps.

Comparison between the pathwise and tau-leaping simulations with adaptive switching for periodic forcing. (a) The distribution of *X* at time *T* = 10. Shown with blue stars is the distribution obtained by tau-leaping algorithm and red circles is the distribution obtained by pathwise algorithm. Five thousand realizations are simulated. (b) Time history of utilized stepsizes by adaptive tau-leaping algorithm. The red horizontal line corresponds to the threshold Δ*t* _{0}. The blue curve corresponds to the stepsizes used in tau-leaping steps while the green dots corresponds to the stepsizes used in pathwise simulation steps.

The distribution of *X* _{3} at time *T* = 10. The symbol * shows the distribution obtained by tau-leaping algorithm and ○ shows the distribution obtained by pathwise simulations.

The distribution of *X* _{3} at time *T* = 10. The symbol * shows the distribution obtained by tau-leaping algorithm and ○ shows the distribution obtained by pathwise simulations.

## Tables

For ODEs (33) and (34) , we use *p*th and *q*th order Runge-Kutta methods, respectively. The time stepsize is chosen to be Δ*t* = 0.01, 0.008, 0.005, 0.004, 0.002, and 0.001. Some parameters are λ = 20, σ = 0.2, and *T* = 1. Two thousand samples are simulated. The numbers shown in the table are the slope by linear fitting compared with the theoretical value in parenthesis.

For ODEs (33) and (34) , we use *p*th and *q*th order Runge-Kutta methods, respectively. The time stepsize is chosen to be Δ*t* = 0.01, 0.008, 0.005, 0.004, 0.002, and 0.001. Some parameters are λ = 20, σ = 0.2, and *T* = 1. Two thousand samples are simulated. The numbers shown in the table are the slope by linear fitting compared with the theoretical value in parenthesis.

*t* _{ r }, *t* _{ f }, *t* _{ g }, and *t* _{ total } are the time for generating random variables, solving the drift ODE (33) , solving the jump ODEs (53) and (54) , and the total computation time, respectively. Mean and Std are the mean value and standard deviation of *X*(*T* = 10). Five thousand simulations are performed. The parameters are the same as those in Fig. 6 .

*t* _{ r }, *t* _{ f }, *t* _{ g }, and *t* _{ total } are the time for generating random variables, solving the drift ODE (33) , solving the jump ODEs (53) and (54) , and the total computation time, respectively. Mean and Std are the mean value and standard deviation of *X*(*T* = 10). Five thousand simulations are performed. The parameters are the same as those in Fig. 6 .

Comparison between the pathwise and tau-leaping simulations with adaptive stepsize for double well potential. The computation time is for 5000 simulations. Mean and Std are the mean value and standard deviation of *X*(10).

Comparison between the pathwise and tau-leaping simulations with adaptive stepsize for double well potential. The computation time is for 5000 simulations. Mean and Std are the mean value and standard deviation of *X*(10).

Comparison between the pathwise and adaptive tau-leaping simulations with different thresholds *r* _{0} for periodic forcing. The computation time is for 5000 simulations. Mean and Std are the mean value and standard deviation of *X*(10).

Comparison between the pathwise and adaptive tau-leaping simulations with different thresholds *r* _{0} for periodic forcing. The computation time is for 5000 simulations. Mean and Std are the mean value and standard deviation of *X*(10).

Comparison between the pathwise and tau-leaping simulations with adaptive stepsize for high dimensional equations. The computation time is for 5000 simulations. In each column of *X* _{ i }, the mean value and standard variation of *X* _{ i }(10) are listed.

Comparison between the pathwise and tau-leaping simulations with adaptive stepsize for high dimensional equations. The computation time is for 5000 simulations. In each column of *X* _{ i }, the mean value and standard variation of *X* _{ i }(10) are listed.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content