From 46075c83fa23c0714918684654c6cb08d67e8481 Mon Sep 17 00:00:00 2001 From: tournierjc Date: Wed, 29 Jan 2020 19:53:02 +0100 Subject: [PATCH] Add Cranck-Nicholson, Method of Lines & TrBDF2 schemes (#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 | 9 ++ .../Methods/Finitedifferences/MixedScheme.cs | 4 +- .../Finitedifferences/ParallelEvolver.cs | 2 +- .../Schemes/CraigSneydScheme.cs | 2 +- .../Schemes/CrankNicolsonScheme.cs | 85 +++++++++++ .../Schemes/DouglasScheme.cs | 2 +- .../Schemes/ExplicitEulerScheme.cs | 5 +- .../Schemes/HundsdorferScheme.cs | 2 +- .../Schemes/ImplicitEulerScheme.cs | 48 +++--- .../Schemes/MethodOfLinesScheme.cs | 87 +++++++++++ .../Schemes/ModifiedCraigSneydScheme.cs | 2 +- .../Finitedifferences/Schemes/TrBDF2Scheme.cs | 143 ++++++++++++++++++ .../Solvers/FdmBackwardSolver.cs | 36 ++++- src/QLNet/Methods/Finitedifferences/TRBDF2.cs | 2 +- 14 files changed, 399 insertions(+), 30 deletions(-) create mode 100644 src/QLNet/Methods/Finitedifferences/Schemes/CrankNicolsonScheme.cs create mode 100644 src/QLNet/Methods/Finitedifferences/Schemes/MethodOfLinesScheme.cs create mode 100644 src/QLNet/Methods/Finitedifferences/Schemes/TrBDF2Scheme.cs diff --git a/src/QLNet.Old/QLNet.Old.csproj b/src/QLNet.Old/QLNet.Old.csproj index e3126426b..9fd82e324 100644 --- a/src/QLNet.Old/QLNet.Old.csproj +++ b/src/QLNet.Old/QLNet.Old.csproj @@ -1092,6 +1092,9 @@ Methods\Finitedifferences\Schemes\CraigSneydScheme.cs + + + Methods\Finitedifferences\Schemes\CrankNicolsonScheme.cs Methods\Finitedifferences\Schemes\DouglasScheme.cs @@ -1104,9 +1107,15 @@ Methods\Finitedifferences\Schemes\ImplicitEulerScheme.cs + + + Methods\Finitedifferences\Schemes\MethodOfLinesScheme.cs Methods\Finitedifferences\Schemes\ModifiedCraigSneydScheme.cs + + + Methods\Finitedifferences\Schemes\TrBDF2Scheme.cs Methods\Finitedifferences\Solvers\Fdm1DimSolver.cs diff --git a/src/QLNet/Methods/Finitedifferences/MixedScheme.cs b/src/QLNet/Methods/Finitedifferences/MixedScheme.cs index cee8ab538..c722154bb 100644 --- a/src/QLNet/Methods/Finitedifferences/MixedScheme.cs +++ b/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> 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; diff --git a/src/QLNet/Methods/Finitedifferences/ParallelEvolver.cs b/src/QLNet/Methods/Finitedifferences/ParallelEvolver.cs index 85de78340..71497c7b7 100644 --- a/src/QLNet/Methods/Finitedifferences/ParallelEvolver.cs +++ b/src/QLNet/Methods/Finitedifferences/ParallelEvolver.cs @@ -62,7 +62,7 @@ public ParallelEvolver(List L, BoundaryConditionSet bcs) } } - public void step(ref object o, double t) + public void step(ref object o, double t, double theta = 1.0) { List a = (List)o; for (int i = 0; i < evolvers_.Count; i++) diff --git a/src/QLNet/Methods/Finitedifferences/Schemes/CraigSneydScheme.cs b/src/QLNet/Methods/Finitedifferences/Schemes/CraigSneydScheme.cs index d521b3314..8f5d069f6 100644 --- a/src/QLNet/Methods/Finitedifferences/Schemes/CraigSneydScheme.cs +++ b/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"); diff --git a/src/QLNet/Methods/Finitedifferences/Schemes/CrankNicolsonScheme.cs b/src/QLNet/Methods/Finitedifferences/Schemes/CrankNicolsonScheme.cs new file mode 100644 index 000000000..19d78ca96 --- /dev/null +++ b/src/QLNet/Methods/Finitedifferences/Schemes/CrankNicolsonScheme.cs @@ -0,0 +1,85 @@ +/* + Copyright (C) 2020 Jean-Camille Tournier (jean-camille.tournier@avivainvestors.com) + + 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 . + + 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> 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>, 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_; + } +} diff --git a/src/QLNet/Methods/Finitedifferences/Schemes/DouglasScheme.cs b/src/QLNet/Methods/Finitedifferences/Schemes/DouglasScheme.cs index 3467c798f..002fe0f70 100644 --- a/src/QLNet/Methods/Finitedifferences/Schemes/DouglasScheme.cs +++ b/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); diff --git a/src/QLNet/Methods/Finitedifferences/Schemes/ExplicitEulerScheme.cs b/src/QLNet/Methods/Finitedifferences/Schemes/ExplicitEulerScheme.cs index a0cd8c371..efb1e5078 100644 --- a/src/QLNet/Methods/Finitedifferences/Schemes/ExplicitEulerScheme.cs +++ b/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) diff --git a/src/QLNet/Methods/Finitedifferences/Schemes/HundsdorferScheme.cs b/src/QLNet/Methods/Finitedifferences/Schemes/HundsdorferScheme.cs index 1a2b8dbd3..a5afb26f2 100644 --- a/src/QLNet/Methods/Finitedifferences/Schemes/HundsdorferScheme.cs +++ b/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"); diff --git a/src/QLNet/Methods/Finitedifferences/Schemes/ImplicitEulerScheme.cs b/src/QLNet/Methods/Finitedifferences/Schemes/ImplicitEulerScheme.cs index 9c8abc34c..f3f71c973 100644 --- a/src/QLNet/Methods/Finitedifferences/Schemes/ImplicitEulerScheme.cs +++ b/src/QLNet/Methods/Finitedifferences/Schemes/ImplicitEulerScheme.cs @@ -49,7 +49,7 @@ 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); @@ -57,25 +57,37 @@ public void step(ref object a, double t) 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_; diff --git a/src/QLNet/Methods/Finitedifferences/Schemes/MethodOfLinesScheme.cs b/src/QLNet/Methods/Finitedifferences/Schemes/MethodOfLinesScheme.cs new file mode 100644 index 000000000..8ea27453d --- /dev/null +++ b/src/QLNet/Methods/Finitedifferences/Schemes/MethodOfLinesScheme.cs @@ -0,0 +1,87 @@ +/* + Copyright (C) 2020 Jean-Camille Tournier (jean-camille.tournier@avivainvestors.com) + + 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 . + + 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> 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>); + } + + #endregion + + protected List apply(double t, List 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 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_; + } +} diff --git a/src/QLNet/Methods/Finitedifferences/Schemes/ModifiedCraigSneydScheme.cs b/src/QLNet/Methods/Finitedifferences/Schemes/ModifiedCraigSneydScheme.cs index 865afca10..97bf1188a 100644 --- a/src/QLNet/Methods/Finitedifferences/Schemes/ModifiedCraigSneydScheme.cs +++ b/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); diff --git a/src/QLNet/Methods/Finitedifferences/Schemes/TrBDF2Scheme.cs b/src/QLNet/Methods/Finitedifferences/Schemes/TrBDF2Scheme.cs new file mode 100644 index 000000000..ea3de6fd3 --- /dev/null +++ b/src/QLNet/Methods/Finitedifferences/Schemes/TrBDF2Scheme.cs @@ -0,0 +1,143 @@ +/* + Copyright (C) 2020 Jean-Camille Tournier (jean-camille.tournier@avivainvestors.com) + + 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 . + + 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 TrBDF2Scheme : IMixedScheme, ISchemeFactory + where TrapezoidalScheme : class, IMixedScheme + { + public enum SolverType { BiCGstab, GMRES } + + public TrBDF2Scheme() + { } + + public TrBDF2Scheme(double alpha, + FdmLinearOpComposite map, + TrapezoidalScheme trapezoidalScheme, + List> bcSet = null, + double relTol = 1E-8, + TrBDF2Scheme.SolverType solverType = TrBDF2Scheme.SolverType.BiCGstab) + { + dt_ = null; + beta_ = null; + iterations_ = 0; + alpha_ = alpha; + map_ = map; + trapezoidalScheme_ = trapezoidalScheme; + bcSet_ = new BoundaryConditionSchemeHelper(bcSet); + relTol_ = relTol; + solverType_ = solverType; + } + + #region ISchemeFactory + + public IMixedScheme factory(object L, object bcs, object[] additionalInputs = null) + { + double? alpha = additionalInputs[0] as double?; + TrapezoidalScheme trapezoidalScheme = additionalInputs[1] as TrapezoidalScheme; + double? relTol = additionalInputs[0] as double?; + TrBDF2Scheme.SolverType? solverType = additionalInputs[2] as TrBDF2Scheme.SolverType?; + return new TrBDF2Scheme(alpha.Value, L as FdmLinearOpComposite, trapezoidalScheme, + bcs as List>, 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"); + double intermediateTimeStep = dt_.Value * alpha_; + + object fStar = a; + trapezoidalScheme_.setStep(intermediateTimeStep); + trapezoidalScheme_.step(ref fStar, t); + + bcSet_.setTime(Math.Max(0.0, t - dt_.Value)); + bcSet_.applyBeforeSolving(map_, a as Vector); + + Vector fStarVec = fStar as Vector; + Vector f = (1.0 / alpha_ * fStarVec - Math.Pow(1.0 - alpha_, 2.0) / alpha_ * (a as Vector)) / (2.0 - alpha_); + + if (map_.size() == 1) + { + a = map_.solve_splitting(0, f, -beta_.Value); + } + else + { + if (solverType_ == SolverType.BiCGstab) + { + BiCGStab.MatrixMult preconditioner = x => map_.preconditioner(x, -beta_.Value); + BiCGStab.MatrixMult applyF = x => this.apply(x); + + BiCGStabResult result = + new BiCGStab(applyF, Math.Max(10, (a as Vector).Count), relTol_, preconditioner).solve(f, f); + + iterations_ += result.Iterations; + a = result.X; + } + else if (solverType_ == SolverType.GMRES) + { + GMRES.MatrixMult preconditioner = x => map_.preconditioner(x, -beta_.Value); + GMRES.MatrixMult applyF = x => this.apply(x); + + GMRESResult result = + new GMRES(applyF, Math.Max(10, (a as Vector).Count) / 10, relTol_, preconditioner).solve(f, f); + + iterations_ += result.Errors.Count; + a = result.X; + } + else + Utils.QL_FAIL("unknown/illegal solver type"); + } + bcSet_.applyAfterSolving(a as Vector); + } + + public void setStep(double dt) + { + dt_ = dt; + beta_ = (1.0 - alpha_) / (2.0 - alpha_) * dt_.Value; + } + + public int numberOfIterations() + { + return iterations_; + } + + public Vector apply(Vector r) + { + return r - beta_.Value * map_.apply(r); + } + + protected double? dt_, beta_; + protected double alpha_, relTol_; + protected int iterations_; + + protected FdmLinearOpComposite map_; + protected TrapezoidalScheme trapezoidalScheme_; + protected BoundaryConditionSchemeHelper bcSet_; + protected SolverType solverType_; + } +} diff --git a/src/QLNet/Methods/Finitedifferences/Solvers/FdmBackwardSolver.cs b/src/QLNet/Methods/Finitedifferences/Solvers/FdmBackwardSolver.cs index a3d0bed4d..6a71ef129 100644 --- a/src/QLNet/Methods/Finitedifferences/Solvers/FdmBackwardSolver.cs +++ b/src/QLNet/Methods/Finitedifferences/Solvers/FdmBackwardSolver.cs @@ -27,7 +27,9 @@ public enum FdmSchemeType { HundsdorferType, DouglasType, CraigSneydType, ModifiedCraigSneydType, - ImplicitEulerType, ExplicitEulerType + ImplicitEulerType, ExplicitEulerType, + MethodOfLinesType, TrBDF2Type, + CrankNicolsonType } public FdmSchemeDesc() { } @@ -67,12 +69,15 @@ public double mu // some default scheme descriptions public FdmSchemeDesc Douglas() { return new FdmSchemeDesc(FdmSchemeType.DouglasType, 0.5, 0.0); } + public FdmSchemeDesc CrankNicolson() { return new FdmSchemeDesc(FdmSchemeType.CrankNicolsonType, 0.5, 0.0); } public FdmSchemeDesc ImplicitEuler() { return new FdmSchemeDesc(FdmSchemeType.ImplicitEulerType, 0.0, 0.0); } public FdmSchemeDesc ExplicitEuler() { return new FdmSchemeDesc(FdmSchemeType.ExplicitEulerType, 0.0, 0.0); } public FdmSchemeDesc CraigSneyd() { return new FdmSchemeDesc(FdmSchemeType.CraigSneydType, 0.5, 0.5); } public FdmSchemeDesc ModifiedCraigSneyd() { return new FdmSchemeDesc(FdmSchemeType.ModifiedCraigSneydType, 1.0 / 3.0, 1.0 / 3.0); } public FdmSchemeDesc Hundsdorfer() { return new FdmSchemeDesc(FdmSchemeType.HundsdorferType, 0.5 + Math.Sqrt(3.0) / 6.0, 0.5); } public FdmSchemeDesc ModifiedHundsdorfer() { return new FdmSchemeDesc(FdmSchemeType.HundsdorferType, 1.0 - Math.Sqrt(2.0) / 2.0, 0.5); } + public FdmSchemeDesc MethodOfLines(double eps = 0.001, double relInitStepSize = 0.01) { return new FdmSchemeDesc(FdmSchemeType.MethodOfLinesType, eps, relInitStepSize); } + public FdmSchemeDesc TrBDF2() { return new FdmSchemeDesc(FdmSchemeType.TrBDF2Type, 2.0 - Math.Sqrt(2.0), 1E-8); } } public class FdmBackwardSolver @@ -126,6 +131,14 @@ FiniteDifferenceModel dampingModel dsModel.rollback(ref a, dampingTo, to, steps, condition_); } break; + case FdmSchemeDesc.FdmSchemeType.CrankNicolsonType: + { + CrankNicolsonScheme cnEvolver = new CrankNicolsonScheme(schemeDesc_.theta, map_, bcSet_); + FiniteDifferenceModel + cnModel = new FiniteDifferenceModel(cnEvolver, condition_.stoppingTimes()); + cnModel.rollback(ref a, dampingTo, to, steps, condition_); + } + break; case FdmSchemeDesc.FdmSchemeType.CraigSneydType: { CraigSneydScheme csEvolver = new CraigSneydScheme(schemeDesc_.theta, schemeDesc_.mu, @@ -161,6 +174,27 @@ FiniteDifferenceModel dampingModel explicitModel.rollback(ref a, dampingTo, to, steps, condition_); } break; + case FdmSchemeDesc.FdmSchemeType.MethodOfLinesType: + { + MethodOfLinesScheme methodOfLines = new MethodOfLinesScheme(schemeDesc_.theta, schemeDesc_.mu, map_, bcSet_); + FiniteDifferenceModel + molModel = new FiniteDifferenceModel(methodOfLines, condition_.stoppingTimes()); + molModel.rollback(ref a, dampingTo, to, steps, condition_); + } + break; + case FdmSchemeDesc.FdmSchemeType.TrBDF2Type: + { + FdmSchemeDesc trDesc = new FdmSchemeDesc().CraigSneyd(); + CraigSneydScheme hsEvolver = new CraigSneydScheme(trDesc.theta, trDesc.mu, map_, bcSet_); + + TrBDF2Scheme trBDF2 = new TrBDF2Scheme( + schemeDesc_.theta, map_, hsEvolver, bcSet_, schemeDesc_.mu); + + FiniteDifferenceModel> + trBDF2Model = new FiniteDifferenceModel>(trBDF2, condition_.stoppingTimes()); + trBDF2Model.rollback(ref a, dampingTo, to, steps, condition_); + } + break; default: Utils.QL_FAIL("Unknown scheme type"); break; diff --git a/src/QLNet/Methods/Finitedifferences/TRBDF2.cs b/src/QLNet/Methods/Finitedifferences/TRBDF2.cs index 3b7a97e63..d2fb3541b 100644 --- a/src/QLNet/Methods/Finitedifferences/TRBDF2.cs +++ b/src/QLNet/Methods/Finitedifferences/TRBDF2.cs @@ -77,7 +77,7 @@ public Trbdf2(Operator L, List> bcs) alpha_ = 2.0 - Math.Sqrt(2.0); } - public void step(ref object a, double t) + public void step(ref object a, double t, double theta = 1.0) { int i; Vector aInit = new Vector((a as Vector).size());