-
Notifications
You must be signed in to change notification settings - Fork 25
Expand file tree
/
Copy patheuler_test.cpp
More file actions
79 lines (64 loc) · 2.32 KB
/
Copy patheuler_test.cpp
File metadata and controls
79 lines (64 loc) · 2.32 KB
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
/*
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Daniel Abele, Martin J. Kuehn
*
* Contact: Martin J. Kuehn <[email protected]>
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include "memilio/math/euler.h"
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <cmath>
#include <numbers>
void init_vectors(std::vector<Eigen::VectorXd>& y, std::vector<Eigen::VectorXd>& sol, size_t n)
{
y = std::vector<Eigen::VectorXd>(n, Eigen::VectorXd::Constant(1, 0));
sol = std::vector<Eigen::VectorXd>(n, Eigen::VectorXd::Constant(1, 0));
}
// Test for y'(t) = cos(t)
void integration_test(std::vector<Eigen::VectorXd>& y, std::vector<Eigen::VectorXd>& sol, size_t n, double dt,
double& err)
{
sol[0][0] = std::sin(0);
sol[n - 1][0] = std::sin((n - 1) * dt);
auto f = [](auto&& /*y*/, auto&& t, auto&& dydt) {
dydt[0] = std::cos(t);
};
double t = 0.;
for (size_t i = 0; i < n - 1; i++) {
sol[i + 1][0] = std::sin((i + 1) * dt);
mio::EulerIntegratorCore<double>().step(f, y[i], t, dt, y[i + 1]);
printf("\n %.8f\t %.8f", y[i + 1][0], sol[i + 1][0]);
// printf("\n approx: %.4e, sol: %.4e, error %.4e", y[i+1][0], sol[i+1][0], err);
err += std::pow(std::abs(y[i + 1][0] - sol[i + 1][0]), 2.0);
}
}
int main()
{
std::vector<Eigen::VectorXd> y;
std::vector<Eigen::VectorXd> sol;
const double pi = std::numbers::pi_v<double>;
size_t n = 10;
double t0 = 0;
double tmax = 2 * pi;
double dt = (tmax - t0) / n;
double err = 0;
printf("\n .%.8f. \n", dt);
init_vectors(y, sol, n);
integration_test(y, sol, n, dt, err);
err = std::sqrt(err) / n;
printf("\nFor n=%d the error is %.4e\n", (int)n, err);
}