jacob/BD_Integratoren
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Konv added

Browse Source
This commit is contained in:
Jacob Holder 2021-12-16 16:20:41 +01:00
parent 6b69a2f41b
commit 82e81eb91e
Signed by: jacob
GPG Key ID: 2194FC747048A7FD
6 changed files with 145 additions and 95 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
Auswertung/harmonic.py
View File

@ -1,38 +1,38 @@
import glob
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sb
import numpy as np import numpy as np
import pandas as pd
def main(): def main():
time = np.arange(0.0, 4.0, 0.01) time = np.arange(0.0, 4.0, 0.01)
xstart = 1.0 xstart = 1.0
ystart = 1.0 ystart = 1.0
x_mean = xstart * np.exp(-time) x_mean = xstart * np.exp(-time)
y_mean = ystart * np.exp(-time) y_mean = ystart * np.exp(-time)
x_sq_mean = xstart*xstart * np.exp(-2*time) + (1 - np.exp(-2*time)) x_sq_mean = xstart * xstart * np.exp(-2 * time) + (1 - np.exp(-2 * time))
y_sq_mean = ystart*ystart * np.exp(-2 * time) + (1 - np.exp(-2 * time)) y_sq_mean = ystart * ystart * np.exp(-2 * time) + (1 - np.exp(-2 * time))
x = pd.read_csv("./data/out/harmonic_force_euler_L0_x.dat") x = pd.read_csv("./data/out/harmonic_force_euler_L0_x.dat")
x_sq = pd.read_csv("./data/out/harmonic_force_euler_L0_x_squared.dat") x_sq = pd.read_csv("./data/out/harmonic_force_euler_L0_x_squared.dat")
msd = pd.read_csv("./data/out/harmonic_force_euler_L0_msd.dat") msd = pd.read_csv("./data/out/harmonic_force_euler_L0_msd.dat")
msd_mean = (1-np.exp(-2*time)) msd_mean = (1 - np.exp(-2 * time))
msd_2_mean = 2 * (1 - 2 * np.exp(- time) + np.exp(-2 * time)) + 2 * (1 - np.exp(-2 * time))
plt.plot(time, x_mean, label="a <x>") plt.plot(time, x_mean, label="a <x>")
plt.plot(x["time"], x["val"], label="<x>") plt.plot(x["time"], x["val"], label="<x>")
plt.plot(time, x_sq_mean, label="a <x²>") plt.plot(time, x_sq_mean, label="a <x²>")
plt.plot(x_sq["time"], x_sq["val"], label="<x²>") plt.plot(x_sq["time"], x_sq["val"], label="<x²>")
plt.plot(time, 2*msd_mean,label="a msd") plt.plot(time, 2 * msd_mean, label="a msd")
plt.plot(time, msd_2_mean, label="a msd")
plt.plot(msd["time"], msd["val"], label="msd") plt.plot(msd["time"], msd["val"], label="msd")
#plt.plot(x["time"], x_sq["val"] - np.power(x["val"], 2.0)) # plt.plot(x["time"], x_sq["val"] - np.power(x["val"], 2.0))
plt.plot(x["time"], 2.0*(x_sq["val"] - np.power(x["val"], 2.0)), label="<x²>-<x>²") plt.plot(x["time"], 2 * (x_sq["val"] - np.power(x["val"], 2.0)), label="<x²>-<x>²")
plt.legend(loc = 1) plt.legend(loc=1)
plt.show() plt.show()

8
Auswertung/main.py
View File

