Ode/example

examples for class banach and classes in ode

examples/num0/wurzelbanach.cc

+#include <iostream>
+#include <vector>
+#include "hdnum.hh"
+
+
+template<class N>
+class WurzelProblem
+{
+public:
+  /** \brief export size_type */
+  typedef std::size_t size_type;
+
+  /** \brief export number_type */
+  typedef N number_type;
+
+  //! constructor stores parameter lambda
+  WurzelProblem (number_type a_)
+    : a(a_)
+  {}
+
+  //! return number of componentes for the model
+  std::size_t size () const
+  {
+    return 1;
+  }
+
+  //! model evaluation
+  void F (const hdnum::Vector<N>& x, hdnum::Vector<N>& result) const
+  {
+    result[0] = x[0]*x[0] - a;
+  }
+
+  //! jacobian evaluation needed for implicit solvers
+  void F_x (const hdnum::Vector<N>& x, hdnum::DenseMatrix<N>& result) const
+  {
+    result[0][0] = number_type(2.0)*x[0];
+  }
+
+private:
+  number_type a;
+};
+
+
+
+int main ()
+{
+  //typedef double Number;                 // numbertype
+  typedef mpf_class Number;
+  mpf_set_default_prec(1024);
+
+  typedef WurzelProblem<Number> Problem; // problem type
+  Problem problem(2.0);                  // an instance of the problem
+
+  hdnum::Banach banach;                         // an Object of class Banach
+  banach.set_maxit(20);                  // set parameters
+  banach.set_verbosity(2);
+  banach.set_reduction(1e-100);
+  banach.set_abslimit(1e-100);
+  banach.set_linesearchsteps(3);
+  banach.set_sigma(0.1);
+
+  hdnum::Vector<Number> u(problem.size());      // vector for storing solvations
+  u[0] = 2.0;                            // value for starting
+  banach.solve(problem,u);               // compute solvation
+  std::cout << "Ergebnis: " << u[0] << std::endl;
+
+  return 0;
+}

examples/num1/beispielDGL.cc

+#include <iostream>
+#include <vector>
+#include "hdnum.hh"
+
+#include "beispielDGL.hh"
+#include "../../src/ode.hh"
+
+int main ()
+{
+  typedef double Number;               // define a number type
+
+  typedef ModelProblem<Number> Model;  // Model type
+  Model model;                   // instantiate model
+
+  typedef hdnum::EE<Model> Solver; // Solver type
+  Solver solver(model);                // instantiate solver
+  solver.set_dt(0.02);                  // set initial time step
+
+  hdnum::Vector<Number> times;           // store time values here
+  hdnum::Vector<hdnum::Vector<Number>> states; // store states here
+  times.push_back(solver.get_time());  // initial time
+  states.push_back(solver.get_state()); // initial state
+	
+
+  while (solver.get_time()<5.0-1e-6) // the time loop
+    {
+      solver.step();                  // advance model by one time step
+      times.push_back(solver.get_time()); // save time
+      states.push_back(solver.get_state()); // and state
+    }
+
+  gnuplot("beispielmp2-ee-0.02.dat",times,states); // output model result
+
+  return 0;
+}

examples/num1/beispielDGL.hh

+/** @brief Example class for a differential equation model
+
+    The model is
+
+		u'(t) = (5 -2// -2 5)*u(t)
+		u(0) = (1 3)
+    \tparam T a type representing time values
+    \tparam N a type representing states and f-values
+*/
+
+template<class T, class N=T>
+class ModelProblem
+{
+public:
+  /** \brief export size_type */
+  typedef std::size_t size_type;
+
+  /** \brief export time_type */
+  typedef T time_type;
+
+  /** \brief export number_type */
+  typedef N number_type;
+
+  //! constructor
+  ModelProblem ()
+  {}
+
+  //! return number of components for the model
+  std::size_t size () const
+  {
+    return 2; // as the dimension of the differential equation ist two
+  }
+
+  //! set initial state including time value
+  void initialize (T& t0, hdnum::Vector<N>& x0) const
+  {
+    t0 = 0;
+    x0[0] = 1.0;
+		x0[1] = 3.0;
+  }
+
+  //! model evaluation
+  void f (const T& t, const hdnum::Vector<N>& x, hdnum::Vector<N>& result) const
+  {
+    result[0] = 5*x[0] - 2*x[1];
+		result[1] = -2*x[0]+5*x[1];
+  }
+
+  //! jacobian evaluation needed for implicit solvers
+  void f_x (const T& t, const hdnum::Vector<N>& x, hdnum::DenseMatrix<N>& result) const
+  {
+		result[0][0] = 5.0;
+		result[0][1] = -2.0;
+		result[1][0] = -2.0;
+		result[1][1] = 5.0;
+  }
+
+private:
+  //
+};