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
|
// Copyright (C) 2002 Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
// First added: 2002
// Last changed: 2006-08-22
//
// This demo solves the harmonic oscillator on
// the time interval (0, 4*pi) and computes the
// error for a set of methods and orders.
#include <dolfin.h>
using namespace dolfin;
class Harmonic : public ODE
{
public:
Harmonic() : ODE(2, 4.0 * DOLFIN_PI), e(0.0) {}
void u0(uBlasVector& u)
{
u(0) = 0.0;
u(1) = 1.0;
}
void f(const uBlasVector& u, real t, uBlasVector& y)
{
y(0) = u(1);
y(1) = - u(0);
}
bool update(const uBlasVector& u, real t, bool end)
{
if ( !end )
return true;
real e0 = u(0) - 0.0;
real e1 = u(1) - 1.0;
e = std::max(std::abs(e0), std::abs(e1));
return true;
}
real error()
{
return e;
}
private:
real e;
};
int main()
{
dolfin_set("ODE fixed time step", true);
dolfin_set("ODE discrete tolerance", 1e-14);
for (int q = 1; q <= 5; q++)
{
dolfin_set("output destination", "silent");
dolfin_set("ODE method", "cg");
dolfin_set("ODE order", q);
Harmonic ode;
ode.solve();
dolfin_set("output destination", "terminal");
message("cG(%d): e = %.3e", q, ode.error());
}
cout << endl;
for (int q = 0; q <= 5; q++)
{
dolfin_set("output destination", "silent");
dolfin_set("ODE method", "dg");
dolfin_set("ODE order", q);
Harmonic ode;
ode.solve();
dolfin_set("output destination", "terminal");
message("dG(%d): e = %.3e", q, ode.error());
}
return 0;
}
|