summaryrefslogtreecommitdiff
path: root/log-fourier_integrator.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'log-fourier_integrator.cpp')
-rw-r--r--log-fourier_integrator.cpp150
1 files changed, 83 insertions, 67 deletions
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 061cb08..7569245 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -4,15 +4,18 @@
#include <iostream>
int main(int argc, char* argv[]) {
+ /* Model parameters */
unsigned p = 2;
unsigned s = 2;
Real λ = 0.5;
Real τ₀ = 0;
+ /* Log-Fourier parameters */
unsigned log2n = 8;
Real Δτ = 0.1;
Real k = 0.1;
+ /* Iteration parameters */
Real ε = 1e-13;
Real γ = 1;
Real βₘₐₓ = 0.7;
@@ -64,90 +67,103 @@ int main(int argc, char* argv[]) {
μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
}
- std::vector<Real> C(N);
- std::vector<Real> R(N);
- std::vector<Complex> Ct(N);
- std::vector<Complex> Rt(N);
+ std::vector<Real> Cₜ(N);
+ std::vector<Real> Rₜ(N);
+ std::vector<Complex> Ĉₜ(N);
+ std::vector<Complex> Ȓₜ(N);
- // start from the exact solution for β = 0
+ /* Start from the exact solution for β = 0 */
for (unsigned n = 0; n < N; n++) {
if (τ₀ > 0) {
- C[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
+ Cₜ[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
} else {
- C[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
+ Cₜ[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
}
- R[n] = exp(-μ * fft.t(n));
+ Rₜ[n] = exp(-μ * fft.t(n));
- Ct[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
- Rt[n] = 1.0 / (μ + 1i * fft.ν(n));
+ Ĉₜ[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
+ Ȓₜ[n] = 1.0 / (μ + 1i * fft.ν(n));
}
Real β = 0;
while (β < βₘₐₓ) {
- Real ΔC = 100;
- while (ΔC > ε) {
- std::vector<Real> RddfC(N);
- std::vector<Real> dfC(N);
- for (unsigned n = 0; n < N; n++) {
- RddfC[n] = R[n] * ddf(λ, p, s, C[n]);
- dfC[n] = df(λ, p, s, C[n]);
+ Real ΔC = 100;
+ Real ΔC₀ = 100;
+ unsigned it = 0;
+ while (ΔC > ε) {
+ std::vector<Real> RddfC(N);
+ for (unsigned n = 0; n < N; n++) {
+ RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]);
+ }
+ std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
+
+ std::vector<Complex> Ȓₜ₊₁(N);
+ std::vector<Complex> Ĉₜ₊₁(N);
+
+ for (unsigned n = 0; n < N; n++) {
+ Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n));
+ Ĉₜ₊₁[n] = - 2 * Γ₀ * Ȓₜ₊₁[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n);
+ }
+
+ std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
+ std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
+
+ ΔC = 0;
+ for (unsigned i = 0; i < N; i++) {
+ ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]);
+ ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]);
+ }
+ ΔC = sqrt(ΔC) / (2*N);
+
+ for (unsigned i = 0; i < N; i++) {
+ Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
+ Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]);
+ Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]);
+ Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
+ }
+
+ μ = γ * μ * Cₜ[0] + (1-γ) * μ;
+
+ if (ΔC < ΔC₀) {
+ ΔC₀ = ΔC;
+ it = 0;
+ γ = std::min(1.001 * γ, 1.0);
+ } else {
+ it++;
+ }
+
+ if (it > 100) {
+ γ = std::max(0.5 * γ, 1e-3);
+ it = 0;
+ ΔC₀ = ΔC;
+ }
+
+ std::cerr << β << " " << μ << " " << ΔC << " " << γ;
+ std::cerr << "\r";
}
- std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
- std::vector<Complex> dfCt = fft.fourier(dfC, true);
- std::vector<Complex> Rtnew(N);
- std::vector<Complex> Ctnew(N);
-
- for (unsigned n = 0; n < N; n++) {
- Rtnew[n] = (1.0 + pow(β, 2) * RddfCt[n] * Rt[n]) / (μ + 1i * fft.ν(n));
- Ctnew[n] = (2 * Γ₀ * std::conj(Rtnew[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ct[n] + dfCt[n] * std::conj(Rtnew[n]))) / (μ + 1i * fft.ν(n));
-// Ctnew[n] = - 2 * Γ₀ * Rtnew[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n);
- }
-
- std::vector<Real> Cnew = fft.inverse(Ctnew);
- std::vector<Real> Rnew = fft.inverse(Rtnew);
-
- ΔC = 0;
- for (unsigned i = 0; i < N; i++) {
- ΔC += std::norm(Ct[i] - Ctnew[i]);
- ΔC += std::norm(Rt[i] - Rtnew[i]);
- }
- ΔC = sqrt(ΔC) / (2*N);
-
- for (unsigned i = 0; i < N; i++) {
- C[i] += γ * (Cnew[i] - C[i]);
- R[i] += γ * (Rnew[i] - R[i]);
- Ct[i] += γ * (Ctnew[i] - Ct[i]);
- Rt[i] += γ * (Rtnew[i] - Rt[i]);
+ /* Integrate the energy using Simpson's rule */
+ Real E = 0;
+ for (unsigned n = 0; n < N/2-1; n++) {
+ Real h₂ₙ = fft.t(2*n+1) - fft.t(2*n);
+ Real h₂ₙ₊₁ = fft.t(2*n+2) - fft.t(2*n+1);
+ Real f₂ₙ = Rₜ[2*n] * df(λ, p, s, Cₜ[2*n]);
+ Real f₂ₙ₊₁ = Rₜ[2*n+1] * df(λ, p, s, Cₜ[2*n+1]);
+ Real f₂ₙ₊₂ = Rₜ[2*n+2] * df(λ, p, s, Cₜ[2*n+2]);
+ E += (h₂ₙ + h₂ₙ₊₁) / 6 * (
+ (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ
+ + pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
+ + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂
+ );
}
+ E *= β;
- μ *= C[0];
-
-// std::cerr << ΔC << std::endl;
- }
-
- /* Integrate the energy using Simpson's rule */
- Real E = 0;
- for (unsigned i = 0; i < N/2-1; i++) {
- Real h₂ᵢ = fft.t(2*i+1) - fft.t(2*i);
- Real h₂ᵢ₊₁ = fft.t(2*i+2) - fft.t(2*i+1);
- Real f₂ᵢ = R[2*i] * df(λ, p, s, C[2*i]);
- Real f₂ᵢ₊₁ = R[2*i+1] * df(λ, p, s, C[2*i+1]);
- Real f₂ᵢ₊₂ = R[2*i+2] * df(λ, p, s, C[2*i+2]);
- E += (h₂ᵢ + h₂ᵢ₊₁) / 6 * (
- (2 - h₂ᵢ₊₁ / h₂ᵢ) * f₂ᵢ
- + pow(h₂ᵢ + h₂ᵢ₊₁, 2) / (h₂ᵢ * h₂ᵢ₊₁) * f₂ᵢ₊₁
- + (2 - h₂ᵢ / h₂ᵢ₊₁) * f₂ᵢ₊₂
- );
- }
- E *= β;
-
- std::cerr << β << " " << μ << " " << Ct[0].real() << " " << E << std::endl;
- β += Δβ;
+ std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl;
+ β += Δβ;
}
for (unsigned i = 0; i < N; i++) {
- std::cout << fft.t(i) << " " << C[i] << std::endl;
+ std::cout << fft.t(i) << " " << Cₜ[i] << std::endl;
}
return 0;