1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
#include "fourier.hpp"
#include <getopt.h>
#include <iostream>
Real energy(const std::vector<Real>& C, std::vector<Real>& R, Real λ, unsigned p, unsigned s, Real Δτ) {
Real I = 0;
for (unsigned σ = 0; σ < C.size(); σ++) {
I += Δτ * df(λ, p, s, C[σ]) * R[σ];
}
return I;
}
int main(int argc, char* argv[]) {
unsigned p = 3;
unsigned s = 4;
Real λ = 0.5;
Real Δτ = 1e-3;
Real τₘₐₓ = 1e3;
Real τ₀ = 0;
Real β = 0.5;
unsigned iterations = 10;
int opt;
while ((opt = getopt(argc, argv, "d:T:t:b:I:")) != -1) {
switch (opt) {
case 'd':
Δτ = atof(optarg);
break;
case 'T':
τₘₐₓ = atof(optarg);
break;
case 't':
τ₀ = atof(optarg);
break;
case 'b':
β = atof(optarg);
break;
case 'I':
iterations = (unsigned)atof(optarg);
break;
default:
exit(1);
}
}
Real Γ₀ = 1;
Real μ = 1;
if (τ₀ > 0) {
μ = (sqrt(1+4*Γ₀*τ₀) - 1) / (2 * τ₀);
}
Real τ = 0;
unsigned N = τₘₐₓ / Δτ + 1;
std::vector<Real> C(N);
std::vector<Real> R(N);
std::vector<Real> Γ(N);
for (unsigned i = 0; i < N; i++) {
Real τ = i * Δτ;
if (τ₀ > 0) {
C[i] = (Γ₀ / μ) * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (1 - pow(μ * τ₀, 2));
Γ[i] = (Γ₀ / τ₀) * exp(-τ / τ₀);
} else {
C[i] = (Γ₀ / μ) * exp(-μ * τ);
}
R[i] = exp(-μ * τ);
}
for (unsigned it = 0; it < iterations; it++) {
/* First step: integrate R from C */
std::vector<Real> R₊(N);
R₊[0] = 1;
for (unsigned i = 1; i < N; i++) {
Real I = 0;
for (unsigned j = 0; j <= i; j++) {
I += R[i - j] * ddf(λ, p, s, C[i - j]) * R[j] * Δτ;
}
Real dR = -μ * R[i] + pow(β, 2) * I;
R₊[i] = R₊[i - 1] + dR * Δτ;
}
/* Second step: integrate C from R */
}
τ = 0;
for (Real Ci : C) {
std::cout << τ << " " << Ci << std::endl;
τ += Δτ;
}
std::cerr << - 2 * β / Γ₀ * energy(C, R, λ, p, s, Δτ) << std::endl;
return 0;
}
|
