00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #pragma once
00014
00015 #ifndef __NONLINEARSOLUTIONPATH__
00016 #define __NONLINEARSOLUTIONPATH__
00017
00018
00019 #include "DynamicallyConstructedPath.h"
00020 #include "alg_types.h"
00021 #include "NonLinearGroupElement.h"
00022
00023
00025
00026 template <typename my_alg_type_IN, typename my_alg_type_OUT>
00027 class NonLinearSolutionPath :
00028 public DynamicallyConstructedPath<my_alg_type_OUT>
00029 {
00030 private:
00032 Path<my_alg_type_IN> & theControl;
00033
00035 const std::vector<POLYLIE<my_alg_type_OUT>> theVectorFields;
00036
00038 const SCA Inf;
00039
00041 const AbstractSolutionPoint<my_alg_type_OUT> Ini;
00042
00043 public:
00045 NonLinearSolutionPath(Path<my_alg_type_IN> & theControlIn, std::vector<POLYLIE<my_alg_type_OUT>> & theVectorFieldsIn, AbstractSolutionPoint<my_alg_type_OUT> & theInitialValue, SCA theStartTime = SCA(0))
00046 : theControl(theControlIn), theVectorFields(theVectorFieldsIn), Inf(theStartTime), Ini(theInitialValue)
00047 {
00048 };
00049
00051 ~NonLinearSolutionPath(void)
00052 {
00053 };
00054
00056 void ComputeChildLieIncrements( Increment<my_alg_type_OUT> & nvLeft, Increment<my_alg_type_OUT> & nvRight, ConstIterator itLeafAbove ) const
00057 {
00058 const dyadic_interval & diAbove = itLeafAbove->first;
00059 const Increment<my_alg_type_OUT> & nvAbove = itLeafAbove->second;
00060
00061
00062 SCA midPoint = SCA(diAbove.inf() + ((diAbove.sup() - diAbove.inf())*0.5));
00063 LIE incLeft = ComputeIncrement( SCA(diAbove.inf()), midPoint );
00064 LIE incRight = nvAbove.LieValue() - incLeft;
00065 nvLeft.LieValue(incLeft);
00066 nvRight.LieValue(incRight);
00067 };
00068
00070 LIE MakeRootLieIncrement( const dyadic_interval &increment ) const
00071 {
00072 return ComputeIncrement( SCA(increment.inf()), SCA(increment.sup()) );
00073 };
00074
00076 void MakeNeighborRootLieIncrement( LIE & ans, const Iterator & OldRoot ) const
00077 {
00079 dyadic_interval diIncrement = OldRoot->first;
00080 diIncrement.flip_interval();
00081 ans = ComputeIncrement( diIncrement.inf(), diIncrement.sup() );
00082 };
00083
00084 private:
00085 LIE ComputeIncrement(const SCA t0, const SCA t1, const int accuracy=11) const
00086 {
00087 if (t0 == Inf)
00088 {
00089
00090 AbstractSolutionPoint<my_alg_type_OUT> solution;
00091 solution = ComputeSolution(Inf, Ini, t1, accuracy);
00092
00093
00094 AbstractSolutionPoint<my_alg_type_OUT> inc;
00095 AbstractSolutionPoint<my_alg_type_OUT> temp(Ini);
00096 inc = solution - temp;
00097
00098 LIE result;
00099
00100 for(unsigned int i = 1; i <= my_alg_type_OUT::myDIM; i++)
00101 {
00102 result += LIE(LET(i),inc[LET(i)]);
00103 }
00104
00105 return result;
00106 }
00107 else
00108 {
00109 if (t0 > Inf)
00110 {
00111
00112
00113
00114 AbstractSolutionPoint<my_alg_type_OUT> a0;
00115 a0 = ComputeSolution(Inf, Ini, t0, accuracy);
00116
00117
00118 AbstractSolutionPoint<my_alg_type_OUT> solution;
00119 solution = ComputeSolution(t0, a0, t1, accuracy);
00120
00121
00122 AbstractSolutionPoint<my_alg_type_OUT> inc;
00123 inc = solution - a0;
00124
00125 LIE result;
00126
00127 for(unsigned int i = 1; i <= my_alg_type_OUT::myDIM; i++)
00128 {
00129 result += LIE(LET(i),inc[LET(i)]);
00130 }
00131
00132 return result;
00133 }
00134 else
00135 {
00136
00137 std::cout << "Error (NonLinearSolutionPath): Time less than start time!" << std::endl;
00138 abort();
00139 return LIE();
00140 }
00141 }
00142 };
00143
00145 AbstractSolutionPoint<my_alg_type_OUT> ComputeSolution(SCA t0, const AbstractSolutionPoint<my_alg_type_OUT> & a0, SCA t1, const unsigned int my_accuracy) const
00146 {
00147 dyadic_interval Width(t1-t0, my_accuracy);
00148 dyadic_interval Begin = dyadic_interval::dyadic_bracket(t0,Width.n);
00149 dyadic_interval End = dyadic_interval::dyadic_bracket(t1,Width.n);
00150
00151 unsigned int NoIncrements = End.k - Begin.k;
00152
00153 SPB::NonLinearGroupElement<my_alg_type_IN,my_alg_type_OUT> CurrentGroupElt;
00154
00155 AbstractSolutionPoint<my_alg_type_OUT> CurrentValue;
00156 CurrentValue = a0;
00157
00158 AbstractSolutionPoint<my_alg_type_OUT> temp;
00159
00160 dyadic_interval CurrentInterval=Begin;
00161
00162 lietovectorfield<my_alg_type_IN,POLYLIE_OUT> l2vf(theVectorFields);
00163
00164 for(unsigned int i=1; i <= NoIncrements; i++)
00165 {
00166 CurrentGroupElt = SPB::NonLinearGroupElement<my_alg_type_IN,my_alg_type_OUT>(theControl.DescribePath(CurrentInterval++), l2vf);
00167
00168 temp = CurrentGroupElt.evaluate(CurrentValue);
00169 CurrentValue = temp;
00170 }
00171
00172 return CurrentValue;
00173 };
00174
00175 };
00176
00177 #endif // __NONLINEARSOLUTIONPATH__