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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
#include <getopt.h>
#include <vector>
#include <array>
#include <cmath>
#include <iostream>
using Real = double;
class Point : public std::array<Real, 2> {
public:
Real τ() const {
return front();
}
Real C() const {
return back();
}
};
unsigned p = 2;
Real f(Real q) {
return 0.5 * pow(q, p);
}
Real df(Real q) {
return 0.5 * p * pow(q, p - 1);
}
Real ddf(Real q) {
return 0.5 * p * (p - 1) * pow(q, p - 2);
}
Real integrate(const std::vector<Point>& Cₜ, double τ₀) {
Real I = 0;
#pragma omp parallel for reduction(+:I)
for (unsigned i = 0; i < Cₜ.size() - 1; i++) {
unsigned j = Cₜ.size() - 1 - i;
Real Δτ = Cₜ[i + 1].τ() - Cₜ[i].τ();
Real C = (Cₜ[j].C() + Cₜ[j - 1].C()) / 2;
Real dC = (Cₜ[i + 1].C() - Cₜ[i].C()) / Δτ;
Real dddC = 0;
if (i > 4 && i < Cₜ.size() && Cₜ.size() > 4) {
dddC += (Cₜ[j].C() - 3 * Cₜ[j+1].C() + 3 * Cₜ[j+2].C() - Cₜ[j+3].C()) / pow(Δτ, 3);
dddC += (Cₜ[j+1].C() - 3 * Cₜ[j+2].C() + 3 * Cₜ[j+3].C() - Cₜ[j+4].C()) / pow(Δτ, 3);
dddC /= 2;
}
I += Δτ * df(C) * (dC - pow(τ₀, 2) * dddC);
}
return I;
}
Real energy(const std::vector<Point>& Ct, Real τ₀) {
Real I = 0;
for (unsigned i = 0; i < Ct.size() - 1; i++) {
Real Δτ = Ct[i + 1].τ() - Ct[i].τ();
Real C = (Ct[i].C() + Ct[i + 1].C()) / 2;
Real dC = (Ct[i + 1].C() - Ct[i].C()) / Δτ;
Real dddC = 0;
if (i > 1 && i < Ct.size() - 2 && Ct.size() > 3) {
dddC = (Ct[i+1].C() - 3 * Ct[i].C() + 3 * Ct[i-1].C() - Ct[i-2].C()) / pow(Δτ, 3);
}
I += Δτ * df(C) * (dC - pow(τ₀, 2) * dddC);
}
return I;
}
int main(int argc, char* argv[]) {
Real Δτ = 1e-3;
Real τₘₐₓ = 1e3;
Real τ₀ = 0;
Real y = 0.5;
int opt;
while ((opt = getopt(argc, argv, "d:T:t:y:")) != -1) {
switch (opt) {
case 'd':
Δτ = atof(optarg);
break;
case 'T':
τₘₐₓ = atof(optarg);
break;
case 't':
τ₀ = atof(optarg);
break;
case 'y':
y = atof(optarg);
break;
default:
exit(1);
}
}
Real z = (sqrt(1 + 2 * τ₀) - 1) / (2 * τ₀);
std::vector<Point> Cₜ;
Cₜ.reserve(τₘₐₓ / Δτ + 1);
Cₜ.push_back({0, 1});
while (Cₜ.back().τ() < τₘₐₓ) {
Real dC = -z * Cₜ.back().C() - 2 * pow(y, 2) * integrate(Cₜ, τ₀);
Cₜ.push_back({Cₜ.back().τ() + Δτ, Cₜ.back().C() + Δτ * dC});
std::cout << Cₜ.back().τ() << " " << Cₜ.back().C() << std::endl;
}
std::cerr << - 2 * y * energy(Cₜ, τ₀) << std::endl;
return 0;
}
|
