Add Cranck-Nicholson, Method of Lines & TrBDF2 schemes (amaggiulli#244)

* Update step method & add new schemes

- Update step method by adding theta parameter.
- Add Cranck-Nicholson, Method of Lines & TrBDF2 schemes

* Add new schemes to FdmBackwardSolver

* Update IMixedScheme interface and dependencies

* Update QLNet Old project file

* Delete FdmBackwardSolver.cs

* Move to Solvers folder

src/QLNet.Old/QLNet.Old.csproj

@@ -1092,6 +1092,9 @@
     </Compile>
     <Compile Include="..\QLNet\Methods\Finitedifferences\Schemes\CraigSneydScheme.cs">
       <Link>Methods\Finitedifferences\Schemes\CraigSneydScheme.cs</Link>
+    </Compile>
+	<Compile Include="..\QLNet\Methods\Finitedifferences\Schemes\CrankNicolsonScheme.cs">
+      <Link>Methods\Finitedifferences\Schemes\CrankNicolsonScheme.cs</Link>
     </Compile>
     <Compile Include="..\QLNet\Methods\Finitedifferences\Schemes\DouglasScheme.cs">
       <Link>Methods\Finitedifferences\Schemes\DouglasScheme.cs</Link>
@@ -1104,9 +1107,15 @@
     </Compile>
     <Compile Include="..\QLNet\Methods\Finitedifferences\Schemes\ImplicitEulerScheme.cs">
       <Link>Methods\Finitedifferences\Schemes\ImplicitEulerScheme.cs</Link>
+    </Compile>
+	<Compile Include="..\QLNet\Methods\Finitedifferences\Schemes\MethodOfLinesScheme.cs">
+      <Link>Methods\Finitedifferences\Schemes\MethodOfLinesScheme.cs</Link>
     </Compile>
     <Compile Include="..\QLNet\Methods\Finitedifferences\Schemes\ModifiedCraigSneydScheme.cs">
       <Link>Methods\Finitedifferences\Schemes\ModifiedCraigSneydScheme.cs</Link>
+    </Compile>
+	<Compile Include="..\QLNet\Methods\Finitedifferences\Schemes\TrBDF2Scheme.cs">
+      <Link>Methods\Finitedifferences\Schemes\TrBDF2Scheme.cs</Link>
     </Compile>
     <Compile Include="..\QLNet\Methods\Finitedifferences\Solvers\Fdm1DimSolver.cs">
       <Link>Methods\Finitedifferences\Solvers\Fdm1DimSolver.cs</Link>

src/QLNet/Methods/Finitedifferences/MixedScheme.cs

@@ -27,7 +27,7 @@ public interface ISchemeFactory
 
    public interface IMixedScheme
    {
-      void step(ref object a, double t);
+      void step(ref object a, double t, double theta = 1.0);
       void setStep(double dt);
    }
 
@@ -55,7 +55,7 @@ public MixedScheme(Operator L, double theta, List<BoundaryCondition<IOperator>>
          bcs_ = bcs;
       }
 
-      public void step(ref object o, double t)
+      public void step(ref object o, double t, double theta = 1.0)
       {
          Vector a = (Vector)o;
 

src/QLNet/Methods/Finitedifferences/ParallelEvolver.cs

@@ -62,7 +62,7 @@ public ParallelEvolver(List<IOperator> L, BoundaryConditionSet bcs)
          }
       }
 