@ -21,11 +21,15 @@ def make_figure(length_index, length, forces, integrators, computes):
#ax.fill_between(file["time"], y1=file["val"] - file["std"], y2=file["val"] + file["std"], alpha=.5) #ax.fill_between(file["time"], y1=file["val"] - file["std"], y2=file["val"] + file["std"], alpha=.5)
plt.ylabel(compute) plt.ylabel(compute)
legend.append(integrator) legend.append(integrator)
if index_integrator == 0 and not np.any(file["target"].isnull()): if index_integrator == 0 and not np.any(file["target"].isnull()):
plt.plot(file["time"], file["target"]) plt.plot(file["time"], file["target"])
legend.append("analytisch") legend.append("analytisch")
plt.legend(title='iterators', labels=legend) plt.legend(title='iterators', labels=legend)
for ax1 in axs:
for ax in ax1:
ax.axvline(1.0)
pad = 5 # in points pad = 5 # in points
for ax, col in zip(axs[0], computes): for ax, col in zip(axs[0], computes):
@ -40,13 +44,13 @@ def make_figure(length_index, length, forces, integrators, computes):
fig.tight_layout() fig.tight_layout()
fig.savefig(f"fig{length_index}.png") fig.savefig(f"fig{length_index}.png")
plt.show()
def main(): def main():
forces = ["zero", "const", "harmonic"] forces = ["zero", "const", "harmonic"]
integrators = ["euler", "heun", "exact", "bdas", "mbd"] integrators = ["euler", "heun", "exact", "bdas", "mbd"]
computes = ["msd", "oaf", "empxx", "empyy", "x", "x_squared"] computes = ["msd", "oaf", "empxx", "empyy", "x", "x_squared"]
for l_i, l in enumerate(["sphere", "L = 1.5", "L = 2.0"]): for l_i, l in enumerate(["sphere"]): #["sphere", "L = 1.5", "L = 2.0"]
make_figure(length=l, length_index=l_i, forces=forces, computes=computes, integrators=integrators) make_figure(length=l, length_index=l_i, forces=forces, computes=computes, integrators=integrators)

4
C++/src/Calculation.cpp
View File

@ -4,15 +4,19 @@
#include "Calculation.h" #include "Calculation.h"
#include <fstream>
#include <utility> #include <utility>
void Calculation::run(size_t steps) { void Calculation::run(size_t steps) {
std::ofstream file ("dump.xyz");
for (size_t step = 0; step < steps; ++step) { for (size_t step = 0; step < steps; ++step) {
m_integrator(rod, sim); m_integrator(rod, sim);
for (auto &comp : computes) { for (auto &comp : computes) {
comp.eval(rod); comp.eval(rod);
} }
file<<rod.getPos()[0]<<" "<<rod.getPos()[0]<<" 0\n";
} }
file.close();
} }
const Rod2d &Calculation::getRod() const { return rod; } const Rod2d &Calculation::getRod() const { return rod; }

73
C++/src/Compute.cpp
View File

