-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathNumericalPendulumProblem.cs
More file actions
43 lines (37 loc) · 1.78 KB
/
NumericalPendulumProblem.cs
File metadata and controls
43 lines (37 loc) · 1.78 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
using minuit2.net;
using minuit2.net.CostFunctions;
namespace ExampleProblems.CustomProblems;
public class NumericalPendulumProblem : IConfiguredProblem
{
public NumericalPendulumProblem()
{
OptimumParameterValues = [1];
ParameterConfigurations = [ParameterConfiguration.Variable("pendulumLength", 2)];
var model = NumericalPendulumModelFor(startAngle: 1.5, startAngleVelocity: 0.0);
var xValues = Values.LinearlySpacedBetween(0, 3, 0.001);
var yValues = model(xValues, OptimumParameterValues.ToArray());
Cost = CostFunction.LeastSquares(xValues, yValues, ["pendulumLength"], model);
}
public ICostFunction Cost { get; }
public IReadOnlyCollection<double> OptimumParameterValues { get; }
public IReadOnlyCollection<ParameterConfiguration> ParameterConfigurations { get; }
private static Func<IReadOnlyList<double>, IReadOnlyList<double>, IReadOnlyList<double>> NumericalPendulumModelFor(
double startAngle,
double startAngleVelocity)
{
return (x, p) =>
{
const double acceleration = 9.81; // standard gravity (m/s2)
var deltaX = x.Zip(x.Skip(1), (lowerX, upperX) => upperX - lowerX).ToArray();
var angle = Enumerable.Repeat(startAngle, x.Count).ToArray();
var angleVelocity = Enumerable.Repeat(startAngleVelocity, x.Count).ToArray();
for (var i = 0; i < deltaX.Length; i++)
{
// ODE system (see e.g. https://de.wikipedia.org/wiki/Mathematisches_Pendel) solved using Euler method
angle[i + 1] = angle[i] + angleVelocity[i] * deltaX[i];
angleVelocity[i + 1] = angleVelocity[i] - acceleration / p[0] * Math.Sin(angle[i]) * deltaX[i];
}
return angle;
};
}
}