-      public void step(ref object o, double t)
+      public void step(ref object o, double t, double theta = 1.0)
       {
          List<Vector> a = (List<Vector>)o;
          for (int i = 0; i < evolvers_.Count; i++)

src/QLNet/Methods/Finitedifferences/Schemes/CraigSneydScheme.cs

@@ -55,7 +55,7 @@ public IMixedScheme factory(object L, object bcs, object[] additionalInputs = nu
 
       #region IMixedScheme interface
 
-      public void step(ref object a, double t)
+      public void step(ref object a, double t, double theta = 1.0)
       {
          Utils.QL_REQUIRE(t - dt_.Value > -1e-8, () => "a step towards negative time given");
 

src/QLNet/Methods/Finitedifferences/Schemes/CrankNicolsonScheme.cs

@@ -0,0 +1,85 @@
+/*
+ Copyright (C) 2020 Jean-Camille Tournier ([email protected])
+
+ This file is part of QLNet Project https://github.com/amaggiulli/qlnet
+
+ QLNet is free software: you can redistribute it and/or modify it
+ under the terms of the QLNet license.  You should have received a
+ copy of the license along with this program; if not, license is
+ available online at <http://qlnet.sourceforge.net/License.html>.
+
+ QLNet is a based on QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+ The QuantLib license is available online at http://quantlib.org/license.shtml.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+using System;
+using System.Collections.Generic;
+using System.Linq;
+
+namespace QLNet
+{
+   /*! In one dimension the Crank-Nicolson scheme is equivalent to the
+       Douglas scheme and in higher dimensions it is usually inferior to
+       operator splitting methods like Craig-Sneyd or Hundsdorfer-Verwer.
+   */
+   public class CrankNicolsonScheme : IMixedScheme, ISchemeFactory
+   {
+      public CrankNicolsonScheme()
+      { }
+
+      public CrankNicolsonScheme(double theta,
+                                 FdmLinearOpComposite map,
+                                 List<BoundaryCondition<FdmLinearOp>> bcSet = null,
+                                 double relTol = 1E-8,
+                                 ImplicitEulerScheme.SolverType solverType = ImplicitEulerScheme.SolverType.BiCGstab)
+      {
+         dt_ = null;
+         theta_ = theta;
+         explicit_ = new ExplicitEulerScheme(map, bcSet);
+         implicit_ = new ImplicitEulerScheme(map, bcSet, relTol, solverType);
+      }
+
+      #region ISchemeFactory
+
+      public IMixedScheme factory(object L, object bcs, object[] additionalInputs = null)
+      {
+         double? theta = additionalInputs[0] as double?;
+         double? relTol = additionalInputs[1] as double?;
+         ImplicitEulerScheme.SolverType? solverType = additionalInputs[2] as ImplicitEulerScheme.SolverType?;
+         return new CrankNicolsonScheme(theta.Value, L as FdmLinearOpComposite,
+                                        bcs as List<BoundaryCondition<FdmLinearOp>>, relTol.Value, solverType.Value);
+      }
+
+      #endregion
+
+      public void step(ref object a, double t, double theta = 1.0)
+      {
+         Utils.QL_REQUIRE(t - dt_ > -1e-8, () => "a step towards negative time given");
+         if (theta_ != 1.0)
+            explicit_.step(ref a, t, 1.0 - theta_);
+
+         if (theta_ != 0.0)
+            implicit_.step(ref a, t, theta_);
+      }
+
+      public void setStep(double dt)
+      {
+         dt_ = dt;
+         explicit_.setStep(dt_.Value);
+         implicit_.setStep(dt_.Value);
+      }
+
+      public int numberOfIterations()
+      {
+         return implicit_.numberOfIterations();
+      }
+      protected double? dt_;
+      protected double theta_;
+      protected ExplicitEulerScheme explicit_;
+      protected ImplicitEulerScheme implicit_;
+   }
+}

src/QLNet/Methods/Finitedifferences/Schemes/DouglasScheme.cs

@@ -53,7 +53,7 @@ public IMixedScheme factory(object L, object bcs, object[] additionalInputs = nu
 
       #region IMixedScheme interface
 
-      public void step(ref object a, double t)
+      public void step(ref object a, double t, double theta = 1.0)
       {
          Utils.QL_REQUIRE(t - dt_.Value > -1e-8, () => "a step towards negative time given");
          map_.setTime(Math.Max(0.0, t - dt_.Value), t);

src/QLNet/Methods/Finitedifferences/Schemes/ExplicitEulerScheme.cs

@@ -42,16 +42,15 @@ public IMixedScheme factory(object L, object bcs, object[] additionalInputs = nu
       #endregion
 
       #region IMixedScheme interface
-      public void step(ref object a, double t)
+      public void step(ref object a, double t, double theta = 1.0)
       {
          Utils.QL_REQUIRE(t - dt_ > -1e-8, () => "a step towards negative time given");
          map_.setTime(Math.Max(0.0, t - dt_.Value), t);
          bcSet_.setTime(Math.Max(0.0, t - dt_.Value));
 
          bcSet_.applyBeforeApplying(map_);
-         a = (a as Vector) + dt_.Value * map_.apply(a as Vector);
+         a = (a as Vector) + (theta * dt_.Value) * map_.apply(a as Vector);
          bcSet_.applyAfterApplying(a as Vector);
-
       }
 
       public void setStep(double dt)

src/QLNet/Methods/Finitedifferences/Schemes/HundsdorferScheme.cs

@@ -55,7 +55,7 @@ public IMixedScheme factory(object L, object bcs, object[] additionalInputs = nu
 
       #region IMixedScheme interface
 
-      public void step(ref object a, double t)
+      public void step(ref object a, double t, double theta = 1.0)
       {
          Utils.QL_REQUIRE(t - dt_.Value > -1e-8, () => "a step towards negative time given");
 

src/QLNet/Methods/Finitedifferences/Schemes/ImplicitEulerScheme.cs

@@ -49,33 +49,45 @@ public IMixedScheme factory(object L, object bcs, object[] additionalFields = nu
       #endregion
 
       #region IMixedScheme interface
-      public void step(ref object a, double t)
+      public void step(ref object a, double t, double theta = 1.0)
       {
          Utils.QL_REQUIRE(t - dt_.Value > -1e-8, () => "a step towards negative time given");
          map_.setTime(Math.Max(0.0, t - dt_.Value), t);
          bcSet_.setTime(Math.Max(0.0, t - dt_.Value));
 
          bcSet_.applyBeforeSolving(map_, a as Vector);
 
-         if (solverType_ == SolverType.BiCGstab)
+         if (map_.size() == 1)
          {
-            BiCGStabResult result =
-               new BiCGStab(this.apply, Math.Max(10, (a as Vector).Count), relTol_, x => map_.preconditioner(x, -dt_.Value)).solve(a as Vector, a as Vector);
-
-            iterations_ += result.Iterations;
-            a = result.X;
+            a = map_.solve_splitting(0, a as Vector, -theta * dt_.Value);
          }
-         else if (solverType_ == SolverType.GMRES)
+         else
          {
-            GMRESResult result =
-               new GMRES(this.apply, Math.Max(10, (a as Vector).Count) / 10, relTol_, x => map_.preconditioner(x, -dt_.Value)).solve(a as Vector, a as Vector);
-
-            iterations_ += result.Errors.Count;
-            a = result.X;
+            if (solverType_ == SolverType.BiCGstab)
+            {
+               BiCGStab.MatrixMult preconditioner = x => map_.preconditioner(x, -theta * dt_.Value);
+               BiCGStab.MatrixMult applyF = x => this.apply(x, theta);
+
+               BiCGStabResult result =
+                  new BiCGStab(applyF, Math.Max(10, (a as Vector).Count), relTol_, preconditioner).solve(a as Vector, a as Vector);
+
+               iterations_ += result.Iterations;
+               a = result.X;
+            }
+            else if (solverType_ == SolverType.GMRES)
+            {
+               GMRES.MatrixMult preconditioner = x => map_.preconditioner(x, -theta * dt_.Value);
+               GMRES.MatrixMult applyF = x => this.apply(x, theta);
+
+               GMRESResult result =
+                  new GMRES(applyF, Math.Max(10, (a as Vector).Count) / 10, relTol_, preconditioner).solve(a as Vector, a as Vector);
+
+               iterations_ += result.Errors.Count;
+               a = result.X;
+            }
+            else
+               Utils.QL_FAIL("unknown/illegal solver type");
          }
-         else
-            Utils.QL_FAIL("unknown/illegal solver type");
-
          bcSet_.applyAfterSolving(a as Vector);
       }
 
@@ -90,9 +102,9 @@ public int numberOfIterations()
          return iterations_;
       }
 
-      protected Vector apply(Vector r)
+      protected Vector apply(Vector r, double theta = 1.0)
       {
-         return r - dt_.Value * map_.apply(r);
+         return r - (theta * dt_.Value) * map_.apply(r);
       }
 
       protected double? dt_;

src/QLNet/Methods/Finitedifferences/Schemes/MethodOfLinesScheme.cs

@@ -0,0 +1,87 @@
+/*
+ Copyright (C) 2020 Jean-Camille Tournier ([email protected])
+
+ This file is part of QLNet Project https://github.com/amaggiulli/qlnet
+
+ QLNet is free software: you can redistribute it and/or modify it
+ under the terms of the QLNet license.  You should have received a
+ copy of the license along with this program; if not, license is
+ available online at <http://qlnet.sourceforge.net/License.html>.
+
+ QLNet is a based on QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+ The QuantLib license is available online at http://quantlib.org/license.shtml.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+using System;
+using System.Collections.Generic;
+using System.Linq;
+
+namespace QLNet
+{
+   /*! In one dimension the Crank-Nicolson scheme is equivalent to the
+       Douglas scheme and in higher dimensions it is usually inferior to
+       operator splitting methods like Craig-Sneyd or Hundsdorfer-Verwer.
+   */
+   public class MethodOfLinesScheme : IMixedScheme, ISchemeFactory
+   {
+      public MethodOfLinesScheme()
+      { }
+
+      public MethodOfLinesScheme(double eps,
+                                 double relInitStepSize,
+                                 FdmLinearOpComposite map,
+                                 List<BoundaryCondition<FdmLinearOp>> bcSet = null)
+      {
+         dt_ = null;
+         eps_ = eps;
+         relInitStepSize_ = relInitStepSize;
+         map_ = map;
+         bcSet_ = new BoundaryConditionSchemeHelper(bcSet);
+      }
+
+      #region ISchemeFactory
+
+      public IMixedScheme factory(object L, object bcs, object[] additionalInputs = null)
+      {
+         double? eps = additionalInputs[0] as double?;
+         double? relInitStepSize = additionalInputs[1] as double?;
+         return new MethodOfLinesScheme(eps.Value, relInitStepSize.Value,
+                                        L as FdmLinearOpComposite, bcs as List<BoundaryCondition<FdmLinearOp>>);
+      }
+
+      #endregion
+
+      protected List<double> apply(double t, List<double> r)
+      {
+         map_.setTime(t, t + 0.0001);
+         bcSet_.applyBeforeApplying(map_);
+
+         Vector dxdt = -1.0 * map_.apply(new Vector(r));
+
+         return dxdt;
+      }
+
+      public void step(ref object a, double t, double theta = 1.0)
+      {
+         Utils.QL_REQUIRE(t - dt_ > -1e-8, () => "a step towards negative time given");
+         List<double> v = new AdaptiveRungeKutta(eps_, relInitStepSize_ * dt_.Value).value(this.apply, a as Vector, t, Math.Max(0.0, t - dt_.Value));
+         Vector y = new Vector(v);
+         bcSet_.applyAfterSolving(y);
+         a = y;
+      }
+
+      public void setStep(double dt)
+      {
+         dt_ = dt;
+      }
+
+      protected double? dt_;
+      protected double eps_, relInitStepSize_;
+      protected FdmLinearOpComposite map_;
+      protected BoundaryConditionSchemeHelper bcSet_;
+   }
+}

src/QLNet/Methods/Finitedifferences/Schemes/ModifiedCraigSneydScheme.cs

@@ -55,7 +55,7 @@ public IMixedScheme factory(object L, object bcs, object[] additionalInputs)
       #endregion
 
       #region IMixedScheme interface
-      public void step(ref object a, double t)
+      public void step(ref object a, double t, double theta = 1.0)
       {
          Utils.QL_REQUIRE(t - dt_.Value > -1e-8, () => "a step towards negative time given");
          map_.setTime(Math.Max(0.0, t - dt_.Value), t);