Skip to content

Commit 0475b42

Browse files
dragunovdenisdragunovdenis
committed
An autonomous problem with known (u(x)=exp(sin(x))) solution + the corresponding unit test + small improvements
1 parent 3654a19 commit 0475b42

File tree

9 files changed

+116
-53
lines changed

9 files changed

+116
-53
lines changed

‎BVP_2012/BVP_2012_07_17/BVP/BVP/BVP.vcxproj‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -225,6 +225,7 @@ copy /Y "..\..\..\..\Math_C++\mpfr_mpir_x86_x64_msvc2010\mpir\dll\x64\Release\mp
225225
<ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">true</ExcludedFromBuild>
226226
</ClInclude>
227227
<ClInclude Include="MultipleShooting\FourDiagonalSweepMethodStruct.h" />
228+
<ClInclude Include="Problems\AutonomousOscillatingProblem.h" />
228229
<ClInclude Include="Problems\NonAutonomousOscillatingProblem.h" />
229230
<ClInclude Include="Problems\NonAutonomousTroeschProblem.h" />
230231
<ClInclude Include="Problems\OscillatingTestProblem.h" />

‎BVP_2012/BVP_2012_07_17/BVP/BVP/BVP.vcxproj.filters‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -125,5 +125,8 @@
125125
<ClInclude Include="Problems\NonAutonomousOscillatingProblem.h">
126126
<Filter>Header Files\Problems</Filter>
127127
</ClInclude>
128+
<ClInclude Include="Problems\AutonomousOscillatingProblem.h">
129+
<Filter>Header Files\Problems</Filter>
130+
</ClInclude>
128131
</ItemGroup>
129132
</Project>

‎BVP_2012/BVP_2012_07_17/BVP/BVP/MultipleShooting/HybridMultipleShootingComponent.h‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -254,7 +254,7 @@ class HybridMultipleShootingComponent
254254
///which, in turn, depend on the number of knots
255255
auto achievablePrec = GetAcheivablePrecision(MD.size());
256256

257-
for (int i = 0; i < 10 && absCorrection > achievablePrec; i++)
257+
for (int i = 0; i < MaxNumberOfNewtonIterations && absCorrection > achievablePrec; i++)
258258
{
259259
if (absCorrection*auxutils::RoughSqrt((T)MD.size()) < 1) // the process is starting to converge
260260
{
@@ -612,8 +612,9 @@ class HybridMultipleShootingComponent
612612
//auxutils::SaveToFile(result.U, "F:\\V.txt");
613613
}
614614
*/
615-
public:
616615
std::vector<std::pair<int, T>> _debugData;
616+
public:
617+
int MaxNumberOfNewtonIterations;
617618

618619
//Method to save debug data in the specified file (as text)
619620
void SaveDebugData(const char* filename)
@@ -634,6 +635,7 @@ class HybridMultipleShootingComponent
634635
{
635636
_problem = &problem;
636637
_precision = 10 * std::numeric_limits<T>::epsilon();
638+
MaxNumberOfNewtonIterations = 10;
637639
}
638640

639641
vector<InitCondition<T>> Run(

‎BVP_2012/BVP_2012_07_17/BVP/BVP/Problems/AutonomousOscillatingProblem.h‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
#ifndef GUARD_AUTONOMOUS_OSCILLATING_PROBLEM
2+
#define GUARD_AUTONOMOUS_OSCILLATING_PROBLEM
3+
4+
#include "ProblemAutonomousAbstract.h"
5+
#include "../Utils/AuxUtils.h"
6+
7+
template <class T>
8+
class AutonomousOscillatingProblem : public ProblemAutonomousAbstract<T>
9+
{
10+
protected:
11+
virtual T Nonlin(const T& u) override
12+
{
13+
T ln = log(u);
14+
return 1 - ln*(ln + 1);
15+
}
16+
17+
///Derivative of nonlinearity (with respect to u)
18+
virtual T dNonlin(const T& u) override
19+
{
20+
return - (1 + 2*log(u))/u;
21+
}
22+
23+
///Second derivative of nonlinearity
24+
virtual T ddNonlin(const T& u) override
25+
{
26+
return (2*log(u) - 1)/auxutils::sqr(u);
27+
}
28+
};
29+
30+
#endif

‎BVP_2012/BVP_2012_07_17/BVP/BVP/Problems/NonAutonomousOscillatingProblem.h‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -14,31 +14,31 @@ class NonAutonomousOscillatingProblem : public ProblemNonAutonomousAbstract<T>
1414
}
1515

1616
///Derivative of nonlinearity (with respect to u)
17-
virtual T dNonlinDu(const T& u, const T& x)
17+
virtual T dNonlinDu(const T& u, const T& x) override
1818
{
1919
return - (sin(x) + 1)/u;
2020
}
2121

2222
///Derivative of nonlinearity (with respect to x)
23-
virtual T dNonlinDx(const T& u, const T& x)
23+
virtual T dNonlinDx(const T& u, const T& x) override
2424
{
2525
return - cos(x)*log(u);
2626
}
2727

2828
///Second derivative of nonlinearity
29-
virtual T ddNonlinDuDu(const T& u, const T& x)
29+
virtual T ddNonlinDuDu(const T& u, const T& x) override
3030
{
3131
return (sin(x) + 1)/auxutils::sqr(u);
3232
}
3333

3434
///Second derivative of nonlinearity
35-
virtual T ddNonlinDuDx(const T& u, const T& x)
35+
virtual T ddNonlinDuDx(const T& u, const T& x) override
3636
{
3737
return -cos(x)/u;
3838
}
3939

4040
///Second derivative of nonlinearity
41-
virtual T ddNonlinDxDx(const T& u, const T& x)
41+
virtual T ddNonlinDxDx(const T& u, const T& x) override
4242
{
4343
return sin(x)*log(u);
4444
}

‎BVP_2012/BVP_2012_07_17/BVP/BVP/Problems/NonAutonomousTroeschProblem.h‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -26,31 +26,31 @@ class NonAutonomousTroeschProblem : public ProblemNonAutonomousAbstract<T>
2626
}
2727