@ -13,7 +13,7 @@
void Compute::evalMSD(const Rod2d &rod2D) { void Compute::evalMSD(const Rod2d &rod2D) {
const auto &new_pos = rod2D.getPos(); const auto &new_pos = rod2D.getPos();
auto old_pos = start_rod.getPos(); auto old_pos = start_rod.getPos();
auto diff = (new_pos - old_pos); auto diff = new_pos - old_pos;
auto msd = diff.dot(diff); auto msd = diff.dot(diff);
agg.feed(msd); agg.feed(msd);
} }
@ -26,7 +26,7 @@ void Compute::evalX(const Rod2d &rod2D) {
void Compute::evalX_squared(const Rod2d &rod2D) { void Compute::evalX_squared(const Rod2d &rod2D) {
const auto &new_pos = rod2D.getPos(); const auto &new_pos = rod2D.getPos();
auto msd = pow(new_pos[0],2.0); auto msd = pow(new_pos[0], 2.0);
agg.feed(msd); agg.feed(msd);
} }
@ -57,37 +57,34 @@ void Compute::eval_empYY(const Rod2d &rod2D) {
void Compute::eval(const Rod2d &rod2D) { void Compute::eval(const Rod2d &rod2D) {
time_step++; time_step++;
if (time_step % every == 0) { if (time_step % every == 0) {
switch (type) { if (resetting or time_step == every) {
case Type::msd: switch (type) {
evalMSD(rod2D); case Type::msd:
break; evalMSD(rod2D);
case Type::oaf: break;
evalOAF(rod2D); case Type::oaf:
break; evalOAF(rod2D);
case Type::empxx: break;
eval_empXX(rod2D); case Type::empxx:
break; eval_empXX(rod2D);
case Type::empyy: break;
eval_empYY(rod2D); case Type::empyy:
break; eval_empYY(rod2D);
case Type::x: break;
if (time_step == every) { case Type::x:
evalX(rod2D); evalX(rod2D);
} break;
break; case Type::x_squared:
case Type::x_squared:
if (time_step == every) {
evalX_squared(rod2D); evalX_squared(rod2D);
} break;
break; }
if (resetting) {
start_rod = rod2D;
}
return;
} }
if (resetting) {
start_rod = rod2D;
}
return;
} }
} }
Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
Simulation &sim) Simulation &sim)
: start_rod(std::move(rod)), : start_rod(std::move(rod)),
@ -97,7 +94,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
time(sim.getMDeltaT() * static_cast<double>(every)) { time(sim.getMDeltaT() * static_cast<double>(every)) {
switch (type) { switch (type) {
case Type::msd: { case Type::msd: {
resetting = true; resetting = false;
type_str = "msd"; type_str = "msd";
switch (force) { switch (force) {
case Force::zero_F: case Force::zero_F:
@ -107,20 +104,20 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
target = 0.0; target = 0.0;
break; break;
case Force::harmonic_F: case Force::harmonic_F:
target = 2 * (1.0 - std::exp(-2 * time)); target = 2 * (1 - 2 * std::exp(- time) + std::exp(-2 * time)) + 2 * (1 - std::exp(-2 * time));
break; break;
} }
break; break;
} }
case Type::oaf: { case Type::oaf: {
resetting = true; resetting = false;
type_str = "oaf"; type_str = "oaf";
target = std::exp(-rod.getDRot() * time); target = std::exp(-rod.getDRot() * time);
break; break;
} }
case Type::empxx: { case Type::empxx: {
resetting = true; resetting = false;
type_str = "empxx"; type_str = "empxx";
const double Dmean = 0.5 * (rod.getDiff().trace()); const double Dmean = 0.5 * (rod.getDiff().trace());
const double u = 4.0; const double u = 4.0;
@ -138,14 +135,15 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
} }
} break; } break;
case Type::empyy: { case Type::empyy: {
resetting = true; resetting = false;
type_str = "empyy"; type_str = "empyy";
const double Dmean = 0.5 * (rod.getDiff().trace()); const double Dmean = 0.5 * (rod.getDiff().trace());
const double u = 4.0; const double u = 4.0;
const double deltaD = rod.getDiff()(1, 1) - rod.getDiff()(0, 0); const double deltaD = rod.getDiff()(1, 1) - rod.getDiff()(0, 0);
switch (force) { switch (force) {
case Force::zero_F: case Force::zero_F:
target = Dmean + deltaD / 2 * (1 - exp(-u * rod.getDRot() * time)) / target = Dmean + deltaD / 2 *
(1 - exp(-u * rod.getDRot() * time)) /
u / rod.getDRot() / time; u / rod.getDRot() / time;
break; break;
case Force::const_F: case Force::const_F:
@ -169,8 +167,6 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
break; break;
} }
} break; } break;
case Type::x_squared: { case Type::x_squared: {
resetting = false; resetting = false;
@ -184,12 +180,15 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
target = 0.0; target = 0.0;
break; break;
case Force::harmonic_F: case Force::harmonic_F:
target = pow(rod.getPos()[0],2) * exp(-2 * time) + target = pow(rod.getPos()[0], 2) * exp(-2 * time) +
1 * (1 - exp(-2 * time)); 1 * (1 - exp(-2 * time));
break; break;
} }
} break; } break;
} }
resetting = false;
} }
const LiveAgg &Compute::getAgg() const { return agg; } const LiveAgg &Compute::getAgg() const { return agg; }

