update documentation

Author: ytdHuang

docs/src/users_guide/steadystate.md

@@ -94,7 +94,7 @@ tlist = LinRange(0, 50, 100)
 
 # monte-carlo
 sol_mc = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=100, progress_bar=false)
-exp_mc = real(sol_mc.expect[1, :])
+exp_mc = real(average_expect(sol_mc)[1, :])
 
 # master eq.
 sol_me = mesolve(H, ψ0, tlist, c_ops, e_ops=e_ops, progress_bar=false)

docs/src/users_guide/time_evolution/mcsolve.md

@@ -70,6 +70,7 @@ c_ops = [sqrt(0.1) * a]
 e_ops = [a' * a, σm' * σm]
 
 sol_500 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops)
+expect_avg = average_expect(sol_500)
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
@@ -78,8 +79,8 @@ ax = Axis(fig[1, 1],
     ylabel = "Expectation values",
     title = "Monte Carlo time evolution (500 trajectories)",
 )
-lines!(ax, times, real(sol_500.expect[1,:]), label = "cavity photon number", linestyle = :solid)
-lines!(ax, times, real(sol_500.expect[2,:]), label = "atom excitation probability", linestyle = :dash)
+lines!(ax, times, real(expect_avg[1,:]), label = "cavity photon number", linestyle = :solid)
+lines!(ax, times, real(expect_avg[2,:]), label = "atom excitation probability", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 
@@ -104,9 +105,9 @@ ax = Axis(fig[1, 1],
     ylabel = "Expectation values",
     title = "Monte Carlo time evolution",
 )
-lines!(ax, times, real(sol_1.expect[1,:]), label = "1 trajectory", linestyle = :dashdot)
-lines!(ax, times, real(sol_10.expect[1,:]), label = "10 trajectories", linestyle = :dash)
-lines!(ax, times, real(sol_100.expect[1,:]), label = "100 trajectories", linestyle = :solid)
+lines!(ax, times, real(average_expect(sol_1)[1,:]), label = "1 trajectory", linestyle = :dashdot)
+lines!(ax, times, real(average_expect(sol_10)[1,:]), label = "10 trajectories", linestyle = :dash)
+lines!(ax, times, real(average_expect(sol_100)[1,:]), label = "100 trajectories", linestyle = :solid)
 
 axislegend(ax, position = :rt)
 

docs/src/users_guide/time_evolution/solution.md

@@ -89,7 +89,23 @@ Some other solvers can have other output.
 
 ## [Multiple trajectories solution](@id doc-TE:Multiple-trajectories-solution)
 
-This part is still under construction, please read the docstrings for the following functions first:
+The solutions are different for solvers which compute multiple trajectories, such as the [`TimeEvolutionMCSol`](@ref) (Monte Carlo) or the [`TimeEvolutionStochasticSol`](@ref) (stochastic methods). The expectation values and states for all trajectories at each time points will be stored in these solutions. The index-order of the elements in fields `states` and `expect` are:
 
-- [`TimeEvolutionMCSol`](@ref)
-- [`TimeEvolutionStochasticSol`](@ref)
+```julia
+sol.states[trajectory][time]
+sol.expect[e_op,trajectory,time]
+```
+
+We also provide the following functions for statistical analysis of multi-trajectory `sol`utions:
+
+| **Functions** | **Description** |
+|:------------|:----------------|
+| [`average_states(sol)`](@ref average_states) | Return the trajectory-averaged result states (as density [`Operator`](@ref)) at each time point. |
+| [`average_expect(sol)`](@ref average_expect) | Return the trajectory-averaged expectation values at each time point. |
+| [`std_expect(sol)`](@ref std_expect) | Return the trajectory-wise standard deviation of the expectation values at each time point. |
+
+Multi-trajectory solutions also keep the random number generator (`rng`) to allow recomputing the results:
+
+```julia
+rng = sol.rng # this can be specified as a keyword argument (`rng = rng`) to the solvers
+```

docs/src/users_guide/time_evolution/stochastic.md

@@ -111,11 +111,13 @@ sse_sol = ssesolve(
 measurement_avg = sum(sse_sol.measurement, dims=2) / size(sse_sol.measurement, 2)
 measurement_avg = dropdims(measurement_avg, dims=2)
 
+expect_avg = average_expect(sse_sol)
+
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time")
 lines!(ax, tlist[1:end-1], real(measurement_avg[1,:]), label = L"J_x", color = :red, linestyle = :solid)
-lines!(ax, tlist, real(sse_sol.expect[1,:]),  label = L"\langle x \rangle", color = :black, linestyle = :solid)
+lines!(ax, tlist, real(expect_avg[1,:]),  label = L"\langle x \rangle", color = :black, linestyle = :solid)
 
 axislegend(ax, position = :rt)
 
@@ -141,14 +143,16 @@ sme_sol = smesolve(
     store_measurement = Val(true),
 )
 
+expect_avg = average_expect(sme_sol)
+
 measurement_avg = sum(sme_sol.measurement, dims=2) / size(sme_sol.measurement, 2)
 measurement_avg = dropdims(measurement_avg, dims=2)
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time")
 lines!(ax, tlist[1:end-1], real(measurement_avg[1,:]), label = L"J_x", color = :red, linestyle = :solid)
-lines!(ax, tlist, real(sme_sol.expect[1,:]),  label = L"\langle x \rangle", color = :black, linestyle = :solid)
+lines!(ax, tlist, real(expect_avg[1,:]),  label = L"\langle x \rangle", color = :black, linestyle = :solid)
 
 axislegend(ax, position = :rt)
 

docs/src/users_guide/time_evolution/time_dependent.md

@@ -132,7 +132,7 @@ e_ops = [
 
 # solve dynamics 
 exp_me = mesolve(H_t, ψ0, tlist, c_ops; e_ops = e_ops, progress_bar = Val(false)).expect
-exp_mc = mcsolve(H_t, ψ0, tlist, c_ops; e_ops = e_ops, ntraj = 100, progress_bar = Val(false)).expect
+exp_mc = average_expect(mcsolve(H_t, ψ0, tlist, c_ops; e_ops = e_ops, ntraj = 100, progress_bar = Val(false)))
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))