2828
///Derivative of nonlinearity (with respect to u)
29-
virtual T dNonlinDu(const T& u, const T& x)
29+
virtual T dNonlinDu(const T& u, const T& x) override
3030
{
3131
return Sinhc<T>::Deriv(u, l)*sinh(x*l);
3232
}
3333

3434
///Derivative of nonlinearity (with respect to x)
35-
virtual T dNonlinDx(const T& u, const T& x)
35+
virtual T dNonlinDx(const T& u, const T& x) override
3636
{
3737
return l*Sinhc<T>::Func(u, l)*cosh(x*l);
3838
}
3939

4040
///Second derivative of nonlinearity
41-
virtual T ddNonlinDuDu(const T& u, const T& x)
41+
virtual T ddNonlinDuDu(const T& u, const T& x) override
4242
{
4343
return Sinhc<T>::DDeriv(u, l)*sinh(x*l);
4444
}
4545

4646
///Second derivative of nonlinearity
47-
virtual T ddNonlinDuDx(const T& u, const T& x)
47+
virtual T ddNonlinDuDx(const T& u, const T& x) override
4848
{
4949
return l*Sinhc<T>::Deriv(u, l)*cosh(x*l);
5050
}
5151

5252
///Second derivative of nonlinearity
53-
virtual T ddNonlinDxDx(const T& u, const T& x)
53+
virtual T ddNonlinDxDx(const T& u, const T& x) override
5454
{
5555
return l*l*Sinhc<T>::Func(u, l)*sinh(x*l);
5656
}

‎BVP_2012/BVP_2012_07_17/BVP/GeneralTest/NonAutonomousProblemsTest.cpp‎

Lines changed: 2 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -44,47 +44,10 @@ namespace GeneralTest
4444

4545
//}
4646