5
C++/src/Rod2d.cpp
View File

@ -20,8 +20,8 @@ Rod2d::Rod2d(double L) : m_pos({1, 1}), m_e({1, 0}) {
const double p = L; const double p = L;
auto D_para = auto D_para =
D0 / (M_Pl * 2.0) * (log(p) - 0.1401 + 1.034 / p - 0.228 / (p * p)); D0 / (M_Pl * 2.0) * (log(p) - 0.1401 + 1.034 / p - 0.228 / (p * p));
auto D_ortho = auto D_ortho = D0 / (M_Pl * 4.0) *
D0 / (M_Pl * 4.0) * (log(p) + 0.8369 + 0.5551 / p - 0.06066 / (p * p)); (log(p) + 0.8369 + 0.5551 / p - 0.06066 / (p * p));
m_D_rot = 3 * D0 / (M_Pl * L * L) * m_D_rot = 3 * D0 / (M_Pl * L * L) *
(log(p) - 0.3512 + 0.7804 / p - 0.09801 / (p * p)); (log(p) - 0.3512 + 0.7804 / p - 0.09801 / (p * p));
m_Diff << D_para, 0, 0, D_ortho; m_Diff << D_para, 0, 0, D_ortho;
@ -46,4 +46,3 @@ double Rod2d::getDRot_Sqrt() const { return m_D_rot_sqrt; }
const Eigen::Vector2d &Rod2d::getE() const { return m_e; } const Eigen::Vector2d &Rod2d::getE() const { return m_e; }
void Rod2d::setE(const Eigen::Vector2d &mE) { m_e = mE; } void Rod2d::setE(const Eigen::Vector2d &mE) { m_e = mE; }

126
C++/src/main.cpp
View File

@ -6,18 +6,21 @@
#include <iostream> #include <iostream>
#include <filesystem> #include <filesystem>
#include <utility>
#include "Calculation.h" #include "Calculation.h"
#include "Integratoren2d.h" #include "Integratoren2d.h"
#include "force_lambdas.h" #include "force_lambdas.h"
constexpr size_t SEED = 1234; #include <fmt/core.h>
constexpr size_t SEED = 1364;
constexpr double stepSize = 0.01; constexpr double stepSize = 0.01;
constexpr size_t n_computes = 100; constexpr size_t n_computes = 10;
constexpr size_t num_integratoren = 5; constexpr size_t num_integratoren = 5;
constexpr size_t num_forces = 3; constexpr size_t num_forces = 3;
constexpr size_t num_length = 3; constexpr size_t num_length = 1;
constexpr size_t delta_compute = 1; constexpr size_t delta_compute = 10;
constexpr size_t numStep = 10000; constexpr size_t numStep = 50000;
inline const std::string header = {"time,val,target,std\n"}; inline const std::string header = {"time,val,target,std,numPoints\n"};
void run(size_t integrator_index, size_t force_index, size_t length_index) { void run(size_t integrator_index, size_t force_index, size_t length_index) {
using type = std::pair<Calculation::inte_force_type, std::string>; using type = std::pair<Calculation::inte_force_type, std::string>;
std::vector<type> vec({type(Integratoren2d::Set1_Euler_f, "euler"), std::vector<type> vec({type(Integratoren2d::Set1_Euler_f, "euler"),
@ -41,7 +44,8 @@ void run(size_t integrator_index, size_t force_index, size_t length_index) {
std::string folder = std::string folder =
"out/" + force_name + name + "_L" + std::to_string(length_index) + "_"; "out/" + force_name + name + "_L" + std::to_string(length_index) + "_";
fs::create_directories("out"); fs::create_directories("out");
double length = std::vector<double>({1.0, 1.5, 2.0})[length_index]; // double length = std::vector<double>({1.0, 1.5, 2.0})[length_index];
double length = 1.0;
{ {
std::vector<std::pair<Compute::Type, size_t>> computes; std::vector<std::pair<Compute::Type, size_t>> computes;
for (size_t i = 1; i < n_computes; ++i) { for (size_t i = 1; i < n_computes; ++i) {
@ -49,6 +53,8 @@ void run(size_t integrator_index, size_t force_index, size_t length_index) {
computes.emplace_back(Compute::Type::oaf, i * delta_compute); computes.emplace_back(Compute::Type::oaf, i * delta_compute);
computes.emplace_back(Compute::Type::empxx, i * delta_compute); computes.emplace_back(Compute::Type::empxx, i * delta_compute);
computes.emplace_back(Compute::Type::empyy, i * delta_compute); computes.emplace_back(Compute::Type::empyy, i * delta_compute);
computes.emplace_back(Compute::Type::x, i * delta_compute);
computes.emplace_back(Compute::Type::x_squared, i * delta_compute);
} }
Calculation euler(integrator, computes, stepSize, SEED, force, Calculation euler(integrator, computes, stepSize, SEED, force,
@ -56,12 +62,18 @@ void run(size_t integrator_index, size_t force_index, size_t length_index) {
for (size_t i = 0; i < numStep; ++i) { for (size_t i = 0; i < numStep; ++i) {
euler.run(n_computes * delta_compute); euler.run(n_computes * delta_compute);
euler.reset();
} }
std::ofstream msdFile(folder + "msd.dat"); std::ofstream msdFile(folder + "msd.dat");
std::ofstream oafFile(folder + "oaf.dat"); std::ofstream oafFile(folder + "oaf.dat");
std::ofstream empxxFile(folder + "empxx.dat"); std::ofstream empxxFile(folder + "empxx.dat");
std::ofstream empyyFile(folder + "empyy.dat"); std::ofstream empyyFile(folder + "empyy.dat");
std::ofstream xFile(folder + "x.dat");
std::ofstream xsqFile(folder + "x_squared.dat");
xFile << header;
xsqFile << header;
msdFile << header; msdFile << header;
oafFile << header; oafFile << header;
empxxFile << header; empxxFile << header;
@ -71,70 +83,102 @@ void run(size_t integrator_index, size_t force_index, size_t length_index) {
if (com.getType() == Compute::Type::msd) { if (com.getType() == Compute::Type::msd) {
msdFile << com.getTime() << ", " << com.getAgg().getMean() msdFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", " << ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl; << com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
} }
if (com.getType() == Compute::Type::oaf) { if (com.getType() == Compute::Type::oaf) {
oafFile << com.getTime() << ", " << com.getAgg().getMean() oafFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", " << ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl; << com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
} }
if (com.getType() == Compute::Type::empxx) { if (com.getType() == Compute::Type::empxx) {
empxxFile << com.getTime() << ", " << com.getAgg().getMean() empxxFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", " << ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl; << com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
} }
if (com.getType() == Compute::Type::empyy) { if (com.getType() == Compute::Type::empyy) {
empyyFile << com.getTime() << ", " << com.getAgg().getMean() empyyFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", " << ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl; << com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
} }
if (com.getType() == Compute::Type::msd) { if (com.getType() == Compute::Type::msd) {
msdFile << com.getTime() << ", " << com.getAgg().getMean() msdFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", " << ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl; << com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
} }
}
}
{
std::vector<std::pair<Compute::Type, size_t>> repeating_computes;
for (size_t i = 1; i < n_computes; ++i) {
repeating_computes.emplace_back(Compute::Type::x,
i * delta_compute);
repeating_computes.emplace_back(Compute::Type::x_squared,
i * delta_compute);
}
Calculation calc_repeat(integrator, repeating_computes, stepSize, SEED,
force, force_type, zero_Torque, 1.0);
for (size_t i = 0; i < numStep; ++i) {
calc_repeat.run(n_computes * delta_compute);
calc_repeat.reset();
}
std::ofstream xFile(folder + "x.dat");
std::ofstream xsqFile(folder + "x_squared.dat");
xFile << header;
xsqFile << header;
for (const auto &com : calc_repeat.getComputes()) {
if (com.getType() == Compute::Type::x) { if (com.getType() == Compute::Type::x) {
xFile << com.getTime() << ", " << com.getAgg().getMean() << ", " xFile << com.getTime() << ", " << com.getAgg().getMean() << ", "
<< com.getTarget() << ", " << com.getAgg().getSEM() << com.getTarget() << ", " << com.getAgg().getSEM()
<< std::endl; << ", " << com.getAgg().getNumPoints() << std::endl;
} }
if (com.getType() == Compute::Type::x_squared) { if (com.getType() == Compute::Type::x_squared) {
xsqFile << com.getTime() << ", " << com.getAgg().getMean() xsqFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", " << ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl; << com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
} }
} }
} }
std::cout << "Finished run " << folder << std::endl; std::cout << "Finished run " << folder << std::endl;
} }
int main() {
run(0, 2, 0); void konv(int delta, std::ofstream &file,
Calculation::inte_force_type integrator, size_t seed) {
std::vector<std::pair<Compute::Type, size_t>> computes;
double step_Size = 1.0 / delta;
computes.emplace_back(Compute::Type::msd, delta);
Calculation euler(integrator, computes, step_Size, seed, harmonic_Force,
Compute::Force::harmonic_F, zero_Torque, 1.0);
for (int j = 0; j < 1000000; ++j) {
for (int i = 0; i < 1000; ++i) {
euler.run(static_cast<size_t>(delta));
euler.reset();
}
std::cout << euler.getComputes()[0].getAgg().getSEM() << " "
<< euler.getComputes()[0].getDifference() << std::endl;
if (euler.getComputes()[0].getAgg().getSEM() * 10 <
euler.getComputes()[0].getDifference()) {
file << delta << " " << euler.getComputes()[0].getAgg().getMean()
<< " " << euler.getComputes()[0].getAgg().getSEM() << " "
<< euler.getComputes()[0].getTarget() << " "
<< euler.getComputes()[0].getDifference() << std::endl;
return;
}
}
file << "timeout: " << delta << " "
<< euler.getComputes()[0].getAgg().getMean() << " "
<< euler.getComputes()[0].getAgg().getSEM() << " "
<< euler.getComputes()[0].getTarget() << " "
<< euler.getComputes()[0].getDifference() << std::endl;
}
Calculation::inte_force_type integrator(int index) {
if (index % 5 == 0) return Integratoren2d::Set1_Euler_f;
if (index % 5 == 1) return Integratoren2d::Set2_Heun_f;
if (index % 5 == 2) return Integratoren2d::Set3_Exact_f;
if (index % 5 == 3) return Integratoren2d::Set4_BDAS_f;
if (index % 5 == 4) return Integratoren2d::Set5_MBD_f;
return Integratoren2d::Set1_Euler_f;
}
void konv_runner(int index) {
std::ofstream file;
file.open(fmt::format("out/{}_konv.out", index));
for (int i = 1; i < 256; i *= 2) {
konv(i, file, integrator(index), static_cast<size_t>(index / 5));
}
file.close();
}
int main(int argc, char *argv[]) {
int index = 0;
if (argc >1) index = std::stoi(argv[1]);
konv_runner(index);
return EXIT_SUCCESS; return EXIT_SUCCESS;
#pragma omp parallel for default(none) #pragma omp parallel for default(none)
for (size_t j = 0; j < num_forces * num_length * num_integratoren; ++j) { for (size_t j = 0; j < num_forces * num_length * num_integratoren; ++j) {