47-
TEST_METHOD(NonAutonomousOscilatingProblemMultimpeShootingDouble)
47+
TEST_METHOD(StandardNonAutonomousOscilatingProblemMultimpeShootingDouble)
4848
{
49-
numType preH = 0.1;
50-
numType finalH = 0.001;
51-
numType targetValue = 0.5804096620;
52-
numType targetArgument = 10;
5349
NonAutonomousOscillatingProblem<numType> problem;
54-
55-
std::function<bool(const InitCondition<numType>&)> checkFunc =
56-
[=](const InitCondition<numType>& ic) { return (abs(ic.Value) <= targetArgument)
57-
&& (abs(ic.Argument) <= targetArgument); };
58-
59-
HybridCannon<numType> cannon(problem, preH, preH/10.0, checkFunc);
60-
61-
std::function<int(const InitCondition<numType>&)> evalFunc =
62-
[=](const InitCondition<numType>& ic) { return sgn(ic.Value - targetValue); };
63-
BisectionComponent<numType> bc(cannon);
64-
Assert::IsTrue(bc.DerivativeBisectionGen(0.0, targetArgument, 1.0, 100, 0.9, 1.05, evalFunc));
65-
66-
auto knots = cannon.GetKnotVectorStreight();
67-
knots[knots.size() - 1].Value = targetValue;
68-
knots[knots.size() - 1].Argument = targetArgument;
69-
70-
HybridMultipleShootingComponent<numType> HMSComp(problem);
71-
72-
bool succeeded;
73-
std::vector<InitCondition<numType>> solution = HMSComp.Run(knots, finalH, succeeded);
74-
75-
numType maxDev = CalcDeviationFromExactSolution<numType>(solution,
76-
[](const numType& u){ return exp(sin(u)); });
77-
78-
numType vaxDistanceBetweenKnots = CalcMaxSquaredDistanceBetweenNeighbourKnots(solution);
79-
Assert::IsTrue(vaxDistanceBetweenKnots < 2*finalH*finalH,
80-
Message("Too big maximal distance between neighbour knots " +
81-
auxutils::ToString(vaxDistanceBetweenKnots)));
82-
Assert::IsTrue(maxDev < finalH*finalH,
83-
Message("Too big deviation, maxDev= " + auxutils::ToString(maxDev)));
84-
Assert::IsTrue(solution.size() > targetArgument/finalH,
85-
Message("Too few knot in the solution vector, " +
86-
auxutils::ToString(solution.size())));
50+
StandardOscillatinProblemMultipleShoothingTest<numType, NonAutonomousOscillatingProblem<numType>>(problem);
8751
}
88-
8952
};
9053
}

‎BVP_2012/BVP_2012_07_17/BVP/GeneralTest/OscillatingProblemTestUnit.cpp‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "UnitTestAux.h"
44
#include "../BVP/Utils/AuxUtils.h"
55
#include "..\BVP\Problems\OscillatingTestProblem.h"
6+
#include "../BVP/Problems/AutonomousOscillatingProblem.h"
67
#include "..\BVP\Cannon\HybridCannon.h"
78
#include "..\BVP\ShootingSimple\BisectionComponent.h"
89
#include "..\BVP\MultipleShooting\HybridMultipleShootingComponent.h"
@@ -103,5 +104,13 @@ namespace GeneralTest
103104
Assert::IsTrue(abs(solution[0].Derivative - 1.2576169315833) < 100*std::numeric_limits<double>::epsilon(),
104105
Message("Derivative mismatch " + auxutils::ToString(solution[0].Derivative)));
105106
}
107+
108+
TEST_METHOD(StandardAutonomousOscilatingProblemMultimpeShootingDouble)
109+
{
110+
AutonomousOscillatingProblem<double> problem;
111+
StandardOscillatinProblemMultipleShoothingTest<double,
112+
AutonomousOscillatingProblem<double>>(problem);
113+
}
114+
106115
};
107116
}

‎BVP_2012/BVP_2012_07_17/BVP/GeneralTest/UnitTestAux.h‎

Lines changed: 57 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -23,19 +23,23 @@
2323
#include <functional>
2424
#include <vector>
2525
#include "../BVP/FunctionApproximation/InitialCondition.h"
26+
#include "../BVP/Problems/ProblemAbstract.h"
2627

2728
namespace UnitTestAux
2829
{
2930
///Method to calculate the deviation of the given set of knots to the given "exact solution"
3031
template<class T>
3132
T CalcDeviationFromExactSolution(std::vector<InitCondition<T>> knots,
32-
std::function<T(const T&)> exactSolution)
33+
std::function<T(const T&)> exactSolution, bool checkDerivative = false)
3334
{
3435
T deviation = 0;
3536
for (std::vector<InitCondition<T>>::const_iterator iter = knots.begin(); iter!= knots.end(); ++iter)
3637
{
3738
InitCondition<T> knot = (*iter);
38-
deviation = max(deviation, abs(exactSolution(knot.Argument) - knot.Value));
39+
if (checkDerivative)
40+
deviation = max(deviation, abs(exactSolution(knot.Argument) - knot.Derivative));
41+
else
42+
deviation = max(deviation, abs(exactSolution(knot.Argument) - knot.Value));
3943
}
4044

4145
return deviation;
@@ -65,7 +69,58 @@ namespace UnitTestAux
6569
return maxDistance;
6670
}
6771

72+
///A "Standard test" for a problem whose solution is exp(sin(x));
73+
///The "problem" can be either autonomous or not
74+
template<class T,class P>
75+
void StandardOscillatinProblemMultipleShoothingTest(P problem)
76+
{
77+
T preH = 0.1;
78+
T finalH = 0.001;
79+
T targetValue = 0.5804096620;
80+
T targetArgument = 10;
81+
82+
std::function<bool(const InitCondition<T>&)> checkFunc =
83+
[=](const InitCondition<T>& ic) { return (abs(ic.Value) <= targetArgument)
84+
&& (abs(ic.Argument) <= targetArgument); };
85+
86+
HybridCannon<T> cannon(problem, preH, preH/10.0, checkFunc);
87+
88+
std::function<int(const InitCondition<T>&)> evalFunc =
89+
[=](const InitCondition<T>& ic) { return sgn(ic.Value - targetValue); };
90+
BisectionComponent<T> bc(cannon);
91+
Assert::IsTrue(bc.DerivativeBisectionGen(0.0, targetArgument, 1.0, 100, 0.95, 1.01, evalFunc));
92+
93+
auto knots = cannon.GetKnotVectorStreight();
94+
knots[knots.size() - 1].Value = targetValue;
95+
knots[knots.size() - 1].Argument = targetArgument;
96+
97+
HybridMultipleShootingComponent<T> HMSComp(problem);
98+
HMSComp.MaxNumberOfNewtonIterations = 20;
6899

100+
bool succeeded;
101+
std::vector<InitCondition<T>> solution = HMSComp.Run(knots, finalH, succeeded);
102+
103+
Assert::IsTrue(succeeded, Message("Multiple shoothing method has failed"));
104+
105+
T maxDev = CalcDeviationFromExactSolution<T>(solution,
106+
[](const T& u){ return exp(sin(u)); });
107+
108+
T maxDerivativeDev = CalcDeviationFromExactSolution<T>(solution,
109+
[](const T& u){ return cos(u)*exp(sin(u)); }, true);
110+
111+
T vaxDistanceBetweenKnots = CalcMaxSquaredDistanceBetweenNeighbourKnots(solution);
112+
Assert::IsTrue(vaxDistanceBetweenKnots < 2*finalH*finalH,
113+
Message("Too big maximal distance between neighbour knots " +
114+
auxutils::ToString(vaxDistanceBetweenKnots)));
115+
Assert::IsTrue(maxDev < finalH*finalH,
116+
Message("Too big deviation, maxDev= " + auxutils::ToString(maxDev)));
117+
Assert::IsTrue(maxDerivativeDev < 2*finalH*finalH,
118+
Message("Too big deviation between derivatives, maxDev= " +
119+
auxutils::ToString(maxDev)));
120+
Assert::IsTrue(solution.size() > targetArgument/finalH,
121+
Message("Too few knot in the solution vector, " +
122+
auxutils::ToString(solution.size())));
123+
}
69124

70125
static wchar_t* Message(const char* text)
71126
{

0 commit comments

Comments
 (0)